植生生育状況解析方法、プログラムおよび植生生育状況解析装置
【課題】人工衛星等の飛翔体から得られるレーダ画像を利用して、多時期の植生の生育状況を正確に解析できるようにする。
【解決手段】飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得し、取得された複数のレーダ画像を地図データベースに保存する。次に、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行った後、前記レーダ画像の特定領域の後方散乱係数を抽出する。さらに、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する。そして、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する。
【解決手段】飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得し、取得された複数のレーダ画像を地図データベースに保存する。次に、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行った後、前記レーダ画像の特定領域の後方散乱係数を抽出する。さらに、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する。そして、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、人工衛星等の飛翔体に搭載され、撮影対象地域の地表面を広範囲にわたり撮影して地表の状況を取得するレーダ装置により得られるレーダ画像データを利用して、植生の生育状況を解析するのに好適な植生生育状況解析方法、プログラムおよび植生生育状況解析装置に係わる。
【背景技術】
【０００２】
従来、農作物の生産高の推定・予測（生育状況解析）は、人工衛星や航空機から取得した光学リモートセンシングデータをもとに、アルゴリズムが構築され、実用化に至ってきた。しかし、日本の場合は、四季や雨期があり年間を通して雨や曇りの日が多いため、天候に左右される光学リモートセンシングでは安定した観測が難しい。
【０００３】
また、日本のように南北に続く複雑な地形で多品種・少量の生産をしている場合は、精密な農業管理が必要である。しかし、これまで圃場単位の生育状況が把握できる高精度なデータとして、光学リモートセンシングデータ以外の入手が困難であるという事情があった。
【０００４】
一方、地表面に照射した電波（マイクロ波パルス）の反射波を測定する能動型センサの例として、合成開口レーダ（SAR：Synthetic Aperture Radar）がある。合成開口レーダは、マイクロ波の特性を利用して、昼夜、天候に関係なく、広域を面的に撮影することができる。また、従来の合成開口レーダは光学リモートセンシングと比較すると解像度（精細度）に問題があったが、最新の衛星搭載型合成開口レーダには、解像度が１ｍ以下のものも存在し、高解像度の画像が得られるようになってきている。
【０００５】
農作物に対する合成開口レーダを用いた既往の研究として、下記のようなものがあった。
・光学リモートセンシングデータと多時期のＳＡＲ画像を併用した水稲作付け面積の推定
・多重波長、多重偏波のＳＡＲ画像からの水稲作付け面積
・多時期のＳＡＲ画像を用いた農作物の種別分類の研究
【０００６】
これらの研究が海外で行われていたり、多数の論文が公表されたりしているが、水稲など農作物の生育予測に直接結びつくものではなかった。
【０００７】
そこで、多時期のＳＡＲ画像を用いた水稲の生育状況について研究が行われている（例えば非特許文献１参照）。
【０００８】
本願の出願人も、多時期のＳＡＲ画像を用いた水稲生育状況に関する研究結果を発表している（非特許文献２を参照。）。非特許文献２に記載された研究により、図１〜図３に示すような解析結果を得た。
【０００９】
図１は、複数時期のレーダ画像例を示すものであり、各観測日ごとの生育と反射特性、水田の様子、水稲被覆状態、レーダ画像を表している。この例では、ＲＡＤＡＲＳＡＴのＣバンドを使用してレーダ画像を撮影している。５月２２日、６月１５日、７月９日、８月２日の４時期のレーダ画像１Ａ，１Ｂ，１Ｃ，１Ｄを作成し、これらのレーダ画像から水稲の生育状況が全体的に変化していることが把握できる。
【００１０】
図２は、水稲被覆率とレーダ後方散乱係数（レーダ装置から送信されたマイクロ波に対する対象物の散乱強度のうち後方への成分）との関係を示したものであり、横軸は水稲被覆率（％）、縦軸はレーダ後方散乱係数（ｄＢ）である。図２に示すように水稲被覆率とレーダ後方散乱係数には強い相関があり、これらの観測点に最小２乗法を適用すると、水稲被覆率とレーダ後方散乱係数との関係を表す回帰直線２が得られる。
【００１１】
図３は、各時期における圃場の区画別レーダ後方散乱係数分布を示したものである。このように、レーダ画像を利用して、時期ごとの圃場の区画別レーダ後方散乱係数分布図３Ａ，３Ｂ，３Ｃ，３Ｄを作成することにより、水稲の生育状況を区画別に把握することができる。
【００１２】
【非特許文献１】Zengyuan Li, Guoqing Sun, Mike Wooding, Yong Pang, Yanfang Dong, Erxue Chen and Bingxiang Tan, “Rice Monitoring Using Envisat Asar Data in China”, Proc. of the 2004 Envisat & ERS Symposium, Salzburg, Austria6-10 September 2004 (ESA SP-572, April 2005)
【非特許文献２】高橋一義、阿部博明、力丸厚、向井幸男、「時系列ＲＡＤＡＲＳＡＴデータを用いた水稲生育分布の把握」、日本写真測量学会 年次学術講演会 東京国際展示場（東京）、２００４年６月１７・１８日
【発明の開示】
【発明が解決しようとする課題】
【００１３】
ところで、非特許文献１には、多時期のＳＡＲ画像を用いた水稲のモニタリングについての研究内容が記載されているが、ＳＡＲ画像の後方散乱係数と植生との関係は明らかになっていない。
【００１４】
一方、非特許文献２に記載された研究により、衛星ＳＡＲデータの後方散乱係数と生育状況との相関が明らかになった。しかしながら、後方散乱係数は撮影時期における大気の状況やシステムの変動などの影響を受けるため、異なる時期に撮影したＳＡＲ画像の後方散乱係数の較正が必要であり、上記研究結果は多時期の水稲生育状況を正確に表してはいなかった。
【００１５】
本発明は、このような状況に鑑みてなされたものであり、人工衛星等の飛翔体から得られるレーダ画像を利用して、多時期の植生の生育状況を正確に解析できるようにするものである。
【課題を解決するための手段】
【００１６】
本発明の一側面の植生生育状況解析方法およびプログラムは、飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得し、取得された複数のレーダ画像を地図データベースに保存する。次に、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行った後、前記レーダ画像の特定領域の後方散乱係数を抽出する。さらに、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する。そして、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する、これらの処理または機能を含むことを特徴とする。
【００１７】
本発明の一側面の植生生育状況解析方法装置は、飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得部と、前記取得部で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持部と、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理部と、前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理部と、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正部と、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正部により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出部とを備えることを特徴とする。
【００１８】
本発明の一側面においては、所定時期のレーダ画像の後方散乱係数を基準として各時期のレーダ画像の後方散乱係数が較正される。そして、較正後の後方散乱係数を基に各時期のレーダ画像が補正されるので、各時期における農作物等の植物の正確な生育状況を示す生育値を算出することができる。
【発明の効果】
【００１９】
以上のように、本発明によれば、飛翔体から得られる多時期のレーダ画像を利用して、多時期の植生の生育状況を正確に解析できる。
【発明を実施するための最良の形態】
【００２０】
以下、本発明を実施するための最良の形態の例について、添付図面を参照しながら説明する。説明は下記項目に従って順に行う。
１．システムの概要
２．植生生育状況解析装置
３．水稲の生育予測処理
４．測定結果等
５．その他
【００２１】
＜１．システムの概要＞
【００２２】
図４は、本発明の植生生育状況解析装置が適用されるシステムの概要を示すものである。
本実施の形態においては、人工衛星にレーダ装置の一例として合成開口レーダ（ＳＡＲ）を搭載し、その合成開口レーダにより（データ解析センタ１２からの指示に基づいて）地表面を撮影する例としてある。以下、人工衛星に搭載した合成開口レーダを「衛星ＳＡＲ」と称す。
【００２３】
人工衛星に搭載された衛星ＳＡＲ１１は、所定の軌道を周回しながら定期的に、もしくはデータ解析センタ１２からの指示に従って随時地表面の圃場１３を撮影し、撮影データ（レーダ画像データ）をデータ解析センタ１２へ送信する。この例では、農作物として水稲が圃場１３に栽培されているとする。
【００２４】
データ解析センタ１２は、生育カレンダに沿って撮影計画を行い、その撮影計画に基づいた撮影指示を含む無線信号を衛星ＳＡＲ１１に送信し、また、衛星ＳＡＲ１１が撮影したレーダ画像データを、アンテナを介して受信する。これに対して後述する植生生育状況解析装置２０で合成開口処理および所定の補正（較正）を行い、補正された多時期の再生画像やレーダ後方散乱係数（単に「後方散乱係数」ともいう。）の分布図等を生成する。
【００２５】
リアルタイムで水稲の生育評価を行うため、本発明では、最初もしくは所定の時期に取得した画像を基準とした、時系列の画像の後方散乱係数値の較正を行う。この後方散乱係数値の較正を行うことにより、撮影時期による変動の影響を考慮した、より実用的な水稲生育予測システムを実現することができる。
【００２６】
時系列の画像の較正方法は、以下の３つの方法が考えられる。
（１） 対象エリアに基準となる反射板（例えばコーナリフレクタ）を設置して、撮影したＳＡＲ画像から反射板の後方散乱係数値を抽出し、これを基準値とする方法
（２） 撮影したＳＡＲ画像の中で、低変動領域を抽出し、その領域の後方散乱係数値を用いて較正するための回帰線（較正線）を求め、これに後方散乱係数分布を当てはめる方法
（３） 上記（１），（２）の方法を併用する方法
【００２７】
＜２．植生生育状況解析装置＞
［植生生育状況解析装置の構成］
図５は、本発明の植生生育状況解析装置の一実施の形態に係る機能ブロック図を示したものである。植生生育状況解析装置は、上記時系列画像の後方散乱係数値の較正を行い、植生の生育状況を予測するためのデータを生成するものである。この植生生育状況解析装置は、データ解析センタ１２に設置されることを想定しているが、衛星ＳＡＲ１１とともに人工衛星に搭載してもよいし、それぞれに分散してもよい。また、植生生育状況解析装置の各機能は、ハードウェアまたはソフトウェアのいずれでも実現することができる。あるいはそれらの組み合わせにより実現してもよい。
【００２８】
本実施の形態の植生生育状況解析装置２０は主に、Ｉ／Ｆ部２１、撮影決定部２２、レーダ画像保持部２４、位置合わせ処理部２５、基準反射板の後方散乱係数処理部２７、後方散乱係数の判別処理部２８、地物の後方散乱係数処理部２９を備える。また、後方散乱係数保持部３０、較正用後方散乱係数算出部３１、較正線算出部３２、較正線による後方散乱係数較正部３３、生育値算出部３４、圃場集約算出部３６、出力部３８を備える。さらに、生育カレンダデータベース（以下、「生育カレンダＤＢ」という。）２３Ａ、撮影計画データベース（以下、「撮影計画ＤＢ」という。）２３Ｂ、地図データベース（以下、「地図ＤＢ」という。）２６、生育モデルデータベース（以下、「生育モデルＤＢ」という。）３５、圃場地図データベース（以下、「圃場地図ＤＢ」という。）３７を備える。
【００２９】
Ｉ／Ｆ部２１は、取得部の一例であり、衛星ＳＡＲ１１で撮影した撮影データを取り込むいわゆるインタフェースとしての機能を備える。Ｉ／Ｆ部２１は、通信部としての機能も備え、撮影決定部２２からの指示を衛星ＳＡＲ１１に送信し、衛星ＳＡＲ１１に撮影を実行させることができる。撮影データは、生育カレンダＤＢ２３Ａに登録された水稲の生育カレンダや、撮影計画ＤＢ２３Ｂに登録された衛星の撮影時期（周期）をもとに、計画的に撮影されるレーダ画像のデータである。
【００３０】
生育カレンダは、例えば水稲の場合であれば、５月の田植え期から始まって分けつ期、最高分けつ期、幼穂形成期、出穂期、登熟期、そして９月の成熟期までの栽培暦である。また、現在打ち上げられている合成開口レーダ搭載の人工衛星には、例えば、観測周期が１１日のものがある。
【００３１】
撮影決定部２２は、これらの生育カレンダと衛星の撮影周期を考慮してレーダ画像をＩ／Ｆ部２１に取り込む時期を決定する。決定した撮影予定日（データ取得日）の一例を以下に示す。
［観測日］ ［起算日］ ［段階］
１）５月 ７日 ０日 田植期
２）５月３１日 ２４日 分けつ期
３）６月２４日 ３８日 最高分けつ期
４）７月１８日 ７２日 幼穂形成期
５）８月１１日 ９６日 出穂期
６）９月 ４日 １２０日 登熟期
【００３２】
なお、撮影決定部２２は、水稲の生育の評価を行う対象面積により、撮影画像の画素サイズを設定するようにしてもよい。例えば、詳細な評価を行う場合は、地上分解能を高く（例えば１ｍ）することが好ましい。
【００３３】
レーダ画像保持部２４は、所定の時期に撮影された時系列のレーダ画像を、地図ＤＢ２６に保存する機能である。地図ＤＢ２６は、例えばハードディスクやフラッシュメモリ等の不揮発性の記憶手段に構築される。
【００３４】
位置合わせ処理部２５は、地図ＤＢ２６に保存されている地図データをもとに、レーダ画像を地図に重ね合わせることができるように、幾何（位置）的な処理を行う機能である。また、位置合わせ後のレーダ画像を、較正線による後方散乱係数較正部３３、基準反射板の後方散乱係数処理部２７および地物の後方散乱係数処理部２９へ供給する。
【００３５】
基準反射板の後方散乱係数処理部２７は、レーダ画像に写っている基準反射板の後方散乱係数を処理（抽出）する機能である。この基準反射板としては、例えば「加美川雄太、大内和夫、“陸域観測技術衛星「だいち」搭載合成開口レーダＰＡＬＳＡＲの校正・検証実験”、平成１８年度卒業論文、高知工科大学工学部」（文献１）に記載されたコーナリフレクタなど、種々の反射板を適用できる。
【００３６】
観測対象エリアに反射板を用いる場合、２時期の反射板の後方散乱係数値をそれぞれ抽出して、基準の画像の後方散乱係数値をもとにした較正を行うための較正線を求め、次の時期や他の時期に撮影した画像を較正する。文献１にも記載されているように、反射板は、地表面に照射するマイクロ波のバンド（波長）や人工衛星の軌道（撮影方向）によって、適宜、適切な面積や形状、個数、方向を調整して設置することが望ましい。
【００３７】
［基準反射板（コーナーリフレクタ）］
図６は、任意の観測対象エリアに、反射板としてコーナリフレクタを設置した場合のレーダ画像と平均後方散乱係数の計測結果を示したものである。図の左側がレーダ画像、右側が後方散乱係数のヒストグラム（頻度分布）を示すグラフである。
【００３８】
図６に示す観測日５月７日のレーダ画像４１に写るコーナリフレクタ４１ａに相当する画素、すなわち後方散乱係数が０ｄＢ以上の画素数は４、平均後方散乱係数が５．１３ｄＢである。同様に、図７に示す観測日８月１１日のレーダ画像４５に写るコーナリフレクタ４５ａは画素数が１４、平均後方散乱係数が５．６１ｄＢである。
【００３９】
このように、時期によってコーナリフレクタの画素数および平均後方散乱係数が異なるので、この情報を後方散乱係数の較正に反映させる。
【００４０】
後方散乱係数の判別処理部２８は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数の利用の可否を判定する機能であり、後方散乱係数が飽和状況の場合、すなわち、レーダ画像の中のコーナリフレクタが設置されている場所以外に、コーナリフレクタにおける後方散乱係数と同じ値の画素が存在する場合は利用しない。後方散乱係数が飽和している場合、正確な後方散乱係数の情報が得られないため精度の良い較正を行うことができないからである。
【００４１】
地物の後方散乱係数処理部２９は、レーダ画像に写っている地物（人工構造物、池、グラウンドなど）の後方散乱係数を抽出する機能である。ここで処理された後方散乱係数をもとにレーダ画像中の低変動領域を決定する。あるいは過去のデータから低変動領域の特徴を示す地物を予め特定しておき、これを利用してもよい。
【００４２】
［複数地点・多時期のレーダ後方散乱係数］
図８、図９は、多地点・多時期で観測したレーダ後方散乱係数及びその経時変化の例を示すものである。なお、図８、図９は、観測結果に平滑処理として後述するメディアンフィルタを適用したものである。
【００４３】
図８では、特徴領域としてショッピングセンタ、仮設住宅、調整池、グラウンド、駐車場（公園）、池（公園）、並びにコーナリフレクタの７地点、また撮影日として５月７日、５月３１日、６月２４日、７月１８日、８月１１日、９月４日において得られたデータを示している。各データは、特徴領域を構成する各画素の平均値をとったものである。
【００４４】
図９は、図８に示したデータを特徴領域のレーダ後方散乱係数の経時変化としてグラフに表したものである。例えば、変化が少なくレーダ後方散乱係数の大きいショッピングセンタや仮設住宅が、レーダ画像の後方散乱係数を較正する際の基準となる低変動領域の候補として挙げられる。また、変化が少なくレーダ後方散乱係数の小さい調整池を低変動領域の候補として挙げることも可能である。これらの低変動領域が写った、最初または所定の時期のレーダ画像が較正時の基準画像として利用される。
【００４５】
図５に示した植生生育状況解析装置の説明に戻る。
後方散乱係数保持部３０は、レーダ画像から抽出された基準反射板および地物の後方散乱係数を保存する機能である。
【００４６】
較正用後方散乱係数算出部３１は、後方散乱係数保持部３０に保存された後方散乱係数に対して、平滑処理および／または局所最小値処理を行い、較正用後方散乱係数を求める機能であり、後方散乱係数平滑処理部３１Ａ、後方散乱係数局所最小値算出部３１Ｂを備える。
【００４７】
後方散乱係数平滑処理部３１Ａは、レーダ画像に対して所定の平滑処理用フィルタ、例えばメディアンフィルタをかけて、スペックルノイズ等のノイズ成分を低減する機能を備える。メディアンフィルタは、ｋ×ｋの局所領域における値を順番に並べ、その真ん中の値を領域中央の画素の出力値とする処理を行うフィルタであり、その処理は式（１）で表すことができる。
y[m，n] = median{x[i,j],(i,j)∈ω}・・・・（１）
ここで、y[m,n]は対象画素の値、x[i,j]は対象画素を中心とする近傍画素の値
【００４８】
後方散乱係数局所最小値算出部３１Ｂは、レーダ画像に対して所定の局所最小値処理用フィルタ（以下、「最小値フィルタ」という。）をかけて、所定領域の最小値を抽出する機能を備える。最小値フィルタは、ｋ×ｋの局所領域における値の最小値を、領域中央の画素の出力値とする処理を行うフィルタであり、その処理は式（２）で表すことができる。
y[m，n] = min{x[i,j],(i,j)∈ω}・・・・（２）
ここで、y[m,n]は対象画素の値、x[i,j]は対象画素を中心とする近傍画素の値
【００４９】
上述のとおり、後方散乱係数平滑処理部３１Ａおよび後方散乱係数局所最小値算出部３１Ｂは、それぞれレーダ画像に発生するノイズ成分の低減および所定領域の最小値を抽出するものである。これらの機能ブロックがない場合、レーダ画像の画質が劣化し、較正の精度が低下することが想定されるが、取得したレーダ画像が要求仕様を満たしてさえいれば必ずしも必要ではなく、また何れか一方のみを用いたり、あるいは処理順序を逆にしたりすることも可能である。
【００５０】
較正線算出部３２は、例えば最初に撮影した画像（第１画像）に対して、順次撮影した画像（第２画像、第３画像・・・）を較正するために、第１画像の較正用後方散乱係数をもとに、第２画像の較正用後方散乱係数との関係式（較正線）を算出する機能である。この較正線算出部３２は、下記較正線による後方散乱係数較正部３３とともに、特許請求の範囲に記載した較正手段を構成する。
【００５１】
較正線による後方散乱係数較正部３３は、後方散乱係数較正手段の一例であり、較正線算出部３２で算出された較正線を用いて、基準画像以外の画像、例えば第２画像、第３画像・・・の全画素の後方散乱係数を処理する機能である。
【００５２】
生育値算出部３４は、生育モデルＤＢ３５に保存された生育モデル（レーダ画像の後方散乱係数と水稲被覆率の相関関係）の情報をもとに、較正線による後方散乱係数較正部３３で求められた後方散乱係数から水稲被覆率を算出する機能である。
【００５３】
圃場集約算出部３６は、圃場地図ＤＢ３７に保存された圃場地図データをもとに、一筆単位の水稲被覆率を算出する機能であり、圃場代表値算出部３６Ａ、圃場代表値差分算出部３６Ｂ、生育判定部３６Ｃを備える。
【００５４】
圃場代表値算出部３６Ａは、１つのレーダ画像を構成する各画素の水稲被覆率から一筆単位の平均値、もしくは、出現頻度が最も高い値（最頻値）を求める機能である。
【００５５】
圃場代表値差分算出部３６Ｂは、基準画像（例えば第１画像）の圃場代表値と参照画像（例えば第２画像）の圃場代表値の差分を計算する機能である。
【００５６】
生育判定部３６Ｃは、圃場代表値および圃場代表値差分の値から生育判定（生育の早い、遅い）を行う機能である。
【００５７】
出力部３８は、水稲被覆率の圃場代表値、圃場代表値差分の値、および、生育判定の結果を、圃場単位の一覧表、もしくは、ランキングした地図（ランキングマップ）などに出力する機能である。
【００５８】
［後方散乱係数の補正の効果］
図１０、図１１は、後方散乱係数の補正の効果を説明するための図である。図１０は、８月１１日におけるあるショッピングセンタ駐車場の場合を示し、図１１は、８月１１日におけるある調整池の場合を示している。
画像５１Ｂは、オリジナルの画像（補正前の画像）５１Ａにメディアンフィルタを実施した後の画像、画像５１Ｃは、メディアンフィルタを実施した後の画像５１Ｂに最小値フィルタを実施した画像である。オリジナルの画像５１Ａはノイズが多いが、メディアンフィルタ処理および最小値フィルタ処理後は、ノイズが除去されるとともに領域中央の画素が周辺の最小値に置き換えられている。この結果、ショッピングセンタ駐車場のように後方散乱係数自体が大きくその変動が少ない領域において、較正用の後方散乱係数を安定的に求めることができる（矢印参照）。
【００５９】
図１１は調整池の場合の例であり、画像５２Ｂは、オリジナルの画像（補正前の画像）５２Ａにメディアンフィルタを実施した後の画像、画像５２Ｃは、メディアンフィルタを実施した後の画像５２Ｂに最小値フィルタを実施した画像である。調整池の場合、図１０に示したショッピングセンタ駐車場の場合と比較して、さらに後方散乱係数自体が小さくその変動が少ない調整池のような領域においても較正用の後方散乱係数を安定的に求めることができる（矢印参照）。
【００６０】
＜３．水稲の生育予測処理＞
［植生生育状況解析装置の動作］
次に、図１２に示すフローチャートを参照して、植生生育状況解析装置２０による水稲の生育状況解析処理を説明する。
図５において、植生生育状況解析装置２０の撮影決定部２２により、適切な撮影時期、場所および地上分解能等を決定し、これらの情報をＩ／Ｆ部２１を介して衛星ＳＡＲ１１へ送信する。Ｉ／Ｆ部２１において、衛星ＳＡＲ１１で撮影されたＳＡＲ画像（レーダ画像）が人工衛星から受信される（ステップＳ１）。受信した各時期のレーダ画像は、レーダ画像保持部２４により地図ＤＢ２６に保存される。
【００６１】
続いて、位置合わせ処理部２５により、地図ＤＢ２６に保存されている最初に撮影したレーダ画像を基準画像として、次または他の時期に撮影したレーダ画像（同一対象エリア）が重なるように位置合わせが行われる（ステップＳ２）。位置合わせ処理では、例えば緯度経度が明確な地上基準点（GCP：Ground Control Point）や道路交差点、建物など画像に表れている点を基準として、地図に重ね合わせることができるように幾何学的な補正が行われる。
【００６２】
位置合わせが終了後、それぞれの場合に応じて、多時期のレーダ画像の後方散乱係数処理および較正が行われる（ステップＳ３）。
【００６３】
第１に、コーナリフレクタ等の基準反射板を利用する場合を説明する。
基準反射板の後方散乱係数処理部２７は、複数時期のレーダ画像中に写るコーナリフレクタの後方散乱係数値をそれぞれ抽出して、基準画像の後方散乱係数値をもとにした較正線を求め、当該較正線に基づいて次の時期や他の時期に撮影したレーダ画像の後方散乱係数値を処理する。このとき、後方散乱係数の判別処理部２８は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数の利用の可否を判定し、後方散乱係数の値が飽和状況の場合は利用しない。一方、後方散乱係数の値が飽和状況でない場合は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数値を、レーダ画像の撮影時期と対応づけて後方散乱係数保持部３０へ供給する。
【００６４】
そして、後方散乱係数保持部３０に保存された後方散乱係数の情報は、較正用後方散乱係数算出部３１の後方散乱係数平滑処理部３２Ａおよび後方散乱係数局所最小値算出部３２Ｂによってそれぞれフィルタ処理されてノイズ成分が低減され、レーダ画像上の代表的な領域内の後方散乱係数値が求められた後、較正線算出部３２に供給される。較正線算出部３２では、基準画像の低変動領域（内の複数の画素）の後方散乱係数と、その基準画像の低変動領域に対応する他のレーダ画像の特定領域（内の複数の画素）の後方散乱係数から較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６５】
第２に、レーダ画像中の低変動領域の後方散乱係数を用いる場合を説明する。
この場合、較正の基準となる後方散乱係数値は取得したレーダ画像の中から抽出する。地物の後方散乱係数処理部２９は、レーダ画像に写っている地物の後方散乱係数値を計算し、後方散乱係数が大きくその変動が少ない画像領域（例えば、人工構造物）、および、後方散乱係数が低くその変動が少ない画像領域（例えば、池、土のグラウンド）を抽出する。これら後方散乱係数の変動が少ない画像領域、すなわち低変動領域の後方散乱係数値を、レーダ画像の撮影時期と対応付けて後方散乱係数保持部３０へ供給する。
【００６６】
そして、後方散乱係数保持部３０に保存された後方散乱係数の情報は、第１の場合と同様に後方散乱係数平滑処理部３１Ａおよび後方散乱係数局所最小値算出部３１Ｂにてフィルタ処理された後、較正線算出部３２に供給される。較正線算出部３２では、レーダ画像に写る基準反射板の後方散乱係数が求められ、基準画像に写る基準反射板（の複数の画素）の後方散乱係数と、基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板（の複数の画素）の後方散乱係数から較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６７】
第３に、第１項目のコーナリフレクタ等の基準反射板および第２項目のレーダ画像中の低変動領域を併用する場合を説明する。
この場合、較正線算出部３２では、基準画像において後方散乱係数の変動が少ない画像領域の状況に応じて、第２項目におけるレーダ画像上の低変動領域の後方散乱係数のデータと、第１項目におけるコーナリフレクタ等の基準反射板を用いて取得したデータとを選択的に用いて較正線が算出される。すなわち、基準画像から低変動領域を抽出し、これが基準反射板である場合にはその領域の後方散乱係数および基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板の後方散乱係数を用いて、低変動領域が人工構造物等である場合にはその領域の後方散乱係数および基準画像の人工構造物等に対応する他のレーダ画像の特定領域の後方散乱係数を用いて較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６８】
そして、生育値算出部３４により生育モデルを基にした水稲の生育解析が行われる（ステップＳ４）。水稲の生育状況の解析は、予め把握している後方散乱係数値と生育情報のモデルを基に、較正処理済みのレーダ画像から計算する。すなわち、生育モデルＤＢ３５に保存された、レーダ画像の後方散乱係数と植生の生育値の相関関係の情報に基づいて、較正されたレーダ画像の後方散乱係数から当該レーダ画像の各画素に写っている植生の生育値が算出される。
【００６９】
水稲の生育解析が終了後、圃場集約算出部３６により解析結果の集計処理が行われる（ステップＳ５）。例えば、圃場代表値算出部３６Ａにより、上記解析結果から、１つのレーダ画像を構成する各画素の水稲被覆率から一筆単位の平均値、もしくは、最頻値が求められた場合、圃場ごとの生育状況が把握できる。
また、圃場代表値差分算出部３６Ｂにより、基準画像（例えば第１画像）の圃場代表値と参照画像（例えば第２画像）の圃場代表値の差分が求められた場合、２つの時期の間で圃場ごとの生育の違いを把握できる。
さらに、生育判定部３６Ｃにより、圃場代表値および圃場代表値差分の値に基づいて生育判定が行われた場合、圃場間での生育の早いか遅いかの違いを把握できる。
【００７０】
そして、出力部３８により水稲生育状況分布の作成が行われる（ステップＳ６）。すなわち、水稲被覆率の圃場代表値、圃場代表値差分の値、および、生育判定の結果を、圃場単位の一覧表、もしくは、ランキングした地図（ランキングマップ）などに出力する処理が行われる。これらの水稲生育状況分布を表した地図は、水稲生育状況分布図と呼ばれる。
【００７１】
水稲生育状況分布図は、ピクセル（画素）単位、圃場単位、圃場の中を細分化した単位など、営農指導を行う単位ごとに情報の集計処理を行うことにより、個々のユーザの実情に合わせた営農指導を行える。
【００７２】
ピクセル（画素）単位から圃場単位に集計するためには、圃場の形状のデータ、例えば、地理情報データを利用する。対象エリアの地理情報データは、地図ＤＢ２６等に保存しておく。
【００７３】
＜４．測定結果等＞
【００７４】
図１３は、所定の時期（この例では５月７日）のレーダ画像を基に、８月１１日に取得したデータを較正（図１２のステップＳ３）した結果を示すものである。図１３Ａは対象時期の較正後のレーダ画像、図１３Ｂは図１３Ａに示したレーダ画像の較正に用いた較正線を示す図である。なお、本例では後方散乱係数が最小となる撮影時期の最初（５月７日）を基準にグラフを作成したが、後述するように後方散乱係数が最大となる時期である８月１１日を基準として標準化することも可能である。
【００７５】
図１３Ｂは、８月１１日のデータを較正するために用いる、５月７日のレーダ画像を基にした較正線である。縦軸が５月７日の後方散乱係数であり、横軸が８月１１日における後方散乱係数である。この例では、較正線（回帰直線）は下記式で表される。ただし、ｙは基準画像（５月７日）の後方散乱係数、ｘは較正する画像（８月１１日）の後方散乱係数、Ｒ２は決定係数である。なお、図１３Ｂに示したデータは、最小値フィルタを適用後のデータである。
ｙ＝０．９８４ｘ−０．１３７
Ｒ２＝０．９７２
【００７６】
同様にして、所定時期（この例では５月７日）のレーダ画像の後方散乱係数を基準として作成した各時期の較正線により、それぞれの時期に取得したレーダ画像の後方散乱係数が較正される。この結果、レーダ画像の後方散乱係数が較正されるので、各時期における水稲（農作物等の植物）の正確な生育状況を示す生育値を算出することができる。
【００７７】
次に、圃場のレーダ後方散乱係数分布図について説明する。
図１４は、ある圃場（水稲作付圃場）についてのレーダ後方散乱係数の分布を示したものであり、複数の圃場を含む対象エリア全体のレーダ後方散乱係数分布図の要部の拡大図である。図１４において、拡大されたレーダ後方散乱係数分布図には圃場７１〜７３が含まれている。このように、レーダ画像の任意の領域に対応する農作物の生育状況について、画素単位で精細に確認することができる。
【００７８】
図１５は、図１４に示した水稲作付圃場のレーダ後方散乱係数分布図に対応する平均後方散乱係数の経時変化を示すグラフである。平均後方散乱係数とは、レーダ画像内のある領域を構成する複数の画素から得られた後方散乱係数の平均値を算出したものである。圃場７１〜７３の水稲は、田植え後から徐々に成長した結果として、後方散乱係数が大きくなり、８月１１日を極大点として後方散乱係数が小さくなっている。これは、出穂および稲の垂下によって後方散乱が減少した結果と推測される。この図１５のグラフからは、例えば圃場７３の水稲が早熟で、圃場７１の水稲は晩熟であることが把握できる。
【００７９】
本発明の観測の対象は、水稲だけでなくトウモロコシや大麦等の生育段階に伴う形態変化に応じて後方散乱係数が変化する農作物、植物が含まれ、地表面の植生について基準画像を用いて補正した複数の時期のレーダ画像や、補正後のレーダ画像から作成したレーダ後方散乱係数分布図（植生の生育状況分布図）を観察することにより、植生の生育状況をリアルタイムかつ正確に把握できる。
【００８０】
上述した実施の形態によれば、所定時期のレーダ画像の後方散乱係数を基準として各時期のレーダ画像の較正線を算出し、算出した各時期の較正線により、それぞれの時期に取得したレーダ画像の後方散乱係数が較正される。そして、較正後の後方散乱係数を基に各時期のレーダ画像が補正されるので、各時期における農作物等の植物の正確な生育状況を知ることができる。
【００８１】
また、人工衛星等の飛翔体からマイクロ波を用いて撮像したレーダ画像は、全天候型であるために、５月（田植え）から８月（出穂）の間の雨期でも農作物等の植物の生育状況を把握することができる。
【００８２】
また、多時期にわたる正確でリアルタイムな生育状況の把握により、農作物等の植物の適時かつ適切な栽培管理、追肥の投入時期の決定等が可能となる。
【００８３】
また、多時期にわたる正確でリアルタイムな生育状況の把握により、病害虫などの発生の予防、生育障害、災害への対処が可能となる。
【００８４】
さらに、官公庁等で実施している現地調査による水稲作況調査が観測手段（リモートセンシング）で可能となる。
【００８５】
＜５．その他＞
【００８６】
なお、上述した植生生育状況解析装置で行われる一連の処理は、ハードウェアにより実行させることもできるし、ソフトウェアにより実行させることもできる。また、これらの処理を実行する機能はハードウェアとソフトウェアの組み合わせによっても実現できることは言うまでもない。一連の処理をソフトウェアにより実行させる場合には、そのソフトウェアを構成するプログラムが専用のハードウェアに組み込まれているコンピュータ、または、各種のプログラムをインストールすることで各種の機能を実行することが可能な、例えば汎用のコンピュータなどに、プログラム記録媒体からインストールされる。
【００８７】
図１６は、上述した一連の処理をプログラムにより実行するコンピュータの構成例を示すブロック図である。このコンピュータ１０１は、例えば一連の処理を実行するために高性能化した専用コンピュータの他、一定の性能を備えるパーソナルコンピュータなどでもよい。
【００８８】
コンピュータ１０１のＣＰＵ（Central Processing Unit）１１１は、ＲＯＭ（Read Only Memory）１１２、または記録部１１８に記録されているプログラムに従って、上記一連の処理の他、各種の処理を実行する。ＲＡＭ（Random Access Memory）１１３には、ＣＰＵ１１１が実行するプログラムやデータなどが適宜記憶される。これらのＣＰＵ１１１、ＲＯＭ１１２、およびＲＡＭ１１３は、バス１１４により相互に接続されている。
【００８９】
ＣＰＵ１１１にはまた、バス１１４を介して入出力インタフェース１１５が接続されている。入出力インタフェース１１５には、キーボード、マウス、マイクロホンなどよりなる入力部１１６、ディスプレイ、スピーカなどよりなる出力部１１７が接続されている。ＣＰＵ１１１は、入力部１１６から入力される指令に対応して各種の処理を実行する。そして、ＣＰＵ１１１は、処理の結果を出力部１１７に出力する。
【００９０】
入出力インタフェース１１５に接続されている記録部１１８は、例えばハードディスクからなり、ＣＰＵ１１１が実行するプログラムや各種のデータを記録する。
【００９１】
通信部１１９は、インターネットやローカルエリアネットワークなどのネットワークを介して外部の装置と通信する。また通信部１１９を介してプログラムを取得し、記録部１１８に記録してもよい。
【００９２】
入出力インタフェース１１５に接続されているドライブ１２０は、磁気ディスク、光ディスク、光磁気ディスク、あるいは半導体メモリなどのリムーバブルメディア１３１が装着されたとき、それらを駆動し、そこに記録されているプログラムやデータなどを取得する。取得されたプログラムやデータは、必要に応じて記録部１１８に転送され、記録される。
【００９３】
コンピュータにインストールされ、コンピュータによって実行可能な状態とされるプログラムを格納するプログラム記録媒体は、図３０に示すように、リムーバブルメディア１３１によりパッケージメディアとして提供される。リムーバブルメディア１３１としては、磁気ディスク（フレキシブルディスクを含む）、光ディスク（ＣＤ−ＲＯＭ（Compact Disc - Read Only Memory），ＤＶＤ（Digital Versatile Disc），光磁気ディスクを含む）、もしくは半導体メモリなどを適用することができる。あるいは、プログラム記録媒体は、プログラムが一時的もしくは永続的に格納（記録）されるＲＯＭ１１２や、記録部１１８を構成するハードディスクなどにより構成される。
【００９４】
このプログラム記録媒体へのプログラムの格納は、必要に応じてルータ、モデムなどのインタフェースである通信部１１９を介して、ローカルエリアネットワーク(LAN：Local Area Network)、インターネット、デジタル衛星放送といった、有線または無線の通信媒体を利用して行われる。
【００９５】
なお、本明細書において、プログラム記録媒体に格納されるプログラムを記述する処理ステップは、記載された順序に沿って時系列的に行われる処理はもちろん、必ずしも時系列的に処理されなくとも、並列的あるいは個別に実行される処理（例えば、並列処理あるいはオブジェクトによる処理）をも含むものである。
【００９６】
また、プログラムは、一つのコンピュータにより処理されるものであってもよいし、複数のコンピュータによって分散処理されるものであってもよい。さらに、プログラムは、遠方のコンピュータに転送されて実行されるものであってもよい。
【００９７】
以上に述べた実施の形態は、本発明を実施するための好適な形態の具体例であるから、技術的に好ましい種々の限定が付されている。ただし、本発明は、以上の実施の形態の説明において特に本発明を限定する旨の記載がない限り、これらの実施の形態に限られるものではない。したがって、例えば、以上の説明で挙げた処理時間、処理順序および各パラメータの数値的条件等は好適例に過ぎず、また、説明に用いた各図における寸法、形状および配置関係等も実施の形態の一例を示す概略的なものである。したがって、本発明は、上述した実施の形態の例に限定されるものではなく、本発明の要旨を逸脱しない範囲において、種々の変形、変更が可能である。
【図面の簡単な説明】
【００９８】
【図１】複数時期のレーダ画像の例を示す図である。
【図２】水稲被覆率とレーダ後方散乱係数の関係を示す図である。
【図３】圃場の区画別レーダ後方散乱係数分布図の例を示す図である。
【図４】本発明の一実施の形態に係るシステムの概要を示す図である。
【図５】本発明の一実施の形態に係る植生生育状況解析装置の内部構成を示すブロック図である。
【図６】コーナリフレクタの説明に供する図（１）である。
【図７】コーナリフレクタの説明に供する図（５）である。
【図８】特徴領域のレーダ後方散乱係数を示す図（メディアンフィルタ適用）である。
【図９】特徴領域のレーダ後方散乱係数の経時変化を示す図（メディアンフィルタ適用）である。
【図１０】後方散乱係数の補正の効果（ショッピングセンタ駐車場）を示す図である。
【図１１】後方散乱係数の補正の効果（調整池）を示す図である。
【図１２】本発明の一実施の形態に係る水稲の生育状況解析処理を示すフローチャートである。
【図１３】Ａは補正後のレーダ画像、Ｂは較正線を示す図（８月１１日）である。
【図１４】圃場のレーダ後方散乱係数分布図の例を示す図である。
【図１５】水稲作付圃場の平均後方散乱係数の経時変化を示す図である。
【図１６】パーソナルコンピュータの内部構成の例を示すブロック図である。
【符号の説明】
【００９９】
２０…植生生育状況解析装置、２１…Ｉ／Ｆ部、２２…撮影決定部、２３Ａ…生育カレンダＤＢ、２３Ｂ…撮影計画ＤＢ、２４…レーダ画像保持部、２５…位置合わせ処理部、２６…地図ＤＢ、２７…基準反射板の後方散乱係数処理部、２８…後方散乱係数の判別処理部、２９…地物の後方散乱係数処理部、３０…後方散乱係数保持部、３１…較正用後方散乱係数算出部、３１Ａ…後方散乱係数平滑処理部、３１Ｂ…後方散乱係数局所最小値算出部、３２…較正線算出部、３３…較正線による後方散乱係数較正部、３４…生育値算出部、３５…生育モデルＤＢ、３６…圃場集約算出部、３６Ａ…圃場代表値算出部、３６Ｂ…圃場代表値差分算出部、３６Ｃ…生育判定部、３７…圃場地図ＤＢ、３８…出力部
【技術分野】
【０００１】
本発明は、人工衛星等の飛翔体に搭載され、撮影対象地域の地表面を広範囲にわたり撮影して地表の状況を取得するレーダ装置により得られるレーダ画像データを利用して、植生の生育状況を解析するのに好適な植生生育状況解析方法、プログラムおよび植生生育状況解析装置に係わる。
【背景技術】
【０００２】
従来、農作物の生産高の推定・予測（生育状況解析）は、人工衛星や航空機から取得した光学リモートセンシングデータをもとに、アルゴリズムが構築され、実用化に至ってきた。しかし、日本の場合は、四季や雨期があり年間を通して雨や曇りの日が多いため、天候に左右される光学リモートセンシングでは安定した観測が難しい。
【０００３】
また、日本のように南北に続く複雑な地形で多品種・少量の生産をしている場合は、精密な農業管理が必要である。しかし、これまで圃場単位の生育状況が把握できる高精度なデータとして、光学リモートセンシングデータ以外の入手が困難であるという事情があった。
【０００４】
一方、地表面に照射した電波（マイクロ波パルス）の反射波を測定する能動型センサの例として、合成開口レーダ（SAR：Synthetic Aperture Radar）がある。合成開口レーダは、マイクロ波の特性を利用して、昼夜、天候に関係なく、広域を面的に撮影することができる。また、従来の合成開口レーダは光学リモートセンシングと比較すると解像度（精細度）に問題があったが、最新の衛星搭載型合成開口レーダには、解像度が１ｍ以下のものも存在し、高解像度の画像が得られるようになってきている。
【０００５】
農作物に対する合成開口レーダを用いた既往の研究として、下記のようなものがあった。
・光学リモートセンシングデータと多時期のＳＡＲ画像を併用した水稲作付け面積の推定
・多重波長、多重偏波のＳＡＲ画像からの水稲作付け面積
・多時期のＳＡＲ画像を用いた農作物の種別分類の研究
【０００６】
これらの研究が海外で行われていたり、多数の論文が公表されたりしているが、水稲など農作物の生育予測に直接結びつくものではなかった。
【０００７】
そこで、多時期のＳＡＲ画像を用いた水稲の生育状況について研究が行われている（例えば非特許文献１参照）。
【０００８】
本願の出願人も、多時期のＳＡＲ画像を用いた水稲生育状況に関する研究結果を発表している（非特許文献２を参照。）。非特許文献２に記載された研究により、図１〜図３に示すような解析結果を得た。
【０００９】
図１は、複数時期のレーダ画像例を示すものであり、各観測日ごとの生育と反射特性、水田の様子、水稲被覆状態、レーダ画像を表している。この例では、ＲＡＤＡＲＳＡＴのＣバンドを使用してレーダ画像を撮影している。５月２２日、６月１５日、７月９日、８月２日の４時期のレーダ画像１Ａ，１Ｂ，１Ｃ，１Ｄを作成し、これらのレーダ画像から水稲の生育状況が全体的に変化していることが把握できる。
【００１０】
図２は、水稲被覆率とレーダ後方散乱係数（レーダ装置から送信されたマイクロ波に対する対象物の散乱強度のうち後方への成分）との関係を示したものであり、横軸は水稲被覆率（％）、縦軸はレーダ後方散乱係数（ｄＢ）である。図２に示すように水稲被覆率とレーダ後方散乱係数には強い相関があり、これらの観測点に最小２乗法を適用すると、水稲被覆率とレーダ後方散乱係数との関係を表す回帰直線２が得られる。
【００１１】
図３は、各時期における圃場の区画別レーダ後方散乱係数分布を示したものである。このように、レーダ画像を利用して、時期ごとの圃場の区画別レーダ後方散乱係数分布図３Ａ，３Ｂ，３Ｃ，３Ｄを作成することにより、水稲の生育状況を区画別に把握することができる。
【００１２】
【非特許文献１】Zengyuan Li, Guoqing Sun, Mike Wooding, Yong Pang, Yanfang Dong, Erxue Chen and Bingxiang Tan, “Rice Monitoring Using Envisat Asar Data in China”, Proc. of the 2004 Envisat & ERS Symposium, Salzburg, Austria6-10 September 2004 (ESA SP-572, April 2005)
【非特許文献２】高橋一義、阿部博明、力丸厚、向井幸男、「時系列ＲＡＤＡＲＳＡＴデータを用いた水稲生育分布の把握」、日本写真測量学会 年次学術講演会 東京国際展示場（東京）、２００４年６月１７・１８日
【発明の開示】
【発明が解決しようとする課題】
【００１３】
ところで、非特許文献１には、多時期のＳＡＲ画像を用いた水稲のモニタリングについての研究内容が記載されているが、ＳＡＲ画像の後方散乱係数と植生との関係は明らかになっていない。
【００１４】
一方、非特許文献２に記載された研究により、衛星ＳＡＲデータの後方散乱係数と生育状況との相関が明らかになった。しかしながら、後方散乱係数は撮影時期における大気の状況やシステムの変動などの影響を受けるため、異なる時期に撮影したＳＡＲ画像の後方散乱係数の較正が必要であり、上記研究結果は多時期の水稲生育状況を正確に表してはいなかった。
【００１５】
本発明は、このような状況に鑑みてなされたものであり、人工衛星等の飛翔体から得られるレーダ画像を利用して、多時期の植生の生育状況を正確に解析できるようにするものである。
【課題を解決するための手段】
【００１６】
本発明の一側面の植生生育状況解析方法およびプログラムは、飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得し、取得された複数のレーダ画像を地図データベースに保存する。次に、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行った後、前記レーダ画像の特定領域の後方散乱係数を抽出する。さらに、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する。そして、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する、これらの処理または機能を含むことを特徴とする。
【００１７】
本発明の一側面の植生生育状況解析方法装置は、飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得部と、前記取得部で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持部と、所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理部と、前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理部と、前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正部と、レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正部により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出部とを備えることを特徴とする。
【００１８】
本発明の一側面においては、所定時期のレーダ画像の後方散乱係数を基準として各時期のレーダ画像の後方散乱係数が較正される。そして、較正後の後方散乱係数を基に各時期のレーダ画像が補正されるので、各時期における農作物等の植物の正確な生育状況を示す生育値を算出することができる。
【発明の効果】
【００１９】
以上のように、本発明によれば、飛翔体から得られる多時期のレーダ画像を利用して、多時期の植生の生育状況を正確に解析できる。
【発明を実施するための最良の形態】
【００２０】
以下、本発明を実施するための最良の形態の例について、添付図面を参照しながら説明する。説明は下記項目に従って順に行う。
１．システムの概要
２．植生生育状況解析装置
３．水稲の生育予測処理
４．測定結果等
５．その他
【００２１】
＜１．システムの概要＞
【００２２】
図４は、本発明の植生生育状況解析装置が適用されるシステムの概要を示すものである。
本実施の形態においては、人工衛星にレーダ装置の一例として合成開口レーダ（ＳＡＲ）を搭載し、その合成開口レーダにより（データ解析センタ１２からの指示に基づいて）地表面を撮影する例としてある。以下、人工衛星に搭載した合成開口レーダを「衛星ＳＡＲ」と称す。
【００２３】
人工衛星に搭載された衛星ＳＡＲ１１は、所定の軌道を周回しながら定期的に、もしくはデータ解析センタ１２からの指示に従って随時地表面の圃場１３を撮影し、撮影データ（レーダ画像データ）をデータ解析センタ１２へ送信する。この例では、農作物として水稲が圃場１３に栽培されているとする。
【００２４】
データ解析センタ１２は、生育カレンダに沿って撮影計画を行い、その撮影計画に基づいた撮影指示を含む無線信号を衛星ＳＡＲ１１に送信し、また、衛星ＳＡＲ１１が撮影したレーダ画像データを、アンテナを介して受信する。これに対して後述する植生生育状況解析装置２０で合成開口処理および所定の補正（較正）を行い、補正された多時期の再生画像やレーダ後方散乱係数（単に「後方散乱係数」ともいう。）の分布図等を生成する。
【００２５】
リアルタイムで水稲の生育評価を行うため、本発明では、最初もしくは所定の時期に取得した画像を基準とした、時系列の画像の後方散乱係数値の較正を行う。この後方散乱係数値の較正を行うことにより、撮影時期による変動の影響を考慮した、より実用的な水稲生育予測システムを実現することができる。
【００２６】
時系列の画像の較正方法は、以下の３つの方法が考えられる。
（１） 対象エリアに基準となる反射板（例えばコーナリフレクタ）を設置して、撮影したＳＡＲ画像から反射板の後方散乱係数値を抽出し、これを基準値とする方法
（２） 撮影したＳＡＲ画像の中で、低変動領域を抽出し、その領域の後方散乱係数値を用いて較正するための回帰線（較正線）を求め、これに後方散乱係数分布を当てはめる方法
（３） 上記（１），（２）の方法を併用する方法
【００２７】
＜２．植生生育状況解析装置＞
［植生生育状況解析装置の構成］
図５は、本発明の植生生育状況解析装置の一実施の形態に係る機能ブロック図を示したものである。植生生育状況解析装置は、上記時系列画像の後方散乱係数値の較正を行い、植生の生育状況を予測するためのデータを生成するものである。この植生生育状況解析装置は、データ解析センタ１２に設置されることを想定しているが、衛星ＳＡＲ１１とともに人工衛星に搭載してもよいし、それぞれに分散してもよい。また、植生生育状況解析装置の各機能は、ハードウェアまたはソフトウェアのいずれでも実現することができる。あるいはそれらの組み合わせにより実現してもよい。
【００２８】
本実施の形態の植生生育状況解析装置２０は主に、Ｉ／Ｆ部２１、撮影決定部２２、レーダ画像保持部２４、位置合わせ処理部２５、基準反射板の後方散乱係数処理部２７、後方散乱係数の判別処理部２８、地物の後方散乱係数処理部２９を備える。また、後方散乱係数保持部３０、較正用後方散乱係数算出部３１、較正線算出部３２、較正線による後方散乱係数較正部３３、生育値算出部３４、圃場集約算出部３６、出力部３８を備える。さらに、生育カレンダデータベース（以下、「生育カレンダＤＢ」という。）２３Ａ、撮影計画データベース（以下、「撮影計画ＤＢ」という。）２３Ｂ、地図データベース（以下、「地図ＤＢ」という。）２６、生育モデルデータベース（以下、「生育モデルＤＢ」という。）３５、圃場地図データベース（以下、「圃場地図ＤＢ」という。）３７を備える。
【００２９】
Ｉ／Ｆ部２１は、取得部の一例であり、衛星ＳＡＲ１１で撮影した撮影データを取り込むいわゆるインタフェースとしての機能を備える。Ｉ／Ｆ部２１は、通信部としての機能も備え、撮影決定部２２からの指示を衛星ＳＡＲ１１に送信し、衛星ＳＡＲ１１に撮影を実行させることができる。撮影データは、生育カレンダＤＢ２３Ａに登録された水稲の生育カレンダや、撮影計画ＤＢ２３Ｂに登録された衛星の撮影時期（周期）をもとに、計画的に撮影されるレーダ画像のデータである。
【００３０】
生育カレンダは、例えば水稲の場合であれば、５月の田植え期から始まって分けつ期、最高分けつ期、幼穂形成期、出穂期、登熟期、そして９月の成熟期までの栽培暦である。また、現在打ち上げられている合成開口レーダ搭載の人工衛星には、例えば、観測周期が１１日のものがある。
【００３１】
撮影決定部２２は、これらの生育カレンダと衛星の撮影周期を考慮してレーダ画像をＩ／Ｆ部２１に取り込む時期を決定する。決定した撮影予定日（データ取得日）の一例を以下に示す。
［観測日］ ［起算日］ ［段階］
１）５月 ７日 ０日 田植期
２）５月３１日 ２４日 分けつ期
３）６月２４日 ３８日 最高分けつ期
４）７月１８日 ７２日 幼穂形成期
５）８月１１日 ９６日 出穂期
６）９月 ４日 １２０日 登熟期
【００３２】
なお、撮影決定部２２は、水稲の生育の評価を行う対象面積により、撮影画像の画素サイズを設定するようにしてもよい。例えば、詳細な評価を行う場合は、地上分解能を高く（例えば１ｍ）することが好ましい。
【００３３】
レーダ画像保持部２４は、所定の時期に撮影された時系列のレーダ画像を、地図ＤＢ２６に保存する機能である。地図ＤＢ２６は、例えばハードディスクやフラッシュメモリ等の不揮発性の記憶手段に構築される。
【００３４】
位置合わせ処理部２５は、地図ＤＢ２６に保存されている地図データをもとに、レーダ画像を地図に重ね合わせることができるように、幾何（位置）的な処理を行う機能である。また、位置合わせ後のレーダ画像を、較正線による後方散乱係数較正部３３、基準反射板の後方散乱係数処理部２７および地物の後方散乱係数処理部２９へ供給する。
【００３５】
基準反射板の後方散乱係数処理部２７は、レーダ画像に写っている基準反射板の後方散乱係数を処理（抽出）する機能である。この基準反射板としては、例えば「加美川雄太、大内和夫、“陸域観測技術衛星「だいち」搭載合成開口レーダＰＡＬＳＡＲの校正・検証実験”、平成１８年度卒業論文、高知工科大学工学部」（文献１）に記載されたコーナリフレクタなど、種々の反射板を適用できる。
【００３６】
観測対象エリアに反射板を用いる場合、２時期の反射板の後方散乱係数値をそれぞれ抽出して、基準の画像の後方散乱係数値をもとにした較正を行うための較正線を求め、次の時期や他の時期に撮影した画像を較正する。文献１にも記載されているように、反射板は、地表面に照射するマイクロ波のバンド（波長）や人工衛星の軌道（撮影方向）によって、適宜、適切な面積や形状、個数、方向を調整して設置することが望ましい。
【００３７】
［基準反射板（コーナーリフレクタ）］
図６は、任意の観測対象エリアに、反射板としてコーナリフレクタを設置した場合のレーダ画像と平均後方散乱係数の計測結果を示したものである。図の左側がレーダ画像、右側が後方散乱係数のヒストグラム（頻度分布）を示すグラフである。
【００３８】
図６に示す観測日５月７日のレーダ画像４１に写るコーナリフレクタ４１ａに相当する画素、すなわち後方散乱係数が０ｄＢ以上の画素数は４、平均後方散乱係数が５．１３ｄＢである。同様に、図７に示す観測日８月１１日のレーダ画像４５に写るコーナリフレクタ４５ａは画素数が１４、平均後方散乱係数が５．６１ｄＢである。
【００３９】
このように、時期によってコーナリフレクタの画素数および平均後方散乱係数が異なるので、この情報を後方散乱係数の較正に反映させる。
【００４０】
後方散乱係数の判別処理部２８は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数の利用の可否を判定する機能であり、後方散乱係数が飽和状況の場合、すなわち、レーダ画像の中のコーナリフレクタが設置されている場所以外に、コーナリフレクタにおける後方散乱係数と同じ値の画素が存在する場合は利用しない。後方散乱係数が飽和している場合、正確な後方散乱係数の情報が得られないため精度の良い較正を行うことができないからである。
【００４１】
地物の後方散乱係数処理部２９は、レーダ画像に写っている地物（人工構造物、池、グラウンドなど）の後方散乱係数を抽出する機能である。ここで処理された後方散乱係数をもとにレーダ画像中の低変動領域を決定する。あるいは過去のデータから低変動領域の特徴を示す地物を予め特定しておき、これを利用してもよい。
【００４２】
［複数地点・多時期のレーダ後方散乱係数］
図８、図９は、多地点・多時期で観測したレーダ後方散乱係数及びその経時変化の例を示すものである。なお、図８、図９は、観測結果に平滑処理として後述するメディアンフィルタを適用したものである。
【００４３】
図８では、特徴領域としてショッピングセンタ、仮設住宅、調整池、グラウンド、駐車場（公園）、池（公園）、並びにコーナリフレクタの７地点、また撮影日として５月７日、５月３１日、６月２４日、７月１８日、８月１１日、９月４日において得られたデータを示している。各データは、特徴領域を構成する各画素の平均値をとったものである。
【００４４】
図９は、図８に示したデータを特徴領域のレーダ後方散乱係数の経時変化としてグラフに表したものである。例えば、変化が少なくレーダ後方散乱係数の大きいショッピングセンタや仮設住宅が、レーダ画像の後方散乱係数を較正する際の基準となる低変動領域の候補として挙げられる。また、変化が少なくレーダ後方散乱係数の小さい調整池を低変動領域の候補として挙げることも可能である。これらの低変動領域が写った、最初または所定の時期のレーダ画像が較正時の基準画像として利用される。
【００４５】
図５に示した植生生育状況解析装置の説明に戻る。
後方散乱係数保持部３０は、レーダ画像から抽出された基準反射板および地物の後方散乱係数を保存する機能である。
【００４６】
較正用後方散乱係数算出部３１は、後方散乱係数保持部３０に保存された後方散乱係数に対して、平滑処理および／または局所最小値処理を行い、較正用後方散乱係数を求める機能であり、後方散乱係数平滑処理部３１Ａ、後方散乱係数局所最小値算出部３１Ｂを備える。
【００４７】
後方散乱係数平滑処理部３１Ａは、レーダ画像に対して所定の平滑処理用フィルタ、例えばメディアンフィルタをかけて、スペックルノイズ等のノイズ成分を低減する機能を備える。メディアンフィルタは、ｋ×ｋの局所領域における値を順番に並べ、その真ん中の値を領域中央の画素の出力値とする処理を行うフィルタであり、その処理は式（１）で表すことができる。
y[m，n] = median{x[i,j],(i,j)∈ω}・・・・（１）
ここで、y[m,n]は対象画素の値、x[i,j]は対象画素を中心とする近傍画素の値
【００４８】
後方散乱係数局所最小値算出部３１Ｂは、レーダ画像に対して所定の局所最小値処理用フィルタ（以下、「最小値フィルタ」という。）をかけて、所定領域の最小値を抽出する機能を備える。最小値フィルタは、ｋ×ｋの局所領域における値の最小値を、領域中央の画素の出力値とする処理を行うフィルタであり、その処理は式（２）で表すことができる。
y[m，n] = min{x[i,j],(i,j)∈ω}・・・・（２）
ここで、y[m,n]は対象画素の値、x[i,j]は対象画素を中心とする近傍画素の値
【００４９】
上述のとおり、後方散乱係数平滑処理部３１Ａおよび後方散乱係数局所最小値算出部３１Ｂは、それぞれレーダ画像に発生するノイズ成分の低減および所定領域の最小値を抽出するものである。これらの機能ブロックがない場合、レーダ画像の画質が劣化し、較正の精度が低下することが想定されるが、取得したレーダ画像が要求仕様を満たしてさえいれば必ずしも必要ではなく、また何れか一方のみを用いたり、あるいは処理順序を逆にしたりすることも可能である。
【００５０】
較正線算出部３２は、例えば最初に撮影した画像（第１画像）に対して、順次撮影した画像（第２画像、第３画像・・・）を較正するために、第１画像の較正用後方散乱係数をもとに、第２画像の較正用後方散乱係数との関係式（較正線）を算出する機能である。この較正線算出部３２は、下記較正線による後方散乱係数較正部３３とともに、特許請求の範囲に記載した較正手段を構成する。
【００５１】
較正線による後方散乱係数較正部３３は、後方散乱係数較正手段の一例であり、較正線算出部３２で算出された較正線を用いて、基準画像以外の画像、例えば第２画像、第３画像・・・の全画素の後方散乱係数を処理する機能である。
【００５２】
生育値算出部３４は、生育モデルＤＢ３５に保存された生育モデル（レーダ画像の後方散乱係数と水稲被覆率の相関関係）の情報をもとに、較正線による後方散乱係数較正部３３で求められた後方散乱係数から水稲被覆率を算出する機能である。
【００５３】
圃場集約算出部３６は、圃場地図ＤＢ３７に保存された圃場地図データをもとに、一筆単位の水稲被覆率を算出する機能であり、圃場代表値算出部３６Ａ、圃場代表値差分算出部３６Ｂ、生育判定部３６Ｃを備える。
【００５４】
圃場代表値算出部３６Ａは、１つのレーダ画像を構成する各画素の水稲被覆率から一筆単位の平均値、もしくは、出現頻度が最も高い値（最頻値）を求める機能である。
【００５５】
圃場代表値差分算出部３６Ｂは、基準画像（例えば第１画像）の圃場代表値と参照画像（例えば第２画像）の圃場代表値の差分を計算する機能である。
【００５６】
生育判定部３６Ｃは、圃場代表値および圃場代表値差分の値から生育判定（生育の早い、遅い）を行う機能である。
【００５７】
出力部３８は、水稲被覆率の圃場代表値、圃場代表値差分の値、および、生育判定の結果を、圃場単位の一覧表、もしくは、ランキングした地図（ランキングマップ）などに出力する機能である。
【００５８】
［後方散乱係数の補正の効果］
図１０、図１１は、後方散乱係数の補正の効果を説明するための図である。図１０は、８月１１日におけるあるショッピングセンタ駐車場の場合を示し、図１１は、８月１１日におけるある調整池の場合を示している。
画像５１Ｂは、オリジナルの画像（補正前の画像）５１Ａにメディアンフィルタを実施した後の画像、画像５１Ｃは、メディアンフィルタを実施した後の画像５１Ｂに最小値フィルタを実施した画像である。オリジナルの画像５１Ａはノイズが多いが、メディアンフィルタ処理および最小値フィルタ処理後は、ノイズが除去されるとともに領域中央の画素が周辺の最小値に置き換えられている。この結果、ショッピングセンタ駐車場のように後方散乱係数自体が大きくその変動が少ない領域において、較正用の後方散乱係数を安定的に求めることができる（矢印参照）。
【００５９】
図１１は調整池の場合の例であり、画像５２Ｂは、オリジナルの画像（補正前の画像）５２Ａにメディアンフィルタを実施した後の画像、画像５２Ｃは、メディアンフィルタを実施した後の画像５２Ｂに最小値フィルタを実施した画像である。調整池の場合、図１０に示したショッピングセンタ駐車場の場合と比較して、さらに後方散乱係数自体が小さくその変動が少ない調整池のような領域においても較正用の後方散乱係数を安定的に求めることができる（矢印参照）。
【００６０】
＜３．水稲の生育予測処理＞
［植生生育状況解析装置の動作］
次に、図１２に示すフローチャートを参照して、植生生育状況解析装置２０による水稲の生育状況解析処理を説明する。
図５において、植生生育状況解析装置２０の撮影決定部２２により、適切な撮影時期、場所および地上分解能等を決定し、これらの情報をＩ／Ｆ部２１を介して衛星ＳＡＲ１１へ送信する。Ｉ／Ｆ部２１において、衛星ＳＡＲ１１で撮影されたＳＡＲ画像（レーダ画像）が人工衛星から受信される（ステップＳ１）。受信した各時期のレーダ画像は、レーダ画像保持部２４により地図ＤＢ２６に保存される。
【００６１】
続いて、位置合わせ処理部２５により、地図ＤＢ２６に保存されている最初に撮影したレーダ画像を基準画像として、次または他の時期に撮影したレーダ画像（同一対象エリア）が重なるように位置合わせが行われる（ステップＳ２）。位置合わせ処理では、例えば緯度経度が明確な地上基準点（GCP：Ground Control Point）や道路交差点、建物など画像に表れている点を基準として、地図に重ね合わせることができるように幾何学的な補正が行われる。
【００６２】
位置合わせが終了後、それぞれの場合に応じて、多時期のレーダ画像の後方散乱係数処理および較正が行われる（ステップＳ３）。
【００６３】
第１に、コーナリフレクタ等の基準反射板を利用する場合を説明する。
基準反射板の後方散乱係数処理部２７は、複数時期のレーダ画像中に写るコーナリフレクタの後方散乱係数値をそれぞれ抽出して、基準画像の後方散乱係数値をもとにした較正線を求め、当該較正線に基づいて次の時期や他の時期に撮影したレーダ画像の後方散乱係数値を処理する。このとき、後方散乱係数の判別処理部２８は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数の利用の可否を判定し、後方散乱係数の値が飽和状況の場合は利用しない。一方、後方散乱係数の値が飽和状況でない場合は、基準反射板の後方散乱係数処理部２７で処理した後方散乱係数値を、レーダ画像の撮影時期と対応づけて後方散乱係数保持部３０へ供給する。
【００６４】
そして、後方散乱係数保持部３０に保存された後方散乱係数の情報は、較正用後方散乱係数算出部３１の後方散乱係数平滑処理部３２Ａおよび後方散乱係数局所最小値算出部３２Ｂによってそれぞれフィルタ処理されてノイズ成分が低減され、レーダ画像上の代表的な領域内の後方散乱係数値が求められた後、較正線算出部３２に供給される。較正線算出部３２では、基準画像の低変動領域（内の複数の画素）の後方散乱係数と、その基準画像の低変動領域に対応する他のレーダ画像の特定領域（内の複数の画素）の後方散乱係数から較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６５】
第２に、レーダ画像中の低変動領域の後方散乱係数を用いる場合を説明する。
この場合、較正の基準となる後方散乱係数値は取得したレーダ画像の中から抽出する。地物の後方散乱係数処理部２９は、レーダ画像に写っている地物の後方散乱係数値を計算し、後方散乱係数が大きくその変動が少ない画像領域（例えば、人工構造物）、および、後方散乱係数が低くその変動が少ない画像領域（例えば、池、土のグラウンド）を抽出する。これら後方散乱係数の変動が少ない画像領域、すなわち低変動領域の後方散乱係数値を、レーダ画像の撮影時期と対応付けて後方散乱係数保持部３０へ供給する。
【００６６】
そして、後方散乱係数保持部３０に保存された後方散乱係数の情報は、第１の場合と同様に後方散乱係数平滑処理部３１Ａおよび後方散乱係数局所最小値算出部３１Ｂにてフィルタ処理された後、較正線算出部３２に供給される。較正線算出部３２では、レーダ画像に写る基準反射板の後方散乱係数が求められ、基準画像に写る基準反射板（の複数の画素）の後方散乱係数と、基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板（の複数の画素）の後方散乱係数から較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６７】
第３に、第１項目のコーナリフレクタ等の基準反射板および第２項目のレーダ画像中の低変動領域を併用する場合を説明する。
この場合、較正線算出部３２では、基準画像において後方散乱係数の変動が少ない画像領域の状況に応じて、第２項目におけるレーダ画像上の低変動領域の後方散乱係数のデータと、第１項目におけるコーナリフレクタ等の基準反射板を用いて取得したデータとを選択的に用いて較正線が算出される。すなわち、基準画像から低変動領域を抽出し、これが基準反射板である場合にはその領域の後方散乱係数および基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板の後方散乱係数を用いて、低変動領域が人工構造物等である場合にはその領域の後方散乱係数および基準画像の人工構造物等に対応する他のレーダ画像の特定領域の後方散乱係数を用いて較正線が算出される。この較正線の情報は、較正線による後方散乱係数較正部３３に供給され、この較正線に基づいて、他のレーダ画像の後方散乱係数が較正される。このような較正処理が多時期のレーダ画像について行われる。
【００６８】
そして、生育値算出部３４により生育モデルを基にした水稲の生育解析が行われる（ステップＳ４）。水稲の生育状況の解析は、予め把握している後方散乱係数値と生育情報のモデルを基に、較正処理済みのレーダ画像から計算する。すなわち、生育モデルＤＢ３５に保存された、レーダ画像の後方散乱係数と植生の生育値の相関関係の情報に基づいて、較正されたレーダ画像の後方散乱係数から当該レーダ画像の各画素に写っている植生の生育値が算出される。
【００６９】
水稲の生育解析が終了後、圃場集約算出部３６により解析結果の集計処理が行われる（ステップＳ５）。例えば、圃場代表値算出部３６Ａにより、上記解析結果から、１つのレーダ画像を構成する各画素の水稲被覆率から一筆単位の平均値、もしくは、最頻値が求められた場合、圃場ごとの生育状況が把握できる。
また、圃場代表値差分算出部３６Ｂにより、基準画像（例えば第１画像）の圃場代表値と参照画像（例えば第２画像）の圃場代表値の差分が求められた場合、２つの時期の間で圃場ごとの生育の違いを把握できる。
さらに、生育判定部３６Ｃにより、圃場代表値および圃場代表値差分の値に基づいて生育判定が行われた場合、圃場間での生育の早いか遅いかの違いを把握できる。
【００７０】
そして、出力部３８により水稲生育状況分布の作成が行われる（ステップＳ６）。すなわち、水稲被覆率の圃場代表値、圃場代表値差分の値、および、生育判定の結果を、圃場単位の一覧表、もしくは、ランキングした地図（ランキングマップ）などに出力する処理が行われる。これらの水稲生育状況分布を表した地図は、水稲生育状況分布図と呼ばれる。
【００７１】
水稲生育状況分布図は、ピクセル（画素）単位、圃場単位、圃場の中を細分化した単位など、営農指導を行う単位ごとに情報の集計処理を行うことにより、個々のユーザの実情に合わせた営農指導を行える。
【００７２】
ピクセル（画素）単位から圃場単位に集計するためには、圃場の形状のデータ、例えば、地理情報データを利用する。対象エリアの地理情報データは、地図ＤＢ２６等に保存しておく。
【００７３】
＜４．測定結果等＞
【００７４】
図１３は、所定の時期（この例では５月７日）のレーダ画像を基に、８月１１日に取得したデータを較正（図１２のステップＳ３）した結果を示すものである。図１３Ａは対象時期の較正後のレーダ画像、図１３Ｂは図１３Ａに示したレーダ画像の較正に用いた較正線を示す図である。なお、本例では後方散乱係数が最小となる撮影時期の最初（５月７日）を基準にグラフを作成したが、後述するように後方散乱係数が最大となる時期である８月１１日を基準として標準化することも可能である。
【００７５】
図１３Ｂは、８月１１日のデータを較正するために用いる、５月７日のレーダ画像を基にした較正線である。縦軸が５月７日の後方散乱係数であり、横軸が８月１１日における後方散乱係数である。この例では、較正線（回帰直線）は下記式で表される。ただし、ｙは基準画像（５月７日）の後方散乱係数、ｘは較正する画像（８月１１日）の後方散乱係数、Ｒ２は決定係数である。なお、図１３Ｂに示したデータは、最小値フィルタを適用後のデータである。
ｙ＝０．９８４ｘ−０．１３７
Ｒ２＝０．９７２
【００７６】
同様にして、所定時期（この例では５月７日）のレーダ画像の後方散乱係数を基準として作成した各時期の較正線により、それぞれの時期に取得したレーダ画像の後方散乱係数が較正される。この結果、レーダ画像の後方散乱係数が較正されるので、各時期における水稲（農作物等の植物）の正確な生育状況を示す生育値を算出することができる。
【００７７】
次に、圃場のレーダ後方散乱係数分布図について説明する。
図１４は、ある圃場（水稲作付圃場）についてのレーダ後方散乱係数の分布を示したものであり、複数の圃場を含む対象エリア全体のレーダ後方散乱係数分布図の要部の拡大図である。図１４において、拡大されたレーダ後方散乱係数分布図には圃場７１〜７３が含まれている。このように、レーダ画像の任意の領域に対応する農作物の生育状況について、画素単位で精細に確認することができる。
【００７８】
図１５は、図１４に示した水稲作付圃場のレーダ後方散乱係数分布図に対応する平均後方散乱係数の経時変化を示すグラフである。平均後方散乱係数とは、レーダ画像内のある領域を構成する複数の画素から得られた後方散乱係数の平均値を算出したものである。圃場７１〜７３の水稲は、田植え後から徐々に成長した結果として、後方散乱係数が大きくなり、８月１１日を極大点として後方散乱係数が小さくなっている。これは、出穂および稲の垂下によって後方散乱が減少した結果と推測される。この図１５のグラフからは、例えば圃場７３の水稲が早熟で、圃場７１の水稲は晩熟であることが把握できる。
【００７９】
本発明の観測の対象は、水稲だけでなくトウモロコシや大麦等の生育段階に伴う形態変化に応じて後方散乱係数が変化する農作物、植物が含まれ、地表面の植生について基準画像を用いて補正した複数の時期のレーダ画像や、補正後のレーダ画像から作成したレーダ後方散乱係数分布図（植生の生育状況分布図）を観察することにより、植生の生育状況をリアルタイムかつ正確に把握できる。
【００８０】
上述した実施の形態によれば、所定時期のレーダ画像の後方散乱係数を基準として各時期のレーダ画像の較正線を算出し、算出した各時期の較正線により、それぞれの時期に取得したレーダ画像の後方散乱係数が較正される。そして、較正後の後方散乱係数を基に各時期のレーダ画像が補正されるので、各時期における農作物等の植物の正確な生育状況を知ることができる。
【００８１】
また、人工衛星等の飛翔体からマイクロ波を用いて撮像したレーダ画像は、全天候型であるために、５月（田植え）から８月（出穂）の間の雨期でも農作物等の植物の生育状況を把握することができる。
【００８２】
また、多時期にわたる正確でリアルタイムな生育状況の把握により、農作物等の植物の適時かつ適切な栽培管理、追肥の投入時期の決定等が可能となる。
【００８３】
また、多時期にわたる正確でリアルタイムな生育状況の把握により、病害虫などの発生の予防、生育障害、災害への対処が可能となる。
【００８４】
さらに、官公庁等で実施している現地調査による水稲作況調査が観測手段（リモートセンシング）で可能となる。
【００８５】
＜５．その他＞
【００８６】
なお、上述した植生生育状況解析装置で行われる一連の処理は、ハードウェアにより実行させることもできるし、ソフトウェアにより実行させることもできる。また、これらの処理を実行する機能はハードウェアとソフトウェアの組み合わせによっても実現できることは言うまでもない。一連の処理をソフトウェアにより実行させる場合には、そのソフトウェアを構成するプログラムが専用のハードウェアに組み込まれているコンピュータ、または、各種のプログラムをインストールすることで各種の機能を実行することが可能な、例えば汎用のコンピュータなどに、プログラム記録媒体からインストールされる。
【００８７】
図１６は、上述した一連の処理をプログラムにより実行するコンピュータの構成例を示すブロック図である。このコンピュータ１０１は、例えば一連の処理を実行するために高性能化した専用コンピュータの他、一定の性能を備えるパーソナルコンピュータなどでもよい。
【００８８】
コンピュータ１０１のＣＰＵ（Central Processing Unit）１１１は、ＲＯＭ（Read Only Memory）１１２、または記録部１１８に記録されているプログラムに従って、上記一連の処理の他、各種の処理を実行する。ＲＡＭ（Random Access Memory）１１３には、ＣＰＵ１１１が実行するプログラムやデータなどが適宜記憶される。これらのＣＰＵ１１１、ＲＯＭ１１２、およびＲＡＭ１１３は、バス１１４により相互に接続されている。
【００８９】
ＣＰＵ１１１にはまた、バス１１４を介して入出力インタフェース１１５が接続されている。入出力インタフェース１１５には、キーボード、マウス、マイクロホンなどよりなる入力部１１６、ディスプレイ、スピーカなどよりなる出力部１１７が接続されている。ＣＰＵ１１１は、入力部１１６から入力される指令に対応して各種の処理を実行する。そして、ＣＰＵ１１１は、処理の結果を出力部１１７に出力する。
【００９０】
入出力インタフェース１１５に接続されている記録部１１８は、例えばハードディスクからなり、ＣＰＵ１１１が実行するプログラムや各種のデータを記録する。
【００９１】
通信部１１９は、インターネットやローカルエリアネットワークなどのネットワークを介して外部の装置と通信する。また通信部１１９を介してプログラムを取得し、記録部１１８に記録してもよい。
【００９２】
入出力インタフェース１１５に接続されているドライブ１２０は、磁気ディスク、光ディスク、光磁気ディスク、あるいは半導体メモリなどのリムーバブルメディア１３１が装着されたとき、それらを駆動し、そこに記録されているプログラムやデータなどを取得する。取得されたプログラムやデータは、必要に応じて記録部１１８に転送され、記録される。
【００９３】
コンピュータにインストールされ、コンピュータによって実行可能な状態とされるプログラムを格納するプログラム記録媒体は、図３０に示すように、リムーバブルメディア１３１によりパッケージメディアとして提供される。リムーバブルメディア１３１としては、磁気ディスク（フレキシブルディスクを含む）、光ディスク（ＣＤ−ＲＯＭ（Compact Disc - Read Only Memory），ＤＶＤ（Digital Versatile Disc），光磁気ディスクを含む）、もしくは半導体メモリなどを適用することができる。あるいは、プログラム記録媒体は、プログラムが一時的もしくは永続的に格納（記録）されるＲＯＭ１１２や、記録部１１８を構成するハードディスクなどにより構成される。
【００９４】
このプログラム記録媒体へのプログラムの格納は、必要に応じてルータ、モデムなどのインタフェースである通信部１１９を介して、ローカルエリアネットワーク(LAN：Local Area Network)、インターネット、デジタル衛星放送といった、有線または無線の通信媒体を利用して行われる。
【００９５】
なお、本明細書において、プログラム記録媒体に格納されるプログラムを記述する処理ステップは、記載された順序に沿って時系列的に行われる処理はもちろん、必ずしも時系列的に処理されなくとも、並列的あるいは個別に実行される処理（例えば、並列処理あるいはオブジェクトによる処理）をも含むものである。
【００９６】
また、プログラムは、一つのコンピュータにより処理されるものであってもよいし、複数のコンピュータによって分散処理されるものであってもよい。さらに、プログラムは、遠方のコンピュータに転送されて実行されるものであってもよい。
【００９７】
以上に述べた実施の形態は、本発明を実施するための好適な形態の具体例であるから、技術的に好ましい種々の限定が付されている。ただし、本発明は、以上の実施の形態の説明において特に本発明を限定する旨の記載がない限り、これらの実施の形態に限られるものではない。したがって、例えば、以上の説明で挙げた処理時間、処理順序および各パラメータの数値的条件等は好適例に過ぎず、また、説明に用いた各図における寸法、形状および配置関係等も実施の形態の一例を示す概略的なものである。したがって、本発明は、上述した実施の形態の例に限定されるものではなく、本発明の要旨を逸脱しない範囲において、種々の変形、変更が可能である。
【図面の簡単な説明】
【００９８】
【図１】複数時期のレーダ画像の例を示す図である。
【図２】水稲被覆率とレーダ後方散乱係数の関係を示す図である。
【図３】圃場の区画別レーダ後方散乱係数分布図の例を示す図である。
【図４】本発明の一実施の形態に係るシステムの概要を示す図である。
【図５】本発明の一実施の形態に係る植生生育状況解析装置の内部構成を示すブロック図である。
【図６】コーナリフレクタの説明に供する図（１）である。
【図７】コーナリフレクタの説明に供する図（５）である。
【図８】特徴領域のレーダ後方散乱係数を示す図（メディアンフィルタ適用）である。
【図９】特徴領域のレーダ後方散乱係数の経時変化を示す図（メディアンフィルタ適用）である。
【図１０】後方散乱係数の補正の効果（ショッピングセンタ駐車場）を示す図である。
【図１１】後方散乱係数の補正の効果（調整池）を示す図である。
【図１２】本発明の一実施の形態に係る水稲の生育状況解析処理を示すフローチャートである。
【図１３】Ａは補正後のレーダ画像、Ｂは較正線を示す図（８月１１日）である。
【図１４】圃場のレーダ後方散乱係数分布図の例を示す図である。
【図１５】水稲作付圃場の平均後方散乱係数の経時変化を示す図である。
【図１６】パーソナルコンピュータの内部構成の例を示すブロック図である。
【符号の説明】
【００９９】
２０…植生生育状況解析装置、２１…Ｉ／Ｆ部、２２…撮影決定部、２３Ａ…生育カレンダＤＢ、２３Ｂ…撮影計画ＤＢ、２４…レーダ画像保持部、２５…位置合わせ処理部、２６…地図ＤＢ、２７…基準反射板の後方散乱係数処理部、２８…後方散乱係数の判別処理部、２９…地物の後方散乱係数処理部、３０…後方散乱係数保持部、３１…較正用後方散乱係数算出部、３１Ａ…後方散乱係数平滑処理部、３１Ｂ…後方散乱係数局所最小値算出部、３２…較正線算出部、３３…較正線による後方散乱係数較正部、３４…生育値算出部、３５…生育モデルＤＢ、３６…圃場集約算出部、３６Ａ…圃場代表値算出部、３６Ｂ…圃場代表値差分算出部、３６Ｃ…生育判定部、３７…圃場地図ＤＢ、３８…出力部
【特許請求の範囲】
【請求項１】
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得ステップと、
前記取得ステップで取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持ステップと、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理ステップと、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理ステップと、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正ステップと、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正ステップにより較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する生育値算出ステップと、
を含む植生生育状況解析方法。
【請求項２】
前記較正ステップは、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数と前記基準画像の特定領域に対応する他のレーダ画像の特定領域の後方散乱係数から、他のレーダ画像の後方散乱係数を較正するための較正線を算出する較正線算出ステップと、
前記較正線算出ステップにより算出された較正線に基づいて、他のレーダ画像の後方散乱係数を較正する後方散乱係数較正ステップと、をさらに含む
請求項１に記載の植生生育状況解析方法。
【請求項３】
前記後方散乱係数処理ステップでは、時期による後方散乱係数の変動の少ない低変動領域をレーダ画像から抽出し、
前記較正線算出ステップでは、前記基準画像の低変動領域の後方散乱係数と、前記基準画像の低変動領域に対応する他のレーダ画像の特定領域の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項４】
前記対象エリアの地表面の所定位置に基準反射板が設けられ、
前記後方散乱係数処理ステップでは、前記レーダ画像に写る基準反射板の後方散乱係数を求め、
前記較正線算出ステップでは、前記基準画像に写る基準反射板の後方散乱係数と、前記基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項５】
前記較正線算出ステップでは、前記基準画像の低変動領域の後方散乱係数および前記基準画像の低変動領域に対応する他のレーダ画像の特定領域の後方散乱係数、並びに、前記基準画像に写る基準反射板の後方散乱係数および前記基準画像に写る基準反射板に対応する前記他のレーダ画像に写る基準反射板の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項６】
前記所定時期は前記複数の時期のうちの最初である
請求項１〜５のいずれかに記載の植生生育状況解析方法。
【請求項７】
前記後方散乱係数処理ステップにより求められた前記レーダ画像に写る基準反射板の後方散乱係数が飽和状況にあるか否かを判定し、前記後方散乱係数が飽和状況にある場合には前記後方散乱係数を前記較正線算出ステップによる較正線の算出に利用しないことを決定する判別処理ステップを、さらに含む
請求項４または５に記載の植生生育状況解析方法。
【請求項８】
圃場地図を保存する圃場地図データベースと、
前記圃場地図データベースに保存された圃場地図を基に、圃場単位の生育値を算出する圃場集約算出ステップと、をさらに含む
請求項１〜７のいずれかに記載の植生生育状況解析方法。
【請求項９】
前記植生の生育値は、水稲の被覆率である
請求項１〜８のいずれかに記載の植生生育状況解析方法。
【請求項１０】
コンピュータを、
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得手段と、
前記取得手段で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持手段と、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理手段と、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理手段と、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正手段と、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正手段により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出手段
として機能させるためのプログラム。
【請求項１１】
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得部と、
前記取得部で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持部と、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理部と、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理部と、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正部と、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正部により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出部と、
を備える植生生育状況解析装置。
【請求項１】
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得ステップと、
前記取得ステップで取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持ステップと、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理ステップと、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理ステップと、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正ステップと、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正ステップにより較正されたレーダ画像の後方散乱係数から当該レーダ画像に写っている植生の生育値を算出する生育値算出ステップと、
を含む植生生育状況解析方法。
【請求項２】
前記較正ステップは、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数と前記基準画像の特定領域に対応する他のレーダ画像の特定領域の後方散乱係数から、他のレーダ画像の後方散乱係数を較正するための較正線を算出する較正線算出ステップと、
前記較正線算出ステップにより算出された較正線に基づいて、他のレーダ画像の後方散乱係数を較正する後方散乱係数較正ステップと、をさらに含む
請求項１に記載の植生生育状況解析方法。
【請求項３】
前記後方散乱係数処理ステップでは、時期による後方散乱係数の変動の少ない低変動領域をレーダ画像から抽出し、
前記較正線算出ステップでは、前記基準画像の低変動領域の後方散乱係数と、前記基準画像の低変動領域に対応する他のレーダ画像の特定領域の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項４】
前記対象エリアの地表面の所定位置に基準反射板が設けられ、
前記後方散乱係数処理ステップでは、前記レーダ画像に写る基準反射板の後方散乱係数を求め、
前記較正線算出ステップでは、前記基準画像に写る基準反射板の後方散乱係数と、前記基準画像に写る基準反射板に対応する他のレーダ画像に写る基準反射板の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項５】
前記較正線算出ステップでは、前記基準画像の低変動領域の後方散乱係数および前記基準画像の低変動領域に対応する他のレーダ画像の特定領域の後方散乱係数、並びに、前記基準画像に写る基準反射板の後方散乱係数および前記基準画像に写る基準反射板に対応する前記他のレーダ画像に写る基準反射板の後方散乱係数から較正線を算出する
請求項２に記載の植生生育状況解析方法。
【請求項６】
前記所定時期は前記複数の時期のうちの最初である
請求項１〜５のいずれかに記載の植生生育状況解析方法。
【請求項７】
前記後方散乱係数処理ステップにより求められた前記レーダ画像に写る基準反射板の後方散乱係数が飽和状況にあるか否かを判定し、前記後方散乱係数が飽和状況にある場合には前記後方散乱係数を前記較正線算出ステップによる較正線の算出に利用しないことを決定する判別処理ステップを、さらに含む
請求項４または５に記載の植生生育状況解析方法。
【請求項８】
圃場地図を保存する圃場地図データベースと、
前記圃場地図データベースに保存された圃場地図を基に、圃場単位の生育値を算出する圃場集約算出ステップと、をさらに含む
請求項１〜７のいずれかに記載の植生生育状況解析方法。
【請求項９】
前記植生の生育値は、水稲の被覆率である
請求項１〜８のいずれかに記載の植生生育状況解析方法。
【請求項１０】
コンピュータを、
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得手段と、
前記取得手段で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持手段と、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理手段と、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理手段と、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正手段と、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正手段により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出手段
として機能させるためのプログラム。
【請求項１１】
飛翔体に搭載されたレーダ装置で複数の時期に撮影された、同一対象エリアの地表面のレーダ画像を取得する取得部と、
前記取得部で取得された複数のレーダ画像を地図データベースに保存するレーダ画像保持部と、
所定時期に撮影されたレーダ画像を基準画像として、他の時期に撮影されたレーダ画像の位置合わせを行う位置合わせ処理部と、
前記レーダ画像の特定領域の後方散乱係数を抽出する後方散乱係数処理部と、
前記地図データベースに保存された複数のレーダ画像のうち、前記基準画像の特定領域の後方散乱係数を基に、他のレーダ画像の後方散乱係数を較正する較正部と、
レーダ画像の後方散乱係数と植生の生育値の相関関係に基づいて、前記較正部により較正されたレーダ画像の後方散乱係数から当該レーダ画像に写る植生の生育値を算出する生育値算出部と、
を備える植生生育状況解析装置。
【図２】
【図４】
【図５】
【図８】
【図９】
【図１２】
【図１５】
【図１６】
【図１】
【図３】
【図６】
【図７】
【図１０】
【図１１】
【図１３】
【図１４】
【図４】
【図５】
【図８】
【図９】
【図１２】
【図１５】
【図１６】
【図１】
【図３】
【図６】
【図７】
【図１０】
【図１１】
【図１３】
【図１４】
【公開番号】特開２０１０−１１７３２７（Ｐ２０１０−１１７３２７Ａ）
【公開日】平成２２年５月２７日（２０１０．５．２７）
【国際特許分類】
【出願番号】特願２００８−２９２７０５（Ｐ２００８−２９２７０５）
【出願日】平成２０年１１月１４日（２００８．１１．１４）
【出願人】（０００１３５７７１）株式会社パスコ (102)
【出願人】（３０４０２１２８８）国立大学法人長岡技術科学大学 (458)
【Ｆターム（参考）】
【公開日】平成２２年５月２７日（２０１０．５．２７）
【国際特許分類】
【出願日】平成２０年１１月１４日（２００８．１１．１４）
【出願人】（０００１３５７７１）株式会社パスコ (102)
【出願人】（３０４０２１２８８）国立大学法人長岡技術科学大学 (458)
【Ｆターム（参考）】
[ Back to top ]