JP5795932B2

JP5795932B2 - 土地状態推定方法

Info

Publication number: JP5795932B2
Application number: JP2011222455A
Authority: JP
Inventors: 洋輝大橋; 俊宏鯨井
Original assignee: Hitachi Ltd
Current assignee: Hitachi Ltd
Priority date: 2011-10-07
Filing date: 2011-10-07
Publication date: 2015-10-14
Anticipated expiration: 2031-10-07
Also published as: JP2013084072A

Description

本発明は、衛星画像を用いて土地の状態を推定する土地状態推定方法に関する。

自然災害や事故などのなんらかのイベントが生じたときに、その土地の状態を広域的に把握したいという需要がある。

これに対し、従来から、対象地の現地調査を行いその結果を集約するような方法がとられている。ところが、これには対象地域の一点一点を実際に訪れて調査をしなければならないために、時間とコストがかかってしまうという問題がある。また、現地に近づくことが危険で、現地調査を行うこと自体が困難であることもある。

これらの問題に対応するために、衛星画像を用いて、その地域の状態を広域的に把握しようという試みがなされている。

例えば、特許文献１には、人工衛星に搭載した合成開口レーダによって災害が起こる前の平常時に、撮影対象の地表面を撮影しレーダ画像データを取得すること、また、災害発生後には回帰軌道日数よりも短い日数内に撮影対象の地表を撮影し、定常時のレーダ画像データと比較することで早急な被災状況の把握に努めること、さらに、災害発生後の復旧・復興時にも定期的に撮影対象の撮影を行うことにより、復旧・復興計画や２次災害の防止に役立てることが記載されている。また、特許文献２には、地表面の変化を正確に抽出することができるレーダ画像処理装置及びレーダ画像処理方法を提供することが記載されている。

国際公開ＷＯ２００８／０１６１５３号公報特開２００８−４６１０７号公報

McLeod, M. and Slavich, P. and Iskandar, T. and Rachman, A. and Hunt, C. and Irhas, C. and Ali, N. and Moore, N.:Soil Salinity Assessment in Tsunami-Affected Areas of Aceh Using the Electromagnetic Induction Method, International Workshop on Post Tsunami Soil Management, 1-2 July 2008 in Bogor, Indonesia

特許文献２では、異なる時刻に取得された画像から、例えば、合成開口レーダの散乱強度の差分、散乱強度の適当な周辺画素領域内の標準偏差の差分、相関指標及び高さ情報の差分などの、複数の指標値を統合的に用いることで、正確に変化の起きた領域を抽出する技術が説明されている。ところが、２時期の画像から変化域を抽出することに焦点が置かれており、３時期以上の画像を利用して２時期の画像からは分からなかった情報を抽出する手法については明記されていない。

特許文献１では、災害発生前の衛星画像と災害発生直後の衛星画像から変化域を抽出し、応急対応時の変化抽出図を作成する方法と、災害発生直後の衛星画像と災害発生後に定期的に撮影されたレーダ画像データとの差分から変化域を抽出し、復旧・復興時の変化抽出図を作成する方法が記載されている。これらの方法では、３時期の衛星画像が用いられてはいるものの、これは応急時と復旧・復興時を分けて考え、それぞれ変化抽出を行ったものであり、実質的に変化抽出を行う際には３時期の衛星画像を統合的に用いているのではなく、２時期の衛星画像から変化抽出を行っている。

これらの手法では、イベント発生前からイベント発生直後、そしてさらにその後の時期を通してある同一の指標がどのように変化していったかを考慮に入れていないため、これらの時期を通して継続的に同一の指標値の変化の様子を見て初めて分かるような土地の状態については推定ができないという問題点がある。

本発明は、これらの問題に対応するために、３時期以上の衛星画像において、各時期の画像の各地点に対応するピクセルごとに計算できる任意の同一指標の変化を用いて、対応する土地の状態を推定できる土地状態推定方法を提供することを目的とする。

本発明に係る土地状態推定方法は、衛星画像データを格納する格納部を有する計算機における土地状態推定方法である。

より具体的には、本発明の土地状態推定方法は、格納された複数の衛星画像データから、所定の土地領域の津波の発性以前の第１時期の第１の衛星画像データ、津波の発生直後の第２時期の第２の衛星画像データ、及び津波の発生から時間が経過した第３時期の第３の衛星画像データを抽出し、前記第１の衛星画像データ、第２の衛星画像データ、及び第３の衛星画像データから、それぞれ画像のピクセル毎に土壌水分量指標値を算出し、第１の衛星画像データから求めた土壌水分量指標値と第２の衛星画像データから求めた土壌水分量指標値との差分もしくは比率を示す第１の変化指標と、第２の衛星画像データから求めた土壌水分量指標値と第３の衛星画像データから求めた土壌水分量指標値との差分もしくは比率を示す第２の変化指標とを算出し、画像のピクセル毎に、第１の変化指標が正の第１の閾値より大であり、かつ第２の変化指標が正の第２の閾値より大である第１のパターン、第１の変化指標が前記正の第１の閾値より大であり、かつ第２の変化指標が負の第３の閾値と前記正の第２の閾値の間である第２のパターン、及び第１の変化指標が前記正の第１の閾値より大であり、かつ第２の変化指標が前記負の第３の閾値より小である第３のパターンのいずれに当るか判別し、判別されたそれぞれのパターンごとに画像を色分けして表示することを特徴とする。

上記第１から第３のパターンは、海水を被った程度に加えて、土地の排水性の程度を、ひいては土地に浸透した塩分がその後に流出した程度を反映する。したがって、それぞれのパターンごとに画像を色分けして表示することで塩害の大小の分布を把握できる画像を提供することができる。

本発明によって、３時期以上の衛星画像を統合的に用いることで初めて分かるような土地の状態を広域的に推定することができるようになり、これにより対策計画の立案支援などを行うことが可能になる。

処理を行う装置の構成図の例である。機能構成図の例である。３時期の衛星画像におけるNDWI値から定性的な推定モデルを用いて塩分濃度推定を行う方法を説明するフローチャートの例である。定性的に塩分濃度推定を行う際に参照する表の例である。３時期の衛星画像におけるNDWI値の３つの変化パターンのみに着目して塩分濃度推定を行い、その結果を地図上に可視化する方法を説明するフローチャートの例である。特に着目すべきNDWI値の変化パターンの例である。定性的に植生被害度推定を行う際に参照する表の例である。特に着目すべきNDVI値の変化パターンの例である。

以下、実施例を、図面を用いて説明する。

本実施例では、津波発生後の沿岸地域における土壌の塩分濃度を推定する塩分濃度推定システムについて説明する。

塩分濃度推定システムにおける注目指標値としては、土壌水分量と相関のある指標、例えばNDWI(Normalized Difference Water Index)を用いる。この他に、例えば合成開口レーダから得られる後方散乱強度などを用いてもよい。NDWIは土壌に含まれる水分量を表す指標であり、具体的には可視域の光(630nm前後)の輝度VISと、中間赤外域の光(1620nm前後)の輝度MIRを用いて、NDWI=(VIS-MIR)/(VIS+MIR)と計算される値である。

津波発生前後や数ヵ月後のNDVI値から推定される土壌水分量の変化から、土壌の排水性が推定できる。土壌の排水性は津波被害を受けた後、どの程度の期間、水がはけずに浸水していたかと関係がある。そして、津波によって海水を被った土地が浸水していた期間の長さと、その地域が被る塩害の程度の間には、相関関係があることが知られている。非特許文献１では2004年12月26日のスマトラ島沖地震によって発生したインド洋大津波によって被害を受けた地域において、長期間浸水していた場所ほど塩分濃度の指標となる土壌EC(Electrical Conductivity)の値が高くなる傾向があることが報告されている。これは、津波の被害を受けても短期間で水がはけるような土地では、土壌の排水性が良いため、降雨などとともに塩分が流れるのに対し、長期間浸水するような土地では、塩分が流れず結果的に塩分濃度が高くなっているためと考えられる。以上のような理由からNDVI値から推定される土壌水分量の変化から、土壌の排水性が推定でき、ひいては塩害の程度が推定できる。

以下に具体的に衛星画像におけるNDWI値から塩分濃度を推定する方法の一例を説明する。まず３時期の衛星画像から定性的な推定モデルを用いて塩分濃度推定を行う方法について述べ（以下これをケースＡと記載する）、次に３時期以上の衛星画像を用いて定量的な推定モデルを用いて塩分濃度推定を行う方法について述べる（以下これをケースＢと記載する）。

ケースＡでは、津波発生前の衛星画像、津波発生直後の衛星画像、津波発生後数カ月の衛星画像１枚の合計３枚の衛星画像を用いる。ケースＡの説明においては、今後これらをそれぞれ画像１・画像２・画像３と呼ぶこととする。

画像１・画像２・画像３それぞれを複数の画像から選択できる場合は、なるべく、指標値に影響を与えるようなノイズが少なく、かつ選択された画像群における前提条件が同一になるようにする必要がある。例えば、塩分濃度推定システムでは指標値として土壌水分量指標値を用いるが、直前に雨が降っている場合は土壌水分量指標値へ一定の影響を及ぼすことが考えられるため、そのような時点で撮影された画像を避けて選択することが好適である。同様に、天候条件や撮影の時間帯による差異が結果に悪影響を与える可能性があるため、なるべく前提条件が同一になるような画像群を選択する、あるいは、例えば反射率のような異なる条件下でも影響を受けない指標に換算してから処理を行うことが好適である。本実施例においては、水田に水がはられている場合に、その部分の土壌水分量指標値が高くなるといった類のノイズも考え得るが、その場合は対象地域の圃場情報を利用して、水田の場所をフィルタリングし、その部分だけ別途処理を行うなどの対応をするとよい。例えば、水田に該当する部分のピクセルに関しては、後に説明する塩分濃度の判定基準となる土壌水分量指標値の閾値を、他の土地とは別の値に設定するなどとする。

図１は土地状態推定システムを実現するための装置の例である。一連の処理の流れとしては、例えば、補助記憶装置105に保存されている衛星画像を主記憶装置104に読み出し、中央演算処理装置103で土地状態推定のための処理を行い、その処理結果を出力制御装置106を通してディスプレイなどの出力装置120へ表示する。

図２は土地状態推定システムを実現するための処理部の機能構成図の例である。まず衛星画像読込部201で衛星画像ストレージ210から処理対象の衛星画像を読み込む。衛星画像ストレージの記憶媒体としては、ハードディスクやDVDなどが例としてあげられるが、これらに限るものではない。次に、指標値計算部202で、各衛星画像から指標値を計算する。また、衛星画像データは、衛星画像が撮影された時刻、又は、衛星画像が土地状態推定システムに取り込まれた時刻を示す情報と対応付けて衛星画像ストレージ２１０に格納される。また、衛星画像読込部201により読み込まれる衛星画像は、ユーザから入力装置１１０を介して指定されてもよいし、ユーザが指定する時刻情報に基づいて指定されてもよい。本実施例の場合は、具体的にはNDWI値を計算する。続いて、指標値→土地状態特徴量変換部203で、指標値変換用参照物220を参照しながら、指標値計算部202で計算された指標値を、土地状態を表す特徴量に変換する。指標値変換用参照物220の例としては、後に説明する図４や図７に記載の変換表や、式１に記載の式とパラメータの組み合わせなどがあげられるが、これらに限るものではなく、指標値から土地状態を表す特徴量に変換をするための参照物ならばその形は問わない。本実施例では、具体的にはNDWI値から塩分濃度へ変換するための参照表や式などを用いる。最後に、結果表示部204で、推定された土地状態特徴量を表示する。この際の表示方法に制限はなく、例えば、表などの形を用いて土地状態特徴量を表示してもよいし、土地状態特徴量を量子化したものを画像の輝度値に割り当てて表示してもよい。

図３はケースＡの処理のフローチャートの例である。注目している衛星画像中にN個のピクセルが含まれているとすると、まずそれらに0からN-1までの番号を一意に振る。ここで、番号づけの順序に制約はなく任意に定めてよい。S101では注目対象のピクセル番号nを0に初期化する。次にS102で画像i(i=1,2,3)のNDWI値NDWIinをそれぞれ計算する。次にS103において、例えば、図４に記載の表のような、３時期のNDWI値を塩分濃度に変換するための表を参照してピクセルnに対応する土地の塩分濃度を決定する。あるいは、必ずしも表を用いなくとも、例えば図４に記載の表を参照して行うのと同等の変換を行う関数を用いて変換するなど、３時期の衛星画像それぞれから計算されるピクセルごとのNDWI値を塩分濃度へ変換する方法であれば、そのデータ構造や表現手段などは問わない。このようにして、n番目のピクセルの塩分濃度が求められたら、次にS104で注目ピクセルが最後の番号であるN-1まで到達しているかを判定し、まだ最後まで達していなければS105でnをn+1に更新し、S102へ戻る。最後のピクセルに達していればS106で処理を終了する。

ここで、図４に示した表の妥当性について述べる。S103で利用したNDWI値を塩分濃度へ変換するための表や関係式などは、３時期のNDWI値と対応する土地の塩分濃度の関係を妥当にあらわしたものでなければならない。例えば図４では、津波発生前はNDWI値が低かったものの津波発生直後にNDWI値が大きくなり、数カ月後もなお高いNDWI値を持つような土地の塩分濃度を大としている。これは、このような土地は、津波発生前後でNDWI値が上昇していることから、確かに津波被害を受け、なおかつ、数ヵ月後もまだNDWI値が高いことから、排水性が低く、塩分が蓄積されやすいと考えられるためである。また、津波発生前はNDWI値が低かったものの津波発生直後にNDWI値が大きくなり、数カ月後には再びNDWI値が下がっているような土地の塩分濃度を中としている。これは、このような土地は、津波発生前後でNDWI値が上昇していることから、確かに津波被害を受けているためある程度塩分濃度は高くなるが、数ヵ月後にはNDWI値が下がっていることから、土地の排水性は高く、塩分が蓄積されにくいと考えられるからである。その他、津波発生直後でもNDWI値が低かったような土地は、そもそも津波被害を受けていないと考えられるので塩分濃度を小としている。また、津波発生直後のNDWI値は高いが津波発生前のNDWI値もまた高いような土地は、もともと水分を多く含んでいる土地であり、津波の被害を受けていないか、受けていたとしても、他の乾いた土地に比べ相対的に津波によってもたらされる海水のために塩分濃度が高まる危険性は低いため、塩分濃度は小としている。

次に、図４の表中に記載の各時期のNDWI値の高低の判断の仕方について、例を示して説明する。閾値Tiを決める際には、例えば、津波発生前の数時期の衛星画像から通常時のNDWI値を計算しておき、その平均から１標準偏差高い値をTiに設定するなどの方法が考えられる。その他に、その地域のNDWI値の平均値を閾値に設定するなど、代替の方法を用いてもよい。また、図４では８つの組み合わせによって該当ピクセルの塩分濃度を大・中・小の3段階に分けるような表を記載しているが、NDWI値から塩分濃度を推定する代替の表を用いてもよい。例えば各時期で閾値を2段階設けるようにしてもよいし、塩分濃度も大・小の2段階、あるいは大・やや大・中・小などの4段階以上に判別してもよい。例えば４段階に判別する場合ならば、津波発生後数ヵ月後の衛星画像におけるNDWI値の閾値としてT3₁及びT3₂（ただし、T3_1<T3₂）の2つを用い、図４において塩分濃度大と判定されているものを、さらにT3₁≦NDWI3n＜T3₂ならばやや大、NDWI3n≧T3₂ならば大と判定するなどの方法が考えられる。これらの閾値T3₁ 及びT3₂は、例えば、津波発生前の数時期の衛星画像から通常時のNDWI値を計算しておき、その平均から１標準偏差高い値をT3₁ 、２標準偏差高い値をT3₂とするなどの方法が考えられる。ここまでは、閾値によってNDWI値の高低を判断する手法を述べたが、各時期のNDWI値を閾値で分けずに、例えば津波発生前の数時期のNDWI値と時期iのNDWI値を用いてスミルノフ・グラブス検定を行い、時期iのNDWI値が外れ値に相当するかどうかを検定するような手法を用いても良い。いずれにしても各時期のNDWI値と塩分濃度との関係が表されていれば、その変換処理の方法については問わない。

以上で説明したように、３時期以上の衛星画像においてNDWIという一つの指標の変化の様子を継続的に調べることで、２時期の衛星画像からだけでは推定が困難な塩分濃度の把握ができる。例えば、津波発生前と津波発生後のような２時期の画像だけを用いて同様にNDWI値の変化を調べても、上記のような土地の排水性といった特徴は把握できず、従って塩分濃度のような特徴量は推定できない。さらに、この処理中に利用した図４に記載の表は、現地調査無しで作成できるために、地球上のあらゆる場所に適応可能であり、早い時期に情報提供することが可能になる。

以下ではさらに、ケースＡにおいて各時期のNDWI値から塩分濃度を推定する際に、3時期におけるNDWI値の変化パターンにおける特定の3パターンだけに着目して処理を行う例（以下これをケースＡ2と記載する）を説明する。少数のパターンだけに着目することにより、あまり重要でない部分の計算を省略できるため、処理速度の向上やメモリの節約につながる。さらに、例えば推定された塩分濃度を地図にマッピングして可視化を行う場合などは、着目している３パターンの変化が生じた部分にそれぞれR,G,Bの各色を割り当て、それ以外の点は白で示すなどとして色付けをすることで、重要な変化の生じた点だけを一見して把握することができ、視認性の向上にもつながる。以下では、塩分濃度推定の結果として各変化パターンに色割り当てを行い、地図上に可視化を行う例を説明する。

図５はケースＡ2の処理のフローチャートの例である。注目している衛星画像中にN個のピクセルが含まれているとすると、まずそれらに0からN-1までの番号を一意に振る。ここで、番号づけの順序に制約はなく任意に定めてよい。S201では注目対象のピクセル番号nを0に初期化する。次にS202で画像１と画像２のNDWIの差分Dn_12，画像２と画像３のNDWIの差分Dn_23を計算する。尚、ここでは変化パターンを抽出するために差分を用いているが、これを例えば比を用いて変化を表現するなどの代替手段を用いてもよい。次にS203において、これらの差分が図６に記載のいずれかの変化パターンに該当するか判定を行う。ここで、画像１と画像２の間での指標値の増加分に関するしきい値をTi12（＞0）、画像２と画像３の間での指標値の増加分に関するしきい値をTi23（＞0）、画像２と画像３の間での指標値の減少分に関するしきい値をTd23（＜0）として、
変化パターン１とは、Dn_12≧Ti12かつDn_23≧Ti23であるもの、
変化パターン２とは、Dn_12≧Ti12かつTd23≦Dn_23＜Ti23であるもの、
変化パターン３とは、Dn_12≧Ti12かつDn_23＜Td23であるものをさす。

いずれかの変化パターンに該当すればS204でその変化パターンに応じた色をそのピクセルに割り当てる。この際、色の割り当てに重複のないように設定する。例えば、255段階で表現したピクセルのRGB値がそれぞれr,g,bであったときにそれを(r,g,b)と表記するものとして、変化パターン１には(255,0,0)を、変化パターン２には(0,255,0)を、変化パターン３には(0,0,255)を割り当てる。いずれの変化パターンにも該当しない場合は、S205において、変化パターン１から３のいずれにも割り当てられていない色、例えば前記色割り当ての場合では(0,0,0)や(255,255,255)などを割り当てる。この際、変化の程度を表現するため、変化量に応じて、色合いに濃淡をつけるなどの工夫をしてもよい。例えばDn_12とDn_23の絶対値の和｜Dn_12｜＋α｜Dn_23｜（ただし変化パターン1及び３に相当する時α=1，変化パターン2に相当する時α=0とする）の値を基に、割り当てられた色の濃淡に差をつける。前記Dn_12とDn_23の絶対値の和を64から191の128段階に量子化した値をSDとして、変化パターン１には(SD,0,0)、変化パターン２には(0,SD,0)、変化パターン３には(0,0,SD)、を割り当てる。いずれの変化パターンにも該当しない場合は、S105において、変化パターン１から３のいずれにも割り当てられていない色、例えば前記色割り当ての場合では(0,0,0)や(255,255,255)などを割り当てる。S206で注目ピクセルが最後の番号まで到達しているかを判定し、まだ最後まで達していなければS207でnをn+1に更新し、S202へ戻る。最後のピクセルに達していればS208で処理を終了する。このような処理を施すことで、各変化パターンに該当する点だけを3時期画像から可視化することができる。

ここで、津波発生前から津波発生直後にかけてNDWI値が増加している変化パターン１、２、及び３に該当する領域は、津波の被害を受け、土壌水分量が高くなっている箇所であることが推定できる。このような土地は塩分濃度が高くなっている危険性があると考えられる。一方でどの変化パターンにも該当しない場所は、大きな津波被害を受けておらず、従って塩分濃度も高くないと推定できる。

津波被害を受けたであろうと推定できる領域の中で、変化パターン３に該当する領域は津波発生数ヵ月後には土壌水分量が低くなっていることから、土壌の排水性が高く、あまり大きな塩害被害を受けていないと考えられる。一方で変化パターン１や２に該当する領域では津波発生数カ月後も依然として土壌水分量が高くとどまっていることから、土壌の排水性が低く、従って多量の塩分の蓄積により大きな塩害被害を受けている可能性があると推定できる。特に変化パターン１に該当するような土地では、津波によって被った水がはけないばかりか、その後の降雨などによってもたらされた水分まではけずに溜まり続けていると考えられ、際立って排水性の悪い土地であることが推定される。このように単に２時期の指標値の差を見るだけでなく、図６に示すような特殊な変化パターンに着目して可視化を行うことによって、津波被害を受けた地域の中でもさらに塩分濃度が高くなっている可能性が高い場所とそうでない場所を推定することができる。

ここまで述べてきたケースＡ及びケースＡ2のような処理を施すことで、塩分濃度の高い地域とあまり高くない地域が定性的に把握できる。このような結果を利用して、実際に被害を受けた地域の調査に訪れる際にどの地域を重点的に調査すべきかを決定することができる。そしてこのような現地調査を行うことでより詳細にNDWI値と塩分濃度の関係をモデル化することができ、定量的に塩分濃度の推定が行える。そこで以下では、そのような現地調査の結果と3時期以上(ここではM時期とする)の衛星画像から、定量的な推定モデルを用いて塩分濃度推定を行うケースＢについて説明する。

多時期の衛星画像におけるNDWI値と塩分濃度の関係を定量的にモデル化する際には、様々なモデルが考えられるが、ここではその一例として重回帰分析などの線形回帰モデルを用いる。もちろん、これをPLS回帰分析やその他非線形なモデルに置き換えても何ら問題はない。重回帰分析とは、ピクセルnに対応する地点の塩分濃度の値をSn、時期i(i=1,2,…,M)の衛星画像におけるピクセルnのNDWI値をNDWIinとしたときに係数aj(j=0,1,…,M)を用いて，
Sn＝a0＋a1NDWI1n＋a2NDWI2n＋・・・＋aMNDWIMn ・・・（式１）
と表し、係数ajを求めるような分析手法である。ここで、対象地域のすべての点について同一の係数ajを用いるのではなく、例えば津波発生前の通常時の衛星画像から土地被覆分類を行っておき、その種類ごとに別の係数を用いるなどの工夫をすると、その土地の性質ごとにモデル化ができるため、より推定の精度を向上させられる。

係数ajを決定するには、時期ごとの衛星画像からNDWI値を計算するとともに、ケースＡ及びケースＡ2で説明した手法を用いて現地調査を行う地点を決定し、そこで計測された値からSnを求め、回帰分析を行えばよい。この際、ケースＡ及びケースＡ2で説明した方法によって求められる土地ごとの定性的な塩分濃度の情報を参考にしながら、前述のように土地被覆分類を行った上で分類ごとに調査を行うのがよい。このように一定数のサンプリング点で計測を行い、係数が決定できれば、同一地域中の他の同分類の土地についても良い精度で塩分濃度の推定ができる。

現地調査では例えば塩分濃度と相関性の高いことが知られているEC値を計測すればSnを求めることができる。その他に塩分濃度と相関性のある、土壌pH値、ESP(Exchangeable Sodium Percentage;交換ナトリウム濃度)などを計測するなどしても良く、調査内容と塩分濃度の求め方については制限をしない。

本実施例の最後に、ケースＡ及びケースＡ2で説明した方法により推定された定性的な塩分濃度をもとに、優先的に現地調査を行う場所を決定する現地調査計画立案支援システムについて説明する。

既に述べたように、津波によって被害を受けた地域の状況を把握する、ひいてはケースＢで述べたような定量的な塩分濃度推定モデルをたてるためには、現地調査が必要となる。ところが、そのような地域の中でランダムに選んだ地点を調査するだけでは、最悪の場合塩分濃度が低い地点ばかりを調査してしまうことになりかねず、十分なサンプル数を得るためには調査のために大きなコストがかかってしまう。

本実施例で説明したような、定性的な推定モデルを用いて、早期に塩分濃度のマップを作ると、目的に合わせて効果的に調査する点を決定することができる。例えば、ケースＢで説明したような定量的なモデルを作るために、塩分濃度の高い場所のデータも中くらいの場所のデータも低い場所のデータも同数ほど必要な場合は、前記マップから塩分濃度の高い場所・中くらいの場所・低い場所を同数程度選べばよい。あるいは、塩分濃度の高い場所を重点的に調査したい場合は、定性的に塩分濃度が高いと推定された場所の中から必要数だけ調査点を選択すればよい。

また、このような調査を行う場合には、土地被覆分類を行った上で、その種類を勘案したうえで調査する場所と数を決定するとさらに効果的である。例えば、種類ごとに同数ずつ調査する、あるいは、例えば例えば水田や畑など、特に注目したい対象がある場合はその種類の土地を重点的に調査するなどとする。

土地被覆分類は、その地域の被覆分類図などが入手できる場合はそれを利用しても良いし、津波が発生する前の通常期の衛星画像から分類を行ってもよい。分類のための方法は、最尤法などの教師付き分類を用いてもよいし、ISODATA法などの教師なし分類を用いてもよい。

実施例１では、塩分濃度推定システムについて説明したが、同様に衛星画像から計算されるある指標値を用いて対応する土地の状態を表す特徴量を推定するという枠組みを、これ以外の例にも適用することができる。そこで本実施例では、その一例として、鉱山開発における重金属漏れ事故後の、周辺の植生の被害度を推定するシステム（以下これを植生被害度推定システムと呼ぶ）について説明をする。

鉱山開発の際には、時として不慮の事故や自然災害などで有害物質が周辺に漏れ、それが環境破壊を起こすことがある。一方で重金属漏れによる土壌汚染は、一般的に肉眼で確認するのが困難であると言われている。

そこで本実施例では、衛星画像から鉱山周辺の植生の活性度を見ることにより、重金属漏れ事故が周辺の環境に与えた被害度合いを推定するシステムについて説明する。

注目指標値としてはNDVI(Normalized Difference Vegetation Index)を用いる。NDVIは植生の活性度合いを表す指標であり、具体的には近赤外光の輝度NIRと赤色光の輝度Rを用いて、NDVI=(NIR-R)/(NIR+R)と計算される値である。

重金属漏れ事故発生前後やその後のNDVI値から推定される植生活性度合いの変化を追うことで、事故の影響により植生活性度が低下した箇所、さらにその中でも数年のうちに回復したものや数年たっても回復しない箇所などの判別ができる。

本実施例を実現するための装置の例や機能構成図の例は、実施例１で説明した図１に記載の装置や図２に記載の機能構成図と同様のものを用いることができるため、この部分の詳細な説明は省略する。ただし、図２における指標値計算部202においては指標値としてNDVI値を計算するようにし、指標値→土地状態特徴量変換部203においてはNDVI値から植生被害度への変換を行うようにする。

３時期の衛星画像から定性的な推定モデルを用いて植生被害度推定を行う場合（以下これをケースＡ’と記載する）、例えば重金属漏れ事故発生前の衛星画像、重金属漏れ事故発生後の衛星画像、重金属漏れ事故発生後数年後の衛星画像１枚の合計３枚の衛星画像を用いる。ただし、植物の活性度が季節によって変動するため、NDVI値も一年ごとに周期的な変動を示すことが知られていることから、衛星画像を選択する際には、なるべく同季節に撮影されたものを選ぶのがよい。あるいは、季節変動要素を打ち消すために、例えば、NDVI値の地域での平均的な変化データを用いて正規化するなどの処理を行った上で植生被害度推定を行うのがよい。

ケースＡ’における基本的な処理の手順は、塩分濃度推定システムのケースＡの処理手順と同様であるため詳細な説明は省略する。図３に示すフローチャートにおいて、NDWIをNDVIに、塩分濃度を植生被害度に置き換えればよい。ただし、S103において参照する表は、例えば図７に記載のもののように、このケースに適合するものに置き換えなければならない。ここで、図７では、重金属漏れ事故発生前はNDVI値が高かったものの重金属漏れ事故発生後にNDVI値が小さくなり、数年後もなおNDVI値が低いままであるような土地の植生被害度を大としている。これは、このような土地は、重金属漏れ事故発生前後でNDVI値が下がっていることから、確かに重金属漏れによる被害を受け、なおかつ、数年後もまだNDVI値が回復しないことから、特にその土地に多量の重金属が漏れ出て被害が大きい、あるいはその土地の植物が重金属汚染への耐性が弱く結果として被害が大きい、と考えられるためである。また、重金属漏れ事故発生前はNDVI値が高かったものの重金属漏れ事故発生直後にNDVI値が小さくなり、数年後には再びNDVI値が上がっているような土地の植生被害度を中としている。これは、このような土地は、重金属漏れ事故発生前後でNDVI値が下がっていることから、確かに重金属漏れの被害を受けているが、数年後にはNDVI値が上がっていることから、比較的漏れ出た重金属の量が少ない、あるいは植物の重金属汚染への耐性が高いといった理由で、被害はある程度抑えられていると考えられるからである。また、重金属漏れ事故発生前はNDVI値が高く、重金属漏れ事故発生直後もNDVI値は高いが、数年後にはNDVI値が低くなっているような土地の植生被害度も中としている。これは、このような土地は、重金属漏れ事故発生前後ではNDVI値に変化はないものの、数年後にはNDVI値が下がっていることから、急速に被害が広がるようなことはなかったものの、一定程度の影響があったと考えられるためである。その他、重金属漏れ事故発生直後でもNDVI値が高かったような土地は、そもそも重金属漏れ事故の被害を受けていない、あるいは植物の重金属汚染への耐性が十分高いものと考えられるので植生被害度を小としている。また、重金属漏れ事故発生直後のNDVI値は低いが重金属漏れ事故発生前のNDVI値もまた低いような土地は、もともと植物が生えていないため、植生への被害度は小としている。

このような表は、現地調査に行くこと無しに作成ができるために、植生のあるところならば地球上のあらゆる場所に適応可能である。

また一般に、植物の反射スペクトルは樹種によって異なり、それに伴ってNDVIの平均的な値も樹種によって異なるので、事故発生前の通常時期の衛星画像を用いて樹種分類を行ったうえで、樹種ごとに別々の参照表を用いると精度が向上し好ましい。

塩分濃度推定システムにおけるケースＡ2で説明した方法のように、3時期におけるNDVI値の変化パターンにおける特定の数パターンだけに着目して処理を行う例（以下これをケースＡ2’と記載する）についても同様に考えることができる。ただしここでは、例えば図８に示すような4つの変化パターンに着目する。変化パターン４に該当する場所は、事故後数年たってもNDVI値が下がり続けていることから、最も深刻な被害のある場所であると推定できる。変化パターン５に該当する場所は、事故後NDVI値が下がり、その後も回復していないことから大きな被害を受けている場所であると推定できる。変化パターン６に該当する場所は、事故後NDVI値が下がっているものの、数年後にはNDVI値が回復していることから、事故の影響はあったものの比較的被害は軽微であった場所であると推定できる。変化パターン７に該当する場所は、すぐには事故の影響が現れなかったものの、その後NDVI値が低下していることから、一定の被害があった場所であると推定できる。このようにして、衛星画像を用いてNDVI値から被害度合いを推定することにより、広域的に被害を把握することができ、例えば汚染対策を行う場所の優先度付けをしたり、対策に必要な費用をおおまかに見積もったりすることが可能となる。

最後に、多時期の衛星画像におけるNDVI値と植生被害度の関係を定量的にモデル化する方法についても説明する。この場合、塩分濃度推定システムにおけるケースＢで行った処理と同様に、例えば線形回帰モデルなどが植生被害度推定システムにも適用できる。モデルの係数推定のために、植生被害度をNDVI値とは別の手段を用いて測定しなければならない。これには、例えばケースＡ’あるいはケースＡ2’の処理結果を利用して、現地調査を行う地点を決定し、土壌をサンプリングして重金属がどの程度含まれているかを調べるなどの方法を用いると良い。

以上のように、例えば実施例１に記載の塩分濃度推定システム及び本実施例において説明した植生被害度推定システムのような形で、請求項１に記載の土地状態推定システムが実現できる。

100 計算機
101 バス
102 入力制御部
103 中央演算処理部
104 主記憶部
105 補助記憶部
106 出力制御部
110 入力装置
120 出力装置

Claims

衛星画像データを格納する格納部を有する計算機における土地状態推定方法であって、
前記格納部に格納された複数の衛星画像データから、所定の土地領域の津波発生前の第１時期の第１の衛星画像データ、津波発生直後の第２時期の第２の衛星画像データ、及び津波発生から時間が経過した第３時期の第３の衛星画像データを抽出し、
前記第１の衛星画像データ、第２の衛星画像データ、及び第３の衛星画像データから、それぞれ画像のピクセル毎に土壌水分量指標値を算出し、
第１の衛星画像データから求めた土壌水分量指標値と第２の衛星画像データから求めた土壌水分量指標値との差分もしくは比率を示す第１の変化指標と、第２の衛星画像データから求めた土壌水分量指標値と第３の衛星画像データから求めた土壌水分量指標値との差分もしくは比率を示す第２の変化指標とを算出し、
画像のピクセル毎に、第１の変化指標が正の第１の閾値より大であり、かつ第２の変化指標が正の第２の閾値より大である第１のパターン、第１の変化指標が前記正の第１の閾値より大であり、かつ第２の変化指標が負の第３の閾値と前記正の第２の閾値の間である第２のパターン、及び第１の変化指標が前記正の第１の閾値より大であり、かつ第２の変化指標が前記負の第３の閾値より小である第３のパターンのいずれに当るか判別し、
判別されたそれぞれのパターンごとに画像を色分けして表示することを特徴とする土地状態推定方法。