JP4003869B2 - 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 - Google Patents
地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 Download PDFInfo
- Publication number
- JP4003869B2 JP4003869B2 JP2002103966A JP2002103966A JP4003869B2 JP 4003869 B2 JP4003869 B2 JP 4003869B2 JP 2002103966 A JP2002103966 A JP 2002103966A JP 2002103966 A JP2002103966 A JP 2002103966A JP 4003869 B2 JP4003869 B2 JP 4003869B2
- Authority
- JP
- Japan
- Prior art keywords
- data
- observation
- time
- coefficient
- observation data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000003672 processing method Methods 0.000 title description 8
- 238000012545 processing Methods 0.000 claims description 114
- 238000000034 method Methods 0.000 claims description 71
- 230000006870 function Effects 0.000 claims description 60
- 230000001932 seasonal effect Effects 0.000 claims description 32
- 230000008859 change Effects 0.000 claims description 20
- 230000000737 periodic effect Effects 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 20
- 230000007774 longterm Effects 0.000 claims description 6
- 230000035945 sensitivity Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 description 17
- 238000010586 diagram Methods 0.000 description 15
- 238000001514 detection method Methods 0.000 description 12
- 238000004422 calculation algorithm Methods 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 230000002265 prevention Effects 0.000 description 9
- 238000001914 filtration Methods 0.000 description 7
- 238000004088 simulation Methods 0.000 description 7
- 238000013480 data collection Methods 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 6
- 238000003860 storage Methods 0.000 description 6
- 230000002123 temporal effect Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 239000000284 extract Substances 0.000 description 5
- DFPOZTRSOAQFIK-UHFFFAOYSA-N S,S-dimethyl-beta-propiothetin Chemical compound C[S+](C)CCC([O-])=O DFPOZTRSOAQFIK-UHFFFAOYSA-N 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 4
- 238000009826 distribution Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 4
- 238000003786 synthesis reaction Methods 0.000 description 4
- 238000012731 temporal analysis Methods 0.000 description 4
- 238000000700 time series analysis Methods 0.000 description 4
- 238000012937 correction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 239000002131 composite material Substances 0.000 description 2
- 238000001035 drying Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000005295 random walk Methods 0.000 description 2
- 238000012502 risk assessment Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- LQEOTHINZCAFAC-UHFFFAOYSA-N (4-methylpiperazin-1-yl)-(3,4,5-trimethoxyphenyl)methanone Chemical compound COC1=C(OC)C(OC)=CC(C(=O)N2CCN(C)CC2)=C1 LQEOTHINZCAFAC-UHFFFAOYSA-N 0.000 description 1
- 239000002028 Biomass Substances 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 238000006424 Flood reaction Methods 0.000 description 1
- 229910000831 Steel Inorganic materials 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 239000000779 smoke Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 239000010959 steel Substances 0.000 description 1
Images
Landscapes
- Image Analysis (AREA)
- Image Processing (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体に係り、特に、高頻度観測衛星データのノイズ除去処理による植生の季節変動観測法に関する。
【従来の技術】
【０００２】
１９９７年から１９９８年にかけてインドネシアで大規模な森林火災が発生し、煙害等で周辺諸国も含めて甚大な影響を及ぼした（文献１，２参照。なお、文献は後述する。以下同様。）。その後も東南アジア地域では森林火災が頻発し、大規模火災の再発生が懸念されている。現在、国際協力事業団（ＪＩＣＡ）の「インドネシア森林火災予防計画プロジェクト」の第二フェーズが進行中であり、その中で、「森林火災早期発見システム」が稼働している。これは、火災発見を衛星データを用いて行い、森林管理を担当する機関に通報するものである(文献２参照）。一方、火災予防の観点からは火災危険度の評価が待たれている。
【０００３】
火災危険度評価の基礎となるのは森林植生図であるが、乾燥に対する森林の感受性を評価する必要があることを考えると、単なる森林区分では不十分で、時間を含んだ情報に基づいて危険度を算出することが期待されている。
【発明が解決しようとする課題】
【０００４】
従来の植生の季節変動観測法のための時系列処理手法の開発に対して、その課題について以下に説明する。
【０００５】
従来のＮＯＡＡ／ＡＶＨＲＲ あるいはＳＰＯＴ／ｖｅｇｅｔａｔｉｏｎ、ＭＯＤＩＳなどの高頻度観測衛星からのデータによって広域での植生観測が行われている(文献３−１３参照)。特に、大規模森林火災や乾燥害などの農業・林業災害は、植生の乾燥の程度に影響されるため、植生の状態を迅速かつ的確に把握する手法の開発は防災面で大きな意味を持っている。しかし、光学センサーでは、雲やヘイズの影響が避けられず、的確な植生分布およびフェノロジーの把握が困難である。特に、熱帯を対象とする場合、数ヶ月間にもわたって雲の影響があり、雨期にはほとんど有効なデータが観測できないこともある。
【０００６】
また、従来より広域での植生分布についても数多くの研究がある(文献３−７参照)。例えば、ＧＬＣＣ（Ｇｌｏｂａｌ Ｌａｎｄ Ｃｏｖｅｒ Ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ ｄａｔａ ｂａｓｅ）があるが、日本にサバンナが出現するなど誤分類が見受けられる(文献７参照)。これも雲の影響を十分に取り除くことができなかったのが原因である。
【０００７】
さらに、これまでにも、高頻度観測衛星のデータから雲の影響を取り除いてＮＤＶＩ（植生指数）の季節変化を抽出する試みがなされている(文献８参照)。例えば、解像度を落として、複数の画素の最大値を取るものや、１０日〜１ヶ月のデータの最大値を取るものなどがある。しかし、これらは、解像度が低かったり、データの取得間隔が一ヶ月であったりと即時性の要求される防災データとしては使いにくい。データの合成期間が長くなれば、時間的な分解能が低下する。また、これらは、植生の変動を一週間から１０日の間隔で把握する必要のある農林業分野では応用できないという課題がある。
【０００８】
また、Ｆｏｕｒｉｅｒ展開類似の方法を用いて、植生指数の季節変化を三角関数で表すこともなども行われている(文献９−１４参照)。この方法にはどの関数を使うかについて恣意性があり、様々な土地被覆を含む広域データの処理には適していない。また、この方法では、関数の当てはめにおいて、処理に用いる三角関数群はあらかじめ決められており、必ずしも処理に適した関数群が選ばれていないという課題がある。
【０００９】
また、１６×１６などのウインドウを設定し、その中での最大値をウインドウ内のすべての画素値に置き換える方法もある。しかし、このような方法では、空間分解能が低下するという課題がある。
【００１０】
本発明は、以上の点に鑑み、空間分解能を保持したまま、任意の時間分解能を実現させながら、雲やヘイズの影響を取り除くために、時系列モデルを用いた衛星データ処理手法ＬＭＦ（Ｌｏｃａｌ Ｍａｘｉｍｕｍ Ｆｉｔｔｉｎｇ：局所最大値フィッティング）を提案することを目的とする。すなわち、本発明は、画素毎の時系列的な特徴を保存し、原データの時間分解能を完全に保持したまま、雲などの影響およびシステムノイズの除去を目的とするものである。
【００１１】
また、本発明は、時間分解能を保持したまま植生などの変動の様子が捉えられるようにすることで、防災、環境、農業分野での広域現状把握、および予測が容易とすることを目的とする。さらに、本発明は、適切な関数を自動的に選択し、ノイズ除去処理を自動で行えるようにすることにより、定常業務の一部分とすることを可能とすることを目的とする。さらに本発明は、水分指数など他の地球高頻度観測衛星のデータにもそのまま転用および応用を可能とすることを目的とする。
【００１２】
【課題を解決するための手段】
本発明の特徴のひとつとしては、例えば、処理対象画素の時系列データについて、前後ｎ時期について、局所最大値を抽出するフィルターを作用させる。そこで得られた局所最大値群に関数当てはめを行い、画素値を推定する。ある時期に、推定値と実際の画素値のずれが、設定したしきい値よりも大きい場合には雲やシステムノイズの影響とみなしてその時期のデータは採用しない。採用されたデータだけを用いて再度、局所最大値の算出と関数当てはめを行った値で画素値を置き換える。必要であればこの処理を繰り返し、結果が収束するまで行う。
【００１３】
本発明の他の特徴としては、例えば、画素毎の時系列的な局所最大値を求め、それらに対して関数の当てはめをくりかえし行うことで、有効なデータだけを用いた推定を行うものである。上記の処理を画像全体に対して行うため、雲やシステムノイズの影響が除去されるという作用を有する。各画素毎の処理は独立しているため、原データの空間分解能は完全に保持される。
【００１４】
また、本発明は、時期ｋ近傍の局所的最大値を、時間遅れを少なくするように抽出する作用を有する。これにより、雲などの影響を受けた時期の画素値が除去されるという作用を有する。
【００１５】
本発明の第１の解決手段によると、
処理部は、観測対象エリア内の各画素に対応して、衛星による植生の季節変動についての時系列の第１の観測データを記憶した第１の観測データファイルから、画素毎の第１の観測データを読み込むステップと、
【００１６】
処理部は、各画素毎に、ある時期又は時点に対応して、該時期又は時点の前後所定期間の局所最大値を抽出することにより時系列の第２の観測データを求め、それを各画素に対応して第２の観測データファイルに記憶するステップと、
【００１７】
処理部は、第２の観測データファイルに記憶された第２の観測データを、長期的非周期項及び季節変化項を含む次式の関数を用い、各第１の係数ｃｉを最小二乗法又は他の推定法により決定することでフィッティングし、各第１の係数ｃｉを第１のパラメータファイルに記憶するステップと、
【００１８】
【数１０】
［ ここで、
ｃｉ：求める第１の係数
ｔ：時期又は時間
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列（例：１，２，３，４，６，１２）
Ｍ：１年間のデータ数（シーン数） ］
【００１９】
処理部は、第１の係数ｃｉで定められた前記関数に従い、第２の観測データから雑音データを除去した第１の推定観測データを求め、それを第１の推定観測データファイルに記憶する、又は出力部に出力、若しくは表示するステップと
を含む地球観測衛星データのノイズ除去方法、これらのステップをコンピュータに実行させるためのノイズ除去処理プログラム、及びノイズ除去処理プログラムを記録した記録媒体が提供される。
【００２０】
【発明の実施の形態】
本実施の形態では、植生の季節変動観測のために、時系列モデルを用いた新しい衛星データ処理法を提案する。また、雲の影響を取り除いた植生の季節変動に関する観測データから、火災危険度を算出した。以下、火災危険度評価の基礎となる植生の季節変動観測のための衛星データ処理法について述べる。なお、植生の季節変動観測データとしては、例えば、植生指数データ、水分指数、中間赤外線バンド、熱赤外線輝度バンド、温度データ等である。また、観測データに含まれる可能性がある雑音データは、例えば、雲、霧、雨、ヘイズ等の大気状態に起因するもの、センサ感度等のシステムに起因するもの等である。
【００２１】
１． ＬＭＦの理論
ＬＭＦ処理は、時系列フィルタリングと関数によるフィッティング手法を組み合わせたアルゴリズムで、対象画素の時系列データから、雲やヘイズなどの影響を除去し、季節変化パターンを抽出するものである。フィッティングに用いる関数を、（１）式に示す。
【００２２】
【数１１】
ここで、最初の二項は長期的非周期項、その他の項は季節変化項である。また、各パラメータは、以下の通りである。
ｃｉ：求めるパラメータ：ｃ０は平均、ｃ１はトレンド係数、それ以外は周期関数のパラメータとなる。
ｔ：利用するデータの時間間隔に基づく時間単位数（１０日合成のＮＯＡＡデータを用いる場合なら、１０日が１単位）
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：このソフトウェアの構造は周期関数の種類には依存しないため、任意の周期関数と置き換えて利用することが可能。
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列：１，２，３，４，６，１２。これによって、１ヶ月周期、２ヶ月周期、３ヶ月周期、４ヶ月周期、６ヶ月周期および１２ヶ月周期の周期関数（ｓｉｎ、ｃｏｓ）が当てはめられる。
Ｍ：１年間のデータ数：１０日間最大値合成のＮＯＡＡデータを用いる場合、１年間３６シーンになるため、Ｍ＝３６である。
なお、この場合、シーン数Ｎと１年間のデータ数Ｍが一致しないので、一般的な離散 Ｆｏｕｒｉｅｒ級数は計算できない。
【００２３】
１．１． シミュレーションによる処理特性の評価(文献１５参照)
植生指数や熱赤外線輝度温度データでは、観測データには雲などの影響がすでに含まれており、真値は不明である。そこで、本実施の形態では、シミュレーション的手法を用いてフィルタ特性の評価を行った。
まず、図１に、シミュレーションによるフィルタ特性の評価図を示す。例えば、（１）式の係数の組ｃに、発生した乱数（Ｓ１１）を代入して、時系列データを１００００組発生させる（Ｓ１２）。解析のパラメータは、一例として、Ｎ＝５４、Ｍ＝３６とした。次に、雲やヘイズの影響を受けた状況をシミュレートするため、ノイズを発生し（Ｓ１３）、そのノイズを元の値の１０％〜７０％のものに置き換えて、ノイズが加わった状況をシミュレートした（Ｓ１４）。ノイズが加算される確率は各時期について０．５とした。そして、ノイズ加算後の時系列データに対して、ノイズ除去処理し（Ｓ１５）、その結果を調べた（Ｓ１６）。結果の解析は発生させた係数の組ｃと、処理後に再生された係数の組ｃ’とについて、単回帰分析を行った。
【００２４】
１．２． 最小二乗法によるフィッティング
一般に、ノイズ加算後のデータをそのまま（１）式でフィッティングした場合、どの項も係数の約７割が保存される。これに、（２）式で表されるような局所最大値を抽出するような時系列フィルタを作用させる。なお、この局所最大値の計算式は一例であり、適宜の式を用いることもできる。
【００２５】
【数１２】
ここで、
ｄｔ’ は時間単位ｔの時の局所最大値 （第２の観測データ）
ｄｔ は時間単位ｔの時のデータ値 （第１の観測データ）
ｗ は局所最大値を求める片側時間単位数（ウインドウサイズ）
である。
【００２６】
（２）式では、ｄｔにおける左右ｗ間のデータの最大値をもとめ、その小さい方を局所最大値としている。この考え方が本発明のひとつの特徴である。本実施の形態では高頻度観測衛星からの観測データｄｔ（第１の観測データ）に対して、（２）式のようなフィルタリング処理を実行して局所最大値ｄｔ’（第２の観測データ）を求め、さらに、（１）式の関数で最小二乗法によりフィッティングして各係数ｃｉを求める。これによって、データに影響を与える雲やシステムノイズを軽減して、季節変化する地表の情報として時間単位ｔに近い時期のデータを自動抽出することができる。
【００２７】
１．３． フィッティング関数の自動選択
上述のような式（２）によるフィルタと最小二乗法によるフィッティングを併用したものでは、定数項、一年周期、半年周期項はほぼ回復できるが、それら以外の短周期関数項の係数については保存されない場合がある。これを改善するため、本発明は、さらに以下のように、（１）式の関数のなかで、必要なもの以外はフィッティングに用いないようにしてもよい。なお、最小二乗近似に対しても、適宜のフィッティング方法を用いることができる。
【００２８】
いま、最小二乗近似を用いて時系列データを関数でフィッティングすることを考える。フィッティングに用いる関数の個数（種類）が多ければ残差は減少するが、フィッティングした結果は不安定（結果がデータに大きく依存する）となる。赤池らによると、赤池の情報量基準（ＡＩＣ）が最小となるようなモデルが最適であるとされている(文献１６参照)。残差が正規分布に従うと仮定されているから、残差の分散をσとするとＡＩＣは、対数尤度とパラメータ数の差で表される（３）式のようになる。
【００２９】
【数１３】
となる。ただし、
N は所定期間のデータ数
j は用いる関数の個数（周期関数のペア数）
σはオリジナルデータとフィッティングによるデータの残差の分散
本実施の形態ではAICの値が最小になるように、周期関数の組み合わせを自動的に選択するようにした。即ち、この値が最小になるようにjを選ぶようにした。なお、正の定数倍は結果を変えないため問題ない。この例においては、シミュレーションの結果、三ヶ月周期より短い周期の項については改善がみられたが、逆に四ヶ月周期の項についての相関係数は低くなる。
【００３０】
図２は、シミュレーション結果の例を示す図である。
１．４． フィッティング処理の繰り返しによる高精度化
本実施の形態ではさらに、処理の高精度化を実現するため、フィッティングを繰り返すアルゴリズムを採用した。最初に、原データに局所最大値抽出フィルタを作用させた後、フィッティングする。その結果（計算値）と、原データとの差をとり、しきい値よりも大きなものを取り除く。データレンジの１５〜２０％をしきい値とすれば、１年および半年に相当する周期を持つ関数の係数は９５％以上、それ以外の短周期成分も８０％が回復できることが明らかになった(文献１５参照)。
【００３１】
１．５． 熱赤外バンドデータの処理
ＮＯＡＡ／ＡＶＨＲＲの熱赤外バンドでは、雲のあるピクセルは温度が低くなる性質を用いれば、ＮＤＶＩと同じ計算方法を用いればよいと考えられる。したがって、本発明は輝度温度についても処理が可能である。
【００３２】
２． カルマンフィルタを用いたＬＭＦアルゴリズムの高精度化
２．１． ＬＭＦ処理の問題点
ＬＭＦは高頻度観測衛星からの植生指数や輝度温度データから雲の影響を軽減し、季節変化パターンを抽出するのに有効な手法である。しかし、処理期間すべてにわたってフィッティング処理を行うため、どうしても平均化された結果が得られてしまう。また、データ収集期間が短い（１年程度）場合には、フィッティングの結果が不安定になりやすく、年々変動の抽出は困難な場合があった。そこで、本実施の形態ではさらに、データ処理に状態空間モデルを採用し、データ処理を行う手法を開発した（文献２２参照）。
【００３３】
また、データを防災に役立てようとするならば、単に現況を把握するだけでなく、短期的にせよ状況の予測をしながら対策を講じなければならない。したがって、時系列データの解析についても、予測が可能な手法を採用する必要がある。状態空間モデルであれば、この点も問題なく適用することができる。本実施の形態では、一例として、状態空間モデルのパラメータの推定には、Ｋａｌｍａｎの提唱した手法であるかルマンフィルタ（Ｋａｌｍａｎ Ｆｉｌｔｅｒｉｎｇ）を応用するが、それに限らず、適宜の推定法を用いることができる。
【００３４】
カルマンフィルタは、人工衛星の軌道・姿勢の予測と制御のように物理的に定式化が可能なものだけでなく、降水量からの洪水予測や経済状況の時間的変化の予想のように、直接に原因と結果の関係が定式化できなくても統計量さえあればモデル構築が可能な点が大きな特徴である(文献２２参照)。さらに、
○ 少数の観測値しか得られなくても複雑な数値モデルを構築できる。
○ 最新の観測値によってパラメータ群を更新・最適化するので、周期的でない要因も考慮できる。
などの優れた特徴を有しているので、エル・ニーニョのように非定常的な擾乱要素の影響も含めたモデルの構築が可能である(文献２２参照)。
【００３５】
２．２． 時系列モデルの設定
高頻度観測衛星データを用いて得られるＮＤＶＩの季節変化プロファイルは、三角関数の組み合わせで表現できることがＬＭＦ処理の結果から明らかとなった。本実施の形態では、この知見を生かし、離散時変スペクトルの推定技術を応用した(文献２３参照)。すなわち、（１）式の係数ｃｋの組を状態空間モデルの状態ベクトルと考え、その値の推定（状態推定）をカルマンフィルタによって行う。
【００３６】
ＬＭＦでは、雲やヘイズの影響のない季節変化プロファイルを、（１）式の形で記述できるとした。ここで、右辺の第二項が非周期的変動を考慮した項である。この式での各係数ｃｋは、（１）式では固定であるが、これを時変として、季節変化プロファイルの変動に対応できるようにした。右辺第一項が時変になったことで、トレンドも表現できることから、第二項は削除できる。したがって、（４）式となる。
【００３７】
【数１４】
ここで、各パラメータは（１）式と同じ、ただしｃ１（トレンド項）は無い。
この式は先のLMFに対して、パラメータｃが全て時変変数としたものである。これによって、衛星が代わり、観測センサの感度が変わるなどのシステム変動（システム雑音の一部）があっても、季節変化を追うことができる。この点が、特に、このソフトウェアの画期的なところである。
（４）式のｃｋ（ｔ）の時間的変化は自己回帰モデルで現される。
【００３８】
【数１５】
ただし、ａｋは自己回帰係数、ｗＫは平均値０、分散がσｋの白色雑音（システム雑音）、ｍは自己回帰モデルの次数。
ここでは、Ｋ＝０のときのモデル、
【００３９】
【数１６】
は、酔歩散乱モデル
【００４０】
【数１７】
と、時刻ｔと時刻ｔ−１の変化率が等しいと仮定したモデル、
【００４１】
【数１８】
を折衷したものを採用した。
酔歩散乱モデルだけだとデータの追従性が低下し、年々変動を表すことができない場合がある。また、（７−２）式のモデルだけを採用すると、フィルタリング結果の値が大きく変動しすぎる場合がある。（６−１）式を採用することによって、センサ感度の違いなどのノイズによるデータの急激な変化を押さえられる。
K＝０以外の時のモデルは（６−２）式である。
【００４２】
【数１９】
ｗｋ（ｔ）をシステム雑音とみなすと、観測される値ｙ（ｔ）は、（８）式の状態空間モデルで表すことができる。
【００４３】
【数２０】
ここで、各パラメータは以下の通りである。
ｃは求めるパラメータ
Ｆは予め定められた行列
Λはシステム雑音が、一度しか作用しないようにするための定数行列。
Ｈは周期関数を表す行列
ｖ（ｔ）は平均０、分散がσvの白色雑音（雲、ヘイズ等の大気の状態な等に起因する観測雑音）
【００４４】
【数２１】
この点も本発明の汎用性を高める特徴である。これによって、例えば、異なる衛星センサを利用せざるを得ない20年間の衛星データから環境変化を一度に分析できるようになった。
また、各状態変数の雑音の分散Ｑ
【００４５】
【数２２】
は真値が不明であるが、この値が計算上必要であるので、このプログラムでは処理パラメータとして与えることにする。
三角関数による時変係数推定問題では、複素形式での表記が一般的であり、計算量も半減できるが(文献２３参照)、この例では、使用する関数を限定しないアルゴリズムを開発するため、敢えて複素形式にはしなかった。
【００４６】
実際にはＮＤＶＩデータや熱赤外領域の輝度温度データでは、雲の影響があれば、値が必ず低くなる性質がある。これは、観測雑音ｖ（ｔ）が有色となる可能性を含んでいる。したがって、この状態空間モデルで求められる状態変数の値は、必ずしも最適な量とは言えずに、準最適なものである可能性があるが、データ処理には時系列局所最大値フィルタを併用しているので、観測雑音の有色性は軽減されているものと考える。
【００４７】
２．３． カルマンフィルタ
次に、時変係数ｃ（ｔ）を、カルマンフィルタを用いて推定することについて説明する。ここで、ｃ（ｔ｜ｔ）を時刻ｔまでの観測値がわかっている場合の値とすると、一期先の予測はｃ（ｔ＋１｜ｔ）となる。
カルマンフィルタは（９）式の通りである。観測値をｙ（ｔ）とする。
【００４８】
【数２３】
ここで、Ｖ（ｔ）は時刻ｔにおけるｃ（ｔ）の共分散行列である。
これらの式を用いれば、Ｈ、Ｆ、Ｑのすべての要素があらかじめ定められているとすれば、ｃ（ｔ｜ｔ−１）やＶ（ｔ｜ｔ−１）が既知のとき、観測値ｙ（ｔ）から他の量を計算することができる。また、ｙ（ｔ）が何らかの理由で欠損になった場合には、カルマンゲインが計算できないので、（１０）式を用いて状態更新を行う。
【００４９】
【数２４】
ここで、システム雑音の分散σｋおよび、観測雑音の分散σｖは未知なので、処理パラメータとみなす。初期値として、
【００５０】
【数２５】
（ただし、ε＞０）
ε０は正の数（０〜１）、Iは単位行列である。
を用いて最初の計算を行うことがしばしば行われる。
ところで、カルマンフィルタのアルゴリズムを用いて、時間をｔ＝Ｔからｔ＝１に時間を遡って処理すればデータの平滑化（ｓｍｏｏｔｈｉｎｇ）も可能である。それらの式は、
【００５１】
【数２６】
である。時期１からＴまでのデータが取得できており、なおかつ、カルマンフィルタにより、ｃ（ｔ｜ｔ）、Ｖ（ｔ｜ｔ）およびＶ（ｔ＋１｜ｔ）などの要素が計算できているとする。（１２）式を、ｔ＝Ｔからｔを１ずつ減じてゆけば、ｔ＝１まで求めることができる。ｃ（ｔ｜Ｔ）は、得られたデータすべてを用いて推定しているので、一般に、ｃ（ｔ｜ｔ）よりも精確な値であることが知られている。
ここで述べた手法をＬＭＦ−ＫＦ（Ｌｏｃａｌ Ｍａｘｉｍｕｍ Ｆｉｔｔｉｎｇ ｗｉｔｈ Ｋａｌｍａｎ Ｆｉｌｔｅｒ ａｌｇｏｒｉｔｈｍ）と以後略記する。
この方法で、ベクトルＱの要素をすべて１０-4程度にすれば、事実上、各状態変数ｃｋは変化しなくなるので、従来のＬＭＦの結果を再現することも可能である。計算量は増大するが、ＬＭＦ−ＫＦはＬＭＦを完全に包含した手法であることがわかる。
また、このＬＭＦ−ＫＦ法の特徴として、各係数についての自己回帰モデルである（５）式の次数ｍを上げれば、長期間データ欠損があってもそれらを推定することが可能であり、また、そのままのアルゴリズムでデータ収集期間の終わりからの数年間先の予測までが可能である。
【００５２】
２．４． 実際のデータ処理手法
画素毎に行われる、実際のデータ処理の流れを以下に説明する。先にも触れたように、カルマンフィルタルゴリズムの適用にあたって、未知量であるシステム雑音の分散σＫおよび観測雑音の分散σｖはパラメータとして与えることにするが、ｃ（１｜０）、Ｖ（１｜０）も必要である。これらの要素を画素毎に推定するのは大変なので、最初に（１）式を適用してＬＭＦにより係数を求めて初期値とする。
【００５３】
【数２７】
（ただし、ε＞０）
ε0は正の数（０〜１）、Iは単位行列である。
本実施の形態では、これらの初期値を用いてフィルタリングと平滑化を一回行って、次回の計算の初期値とする。
【００５４】
【数２８】
その後、フィルタリングと平滑化を繰り返し、得られたｙ（ｔ｜Ｔ）がほとんど変化しなくなったところで値を確定する。（ｔ ＝ １，２，３，…． ，Ｍ（最大シーン数））
このアルゴリズム（ＬＭＦ−ＫＦ）は、ＬＭＦと同様に、画素毎に処理が行われるため、ＬＭＦの時系列解析のサブルーチンを入れ替えるだけでよい。もちろん、プログラム内でのデータ構造も同じなので、並列化が容易である。
【００５５】
３．ハードウェア
図３に、命令列スケジュール処理を実行するための装置の構成図を示す。この装置は、処理部（ＣＰＵ）１、入力部２、出力部３、記憶部４を備える。処理部１は、ＣＰＵであり、命令列スケジュールプログラムによりスケジュールを実行する。記憶部４は、計算パラメータファイル４１、第１の観測データファイル４２、第２の観測データファイルファイル４３、第１の推定観測データァイル４４、第２の推定観測データファイル４５、第１のパラメータファイル４６、第２のパラメータファイル４７を有する。
【００５６】
４．フローチャート
図４に、実際のデータ処理のフローチャートを示す。なお、破線枠内の処理は、画素毎に行われる。
処理部１は、入力部２又は記憶部４の計算パラメータファイル４１から計算パラメータを読み込む（Ｓ１０１）。ここでは例えば、計算パラメータとして、画像数、像ライン数、画像コラム数、関数セット種類数、欠損シーン数、カルマンフィルタ使用時の関数セット数、局所最大値算出時の片側ウインドウサイズ雲検出閾値、カルマンフィルタ処理時の最大繰り返し回数、カルマンフィルタ収束判定値、システムノイズの分散パラメータなどがある。
【００５７】
処理部１は、衛星による画素データ毎の植生の季節変動についての時系列な観測データを記憶した記憶部４の第１の観測データファイル４２から、画素データ毎の第１の観測データを読み込む（Ｓ１０３）。この観測データは、例えば、植生指数データ、水分指数、中間赤外線バンド、熱赤外線輝度バンド、温度データ等がある。処理部１は、記憶部４の第１の観測データファイル４２から、１ライン又は１カラムのデータを読み込み、読み込まれたライン又はカラムのデータを並列演算計算機のパラメータに設定してもよい。
【００５８】
処理部１は、式（２）を用いて各画素毎に、ある時期の前後所定期間の局所最大値を抽出し、第２の観測データを得て、第２の観測データファイル４３に画素毎に記憶する（Ｓ１０５）。
【００５９】
処理部１は、第２の観測データファイル４３から第２の観測データを読み出し、これに基づき式（１）の関数にフィッティング処理を実行して各係数ｃｋを最小二乗法により決定し、その係数を第１パラメータファイル４６に記憶する（Ｓ１０７）。なお、処理部１は、式（３）で示される赤池の情報量基準ＡＩＣが最小となる周期関数の組み合わせを選択するステップをさらに含むことができる。次のような各ステップが実行される。即ち、
・処理部１は、全ての関数の組み合わせを発生させる第１ステップと、
・処理部１は、各組み合わせに対して最小二乗法を実行してフィッティングに用いられる係数の関数を求める第２ステップと、
・処理部１は、最小二乗法により求めた結果について、それぞれＡＩＣを求める第３ステップと、
・処理部１は、ＡＩＣの最小の組み合わせを使用するフィッティングに用いられる関数として選択する第４ステップと
を含む。
【００６０】
また、処理部１は例えば、ステップＳ１０５の前又は後若しくはステップＳ１０７の前又は後等に、欠測値の補間をしてもよい。さらに欠測値の補間後に、補間修正をしてもよい。さらに、処理部１は、原データとフィッティングされた関数による計算値の差を取り、差がしきい値よりも大きなものを取り除いても良い。また、処理部１は採用すべきデータのみを用いてフィッティング、フィルタリング（カルマンフィルタ）等の各処理を繰り返すことができる。
【００６１】
つぎに、処理部１は、定められた係数ｃｋに従い、関数によりノイズが除去された推定値を第１推定観測データファイル４４に記憶する、又は出力部３に出力若しくは表示する。また、処理部１は、第１パラメータファイル４６から係数を読み出し、カルマンフィルタ処理の初期値を式（１３）のように設定する（Ｓ１０９）。
【００６２】
処理部１は、ｙ（ｔ｜Ｔ）ｏｌｄ とＫの値を、ぞれぞれ初期値として０（ｙ（ｔ｜Ｔ）ｏｌｄ＝０，Ｋ＝０）とする（Ｓ１１１）。つぎに、処理部１は、式（９）に基づいて、カルマンフィルタ処理を行う（Ｓ１１３）。さらに、処理部１は、式（１２）を用いて、データの平滑化する（Ｓ１１５）。
【００６３】
処理部１は、Ｋの値を１繰り上げる（Ｋ＝Ｋ＋１） （Ｓ１１７）。さらに、処理部１は、時期Ｋに、推定値と実際の画素値のずれを求め、この値によって次に行う処理を判断する。即ち、カルマンフィルタ処理が収束してなく、このずれが予め設定したしきい値εよりも大きい場合には、次にステップＳ１２１の処理を行い、一方収束して、しきい値εよりも小さい場合には次にステップＳ１２５の処理を行う（Ｓ１１９）。
【００６４】
即ち、収束していない場合、処理部１は、計算パラメータｃ（１｜０）とＶ（１｜０）の初期値を更新し（ｃ（１｜０）＝ｃ（１｜Ｔ）、Ｖ（１｜０）＝Ｖ（１｜Ｔ））（Ｓ１２１）、さらに、ｙ（ｔ｜Ｔ）ｏｌｄの値をｙ（ｔ｜Ｔ）ｏｌｄ＝ｙ（ｔ｜Ｔ）とする（Ｓ１２３）。一方、収束した場合には、処理部１は、計算結果である推定された係数ｃｉ（ｔ）を第２のパラメータファイル４７に記憶し、ノイズを除去した推定された観測値を第２の推定観測データファイル４５に記憶する、又は、出力部３に出力若しくは表示する（Ｓ１２５）。
【００６５】
５． 実際の衛星データ処理例
ＬＭＦ処理を行うためには、地理補正、センサー補正済みであることが必要である。すべてのシーンで、それぞれの画素の地理座標が揃っている必要がある。そのようなデータは、農林水産研究計算センター衛星データアーカイブシステムであるＳＩＤａＢ(文献１７，１８参照)、Ｐａｔｈｆｉｎｄｅｒ プロジェクト(文献１９参照)、ＵＳＧＳ(文献２０参照)などのｗｅｂサイトからネットワークを通じて入手できる。また、ＳＰＯＴ／ｖｅｇｅｔａｔｉｏｎであれば、切り出し範囲と投影法が揃っていれば、ほとんど前処理なしでＬＭＦ処理することが可能である。
【００６６】
計算上、ＬＭＦ処理を行うには、最低１年分（１０日間合成で３６シーン）が必要である。特に雲の影響が強い熱帯では、１年半（５４シーン）以下では、結果が期待できない。３年以上のデータの蓄積があれば、一ヶ月程度のシーン欠損があってもほぼ問題なく処理することが可能である。計算はかなり膨大なものになるので、並列機などの使用が望ましい。（並列機を用いての処理については付記を参照のこと）
【００６７】
５．１． Ｐａｔｈｆｉｎｄｅｒ データ
図５に、ＬＭＦ処理前後のＮＤＶＩデータ図を示す。この図は、Ｐａｔｈｆｉｎｄｅｒデータ（８ｋｍ ＮＯＡＡ／ＡＶＨＲＲ １０日間合成）を用いてＮＤＶＩのＬＭＦ処理を行った例である。対象画素としては、北緯３６．３１度、東経１３６．９２度（岐阜県大野郡白川村）を選んだ。データ収集期間は９５年１月上旬〜９９年１２月下旬の１８０シーンを用いた。
【００６８】
Ｐａｔｈｆｉｎｄｅｒ データには、画素毎に、データ取得時の雲の有無についてのフラグ（ＣＬＡＶＲ データ）が添付されている。ＣＬＡＶＲの値が、１〜１１ は“ｃｌｏｕｄｙ”、１２〜２２は”ｍｉｘｅｄ”、２３〜３１は”ｃｌｅａｒ”にカテゴリが分かれている(文献２１参照)。雪の影響の無い、９５年５月〜９月のＮＤＶＩデータについて、日本付近を例に取り、ＬＭＦ処理の精度について検証した。
【００６９】
図６に、ＬＭＦ処理前後のデータの相関図を示す。この図は、処理データのＮＤＶＩを横軸にとり、ＬＭＦ処理済みデータを縦軸としたものである。処理前のＮＤＶＩについて０．０１毎に区間を区切り、その区間に属する同じ時期、同一画素のＬＭＦ処理後のＮＤＶＩの平均値と標準偏差をプロットした。雲の影響のない“Ｃｌｅａｒ”カテゴリに属する画素では、ＬＭＦ処理前後の値に高い相関が見られるのに対し、部分的に雲が混ざった“Ｍｉｘｅｄ”や雲に覆われた”Ｃｌｏｕｄｙ”カテゴリに属する画素ではＬＭＦ処理後の値が高くなっている。
【００７０】
雲に覆われた画素のＮＤＶＩの真値を知ることは不可能であるが、少なくとも、雲の影響の無い画素については、ＬＭＦ処理の前後でほぼ値が保存されている。
【００７１】
熱赤外バンドデータの処理：
図７に、ＡＶＨＲＲ／Ｃｈ４のデータをＬＭＦ処理したデータ図を示す。この図によると、ＮＤＶＩの時と同様、雲などの影響が軽減されていることがわかる。
最小値画像：
これまで、ＮＤＶＩ最大値画像は作られてきたが、ＮＤＶＩ最小値画像は雲の影響を受けた画素ばかりを集めてしまうために、意味づけをすることはできなかった。ＬＭＦ処理により雲の影響を除去することではじめて、意味のある最小値画像が得られたことになる。
ここで得られたＮＤＶＩ最小値画像は、年間を通じて植物が多く存在するか否かを示している。また、ＮＤＶＩだけでなく、熱赤外バンドの輝度温度に関しても、ＬＭＦ処理済みデータであれば、最小値の意味づけが可能である。
【００７２】
５．２． ＵＳＧＳ １ｋｍ 全球ＡＶＨＲＲデータ
米国地質調査所（ＵＳＧＳ）にアーカイブされている１ｋｍ解像度の ＡＶＨＲＲ のＮＤＶＩについて、タイ周辺のデータの処理前と処理後のものを比較したところ、タイは１０月〜３月までが乾期、４月〜９月までが雨期であるが、雨期の雲の影響を抑えて植生が繁茂する様子がよく捉えられている。また、データの欠損や原データにおけるパス間の質の違いも吸収できていることがわかる。また、１画素毎に、ＮＤＶＩの最小値と最大値、季節変化パターンが得られるので、植生タイプの分類だけでなく、植生の乾燥度や植生由来の地上可燃物量についての知見を得ることができる。
【００７３】
５．３． 実際のデータ処理例：Ｐａｔｈｆｉｎｄｅｒデータの処理
図８は、ＰａｔｈｆｉｎｄｅｒデータのＮＤＶＩをＬＭＦ−ＫＦ処理した結果を示す図である。ＬＭＦとＬＭＦ−ＫＦの結果を比較すると、ＬＭＦ処理の結果が９５年〜９９年の５年間の平均的な季節変化を表しているのに対し、ＬＭＦ−ＫＦでは、雲の影響を取り除きながら、年々変動にも追従している。注目すべきは、９６年早春は値の立ち上がりが遅く、９８年では逆に値の立ち上がりが早いのがはっきり現れていることである。過去の気象データによると、９６年は北陸地方の積雪量は例年と比較して多かったが、９８年早春はエル・ニーニョ現象によって暖冬となり、積雪量が異常に少なかったことが報告されている。
【００７４】
また、図９は、ＡＶＨＲＲ／Ｃｈ４輝度温度データをＬＭＦ−ＫＦ処理した結果を示す図である。ＮＤＶＩデータと同じ傾向を示しており、ＬＭＦ−ＫＦ では９６年、９８年の特徴がはっきりと現れている。
【００７５】
本実施の形態では、ボルネオ周辺のＮＤＶＩデータについて、９５年〜９９年の５年間のＬＭＦ−ＫＦ処理を行い、１年毎にカラー合成をしたものを求めたところ、９５年〜９９年になるにしたがって、年間を通じてＮＤＶＩの高い場所が減少し、代わって、植生指数の低い赤色の部分が増加している。また、季節変動の大きいことを示す白っぽい部分も増えている。９８年以降、ボルネオ島東部に植生の少ない地域が増えている。これは９８年にエル・ニーニョ現象のため大規模森林火災が発生したが、その延焼地域を表している。
【００７６】
ＬＭＦ処理を用いれば年々変動の解析も可能であるが、データを期間毎に切り出して個別に処理する必要があった。さらには、処理の性質としてデータの収集期間が短くなると結果に雲の影響が含まれやすくなる（モデルが不安定化する）ことや、データ収集期間のつなぎ目での計算値の整合性がとれない場合があった。
【００７７】
しかし、ＬＭＦ−ＫＦでは、長期間のデータを一度に処理しても、それぞれの年の特徴が抽出できることが示された。データの収集期間が長いので、モデルの不安定性が軽減されるためであると考えられる。
【００７８】
６． 時系列データ処理法のまとめ
今回開発した、時系列フィルタと状態空間モデルを用いた、高頻度観測衛星データ処理法では、データ処理期間全体の平均的な季節変化を捉えるだけでなく、雲やヘイズ、システムノイズなどの影響を抑えながら、年々変動の抽出も可能であることがわかった。
【００７９】
この方法により、植生指数だけでなく、熱赤外バンドの輝度温度データについても衛星データから季節変動プロファイルを抽出できるようになったが、これにより、一週間〜１０日毎の時間間隔で全球規模の植生環境の変動を監視することがはじめて可能となった。特に、大規模森林火災や乾燥害などの農業・林業災害は植生の乾燥の程度に影響されるため、植生の状態を迅速かつ的確に把握する手法の開発は防災面で大きな意味を持っている。
【００８０】
植生指数や輝度温度の季節変動データから、森林火災危険度を算出できるが、それについては次節で述べる。
また、植生指数の季節変動パターンの分類によって、広域での土地被覆分類も可能である。
【００８１】
ＬＭＦおよび ＬＭＦ−ＫＦは、画素毎に処理が独立しているため、並列機に適している。したがって、ＭＯＤＩＳなど、従来より解像度の高い高頻度観測衛星データについても、時間をかけずに処理が行える利点がある。
【００８２】
７． 輝度温度データを用いた火災検出
つぎに、ＬＭＦおよびＬＭＦ−ＫＦによって処理されたデータの活用、特に輝度温度データを用いた火災検出の高精度化について言及する。
ＤＭＳＰ／ＯＬＳでの雲域の自動推定：
【００８３】
ＤＭＳＰ／ＯＬＳ での雲域除去は、特定のしきい値を定めるか、ＮＯＡＡ／ＮＧＤＣから配信される全球地表面温度データと比較することで行っている。前者の方法では、しきい値の設定に恣意性があったり、雲除去が適切に行われないなど、精度において問題がある。一方、配信される全球地表面温度データは、配信時刻が衛星データと同時でなく１日程度遅れており、リアルタイム処理には間に合わない。また、地上観測での地表面温度と衛星で観測される熱赤外バンドの
輝度温度は必ずしも同じ物理量ではない。
【００８４】
ＤＭＳＰ／ＯＬＳの熱赤外輝度温度をＬＭＦ処理しておけば、各時期、各画素での雲の影響のない値が得られるので、雲の判定は単純にしきい値を設定すれば足りる。また、観測波長の問題はあるが、次節に示す方法により、熱によっても火災検出が可能になる可能性がある。
ＮＯＡＡ／ＡＶＨＲＲ での火災検出アルゴリズムの単純化：
【００８５】
現在、ＮＯＡＡ／ＡＶＨＲＲ でのホットスポット検出は、単純に輝度温度にしきい値を設定するものや、バンド間演算をしたもので判断するものが多い(文献２４参照)。この方法では、どうしても、裸地や火口、人口構造物（石油および製鉄プラント）などの火災でないものも検出してしまう欠点がある。これらを回避するため、複雑な火災検出アルゴリズムも開発されているが２４）、演算量の点から、リアルタイム処理には適していない。さらには、土地被覆によって火災検出に用いるしきい値だけでなく、アルゴリズムも変更する必要があるとの指摘もなされている。これを実現することを考えた場合、特定の狭い範囲で火災検出システムを構築することは可能であっても、国レベルの広域を対象にすることは困難である。
【００８６】
ここで、ＬＭＦ処理済みの、雲無しの輝度温度データがあれば、それらは、対象となる画素の土地被覆や緯度、標高などの影響がすでに含まれているので、単純にしきい値を設定するだけでも火災精度の向上が期待できる。また、次節以下で述べる方法により、画素内の火災割合を算出することも可能となる。
衛星データによる小規模火災の検出：
【００８７】
Ｐｒｉｎｃｅ らが、静止気象衛星 ＧＯＥＳ を用いてブラジルの森林火災監視を行っている(文献２５参照)。一般に、静止気象衛星は高度が非常に高いため、解像度は低くなってしまうが、画素内の火災割合を算出することで、ＮＯＡＡ／ＡＶＨＲＲ と同程度の火災検知を可能にしている。ここで、熱赤外バンドｉでの、温度Ｔの黒体の分光放射輝度をＬｉ（Ｔ）、画素内の火災割合をｐとすると、
Li(Tobs)＝ｐ・Li(Tfire)＋(1−ｐ)・Lｉ(Tｂ) （１５）
の式が成り立つ。ただし、Ｔｏｂｓは観測されたものの温度、Ｔｆｉｒｅは火災温度、Ｔｂは画素内の火災のない部分の温度である。この式には、未知数がｐ、Ｔｆｉｒｅ、Ｔｂの３個あるので、熱赤外バンドは３つ必要である。しかし、Ｐｒｉｎｃｅらは、Ｔｆｉｒｅ＝４００（Ｋ）、Ｔｂを周囲の比較的低温な画素の平均と仮定して火災割合を算出している。しかし、これでは、火災のない場合、対象となる画素と、その周囲の画素の温度が同じであることを仮定している。換言すれば、対象となる画素と、その周囲の画素の土地被覆や標高が同一であることになるが、この仮定が成り立つ場所はかなり限られてくる。ここで、Ｔｂとして、対象となる画素の同時期のＬＭＦ処理された値を使えば、ＤＭＳＰ／ＯＬＳのように熱赤外バンドが一つしかなくても、火災割合ｐを算出することが可能になる。また、ＮＯＡＡ／ＡＶＨＲＲのように、複数の熱赤外バンドがあれば、火災割合と火災温度の両方を算出することが可能である。
【００８８】
ＮＯＡＡ／ＡＶＨＲＲでは、センサーが飽和してしまう温度が低いため、火災温度の推定は困難であるが、ＭＯＤＩＳ等の新型センサーでは、火災温度による消火優先度の設定なども考えられる。
この方法は、ＧＭＳ−５の後継機（ＭＴＳＡＴ）のデータにも応用が可能であり、毎時の火災情報提供への応用が期待できる。
【００８９】
８．付記及び参考文献
本発明の地球観測衛星データのノイズ除去方法又は地球観測衛星データのノイズ除去装置・システムは、その各手順をコンピュータに実行させるための地球観測衛星データのノイズ除去処理プログラム、地球観測衛星データのノイズ除去処理プログラムを記録したコンピュータ読み取り可能な記録媒体、地球観測衛星データのノイズ除去処理プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
以下に、参考文献を示す。
【００９０】
1) インドネシア森林火災予防プロジェクト、“インドネシア国立公園森林火災調査報告書”、国際協力事業団、1998
2) 沢田治雄、“大規模森林火災へのリモートセンシングの利用”、森林科学、24、22-28、1998
3) 本多嘉明、村井俊治、加藤喜久雄、“世界植生モニタリング」、日本写真測量学会誌「写真測量とリモートセンシング」、31、1、4-14、1992
4) Defries,R. S. and J. R. G. Townshed, “NDVI-derived land cover classifications at a global scale”, Int. J. Remote Sens., 15, 3567-3586, 1994.
5) Myneni, R. B., C. J. Tucker, G. Asrar and C. D. Keeling,“Interannual variations in satellite-sensed vegetation index data from 1981 to 1991”, J. Geophys. Res., 103, D6, 6145-6160, 1998.
6) Los, S. O., C. O. Justice and C. J. Tucker, “A global 1 deg. by 1 deg. NDVI data set for climate studies derived from the GIMMS continental NDVI Data.”, Int. J. Remote Sens., 15, 3493-3518, 1994.
7) GLCC web site : http://edcdaac.usgs.gov/glcc/globdoc1_2.html
8) Viovy, N., O. Arino and A. S. Belward, “The Best Index Slope Extraction (BISE): A method for reducing noise in NDVI time-series”, , Int. J. Remote Sens., 13, 1585-1590, 1992.
9) 李 雲慶、大沼一彦、安田嘉純、“時系列植生指数データのフーリエ解析”、日本写真測量学会誌「写真測量とリモートセンシング」、31、4、4-14、1992
10) Sellers, P. J., C. J. Tucker, G. J. Collatz, S. O. Los, C. O. Justice, D. A. Dazlich and D. A. Randall, “A global 1 deg. by 1 deg. NDVI data set for climate studies. Part 2: The generation of global fields of terrestrial biophysical parameters from the NDVI”, Int. J. Remote Sens., 15, 3519-3545, 1994.
11) Verhoef, W., M. Menenti and S. Azzali, “A colour composite of NOAA-AVHRR-NDVI based on time series analysis (1981-1992)”, Int. J. Remote Sens., 17231-225, 1996.
12) Ricotta, C., G. Avena and A. De Palma, “Mapping and monitoring net primary productivity with AVHRR NDVI time-series: statistical equivalence of cumulative vegetation indices”, ISPRS J. Photogramm. Remote Sens., 54, 325-331, 1999.
13) Azzali, S. and M. Menenti, “Mapping vegetation-soil-climate complexes in southern Africa using temporal Fourier analysis of NOAA-AVHRR NDVI data”, , Int. J. Remote Sens., 21, 973-996, 2000.
14) Roerink, G. J., M. Menenti and W. Verhoef, “Reconstructing cloudfree NDVI composites using Fourier analysis of time series”, Int. J. Remote Sens., 21, 1911-1917, 2000.
15) 澤田義人、沢田治雄、斎藤英樹、“高頻度観測衛星による植生の季節変動観測”、日本写真測量学会講演論文集、2000秋季、209-212、2000
16) 北川源四郎、“FORTRAN77 時系列解析プログラミング”、岩波書店、第５章、1993
17) http://www.affrc.go.jp/Agropedia/
18) 児玉正文、宋 献方、“農林水産リモートセンシングデータベースシステムの構築”、日本リモートセンシング学会第２８回学術講演会論文集、259-260、2000
19) Pathfinder web site:
http://daac.gsfc.nasa.gov/CAMPAIGN_DOCS/LAND_BIO/GLBDST_main.html
20) http://edcdaac.usgs.gov/1KM/1kmhomepage.html
21) Agbu, P. A. and M. E. James, “NOAA/NASA Pathfinder AVHRR Land Data Set User's Manual. Goddard Distributed Active Archive Center”, NASA Goddard Space Flight Center, Greenbelt, 1994.
http://daac.gsfc.nasa.gov/CAMPAIGN_DOCS/LAND_BIO/AVHRR_GD.pdf
22) 北川源四郎、“FORTRAN77 時系列解析プログラミング”、岩波書店、第9章、1993
23) 西山清、“周波数が既知である成分のみを含む時変スペクトルの一推定法”、電子情報通信学会論文誌、J74-A、6、916-918、1991
24) Martin, M. P., P. Ceccato, S. Flasse and I. Downey, “Fire detection and fire growth monitoring using satellite data”, Remote Sensing of Large Wildfires in the European Mediterranean Basin, E. Chuvieco (ed.), Splinger-Verlag, chapt. 6, 1999.
25) Prins, E. M. and W. P. Menzel, “Geostationary satellite detection of biomass burning is South America”, Int. J. Remote Sens., 13, 2783-2799, 1992.
【００９１】
【発明の効果】
本発明によると、高頻度観測衛星からのデータの画素毎の時系列的な特徴を保存し、時間分解能を完全に保持したまま、雲などの影響およびシステムノイズの除去をすることができる。
また、本発明によると、時間分解能を保持したまま植生などの変動の様子が捉えられるため、防災、環境、農業分野での広域現状把握、および予測が容易となる。さらに、本発明によると、適切な関数を自動的に選択し、ノイズ除去処理を自動で行えるものであり、定常業務の一部分とすることが可能である。さらに本発明によると、水分指数など他の地球高頻度観測衛星のデータにもそのまま転用および応用ができる。
【図面の簡単な説明】
【図１】シミュレーションによるフィルタ特性の評価図
【図２】シミュレーション結果の例を示す図
【図３】命令列スケジュール処理を実行するための装置の構成図
【図４】データ処理のフローチャート
【図５】ＬＭＦ処理前後のＮＤＶＩデータ図
【図６】ＬＭＦ処理前後のデータの相関図
【図７】ＡＶＨＲＲ／Ｃｈ４のデータをＬＭＦ処理したデータ図
【図８】ＰａｔｈｆｉｎｄｅｒデータのＮＤＶＩをＬＭＦ−ＫＦ処理した結果を示す図
【図９】ＡＶＨＲＲ／Ｃｈ４輝度温度データをＬＭＦ−ＫＦ処理した結果を示す図
【符号の説明】
１ 処理部
２ 入力部
３ 出力部
４ 記憶部
４１ 計算パラメータファイル
４２ 第１の観測データファイル
４３ 第２の観測データファイル
４４ 第１の推定観測データファイル
４５ 第２の推定観測データファイル
４６ 第１のパラメータファイル
４７ 第２のパラメータファイル
【発明の属する技術分野】
本発明は、衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体に係り、特に、高頻度観測衛星データのノイズ除去処理による植生の季節変動観測法に関する。
【従来の技術】
【０００２】
１９９７年から１９９８年にかけてインドネシアで大規模な森林火災が発生し、煙害等で周辺諸国も含めて甚大な影響を及ぼした（文献１，２参照。なお、文献は後述する。以下同様。）。その後も東南アジア地域では森林火災が頻発し、大規模火災の再発生が懸念されている。現在、国際協力事業団（ＪＩＣＡ）の「インドネシア森林火災予防計画プロジェクト」の第二フェーズが進行中であり、その中で、「森林火災早期発見システム」が稼働している。これは、火災発見を衛星データを用いて行い、森林管理を担当する機関に通報するものである(文献２参照）。一方、火災予防の観点からは火災危険度の評価が待たれている。
【０００３】
火災危険度評価の基礎となるのは森林植生図であるが、乾燥に対する森林の感受性を評価する必要があることを考えると、単なる森林区分では不十分で、時間を含んだ情報に基づいて危険度を算出することが期待されている。
【発明が解決しようとする課題】
【０００４】
従来の植生の季節変動観測法のための時系列処理手法の開発に対して、その課題について以下に説明する。
【０００５】
従来のＮＯＡＡ／ＡＶＨＲＲ あるいはＳＰＯＴ／ｖｅｇｅｔａｔｉｏｎ、ＭＯＤＩＳなどの高頻度観測衛星からのデータによって広域での植生観測が行われている(文献３−１３参照)。特に、大規模森林火災や乾燥害などの農業・林業災害は、植生の乾燥の程度に影響されるため、植生の状態を迅速かつ的確に把握する手法の開発は防災面で大きな意味を持っている。しかし、光学センサーでは、雲やヘイズの影響が避けられず、的確な植生分布およびフェノロジーの把握が困難である。特に、熱帯を対象とする場合、数ヶ月間にもわたって雲の影響があり、雨期にはほとんど有効なデータが観測できないこともある。
【０００６】
また、従来より広域での植生分布についても数多くの研究がある(文献３−７参照)。例えば、ＧＬＣＣ（Ｇｌｏｂａｌ Ｌａｎｄ Ｃｏｖｅｒ Ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ ｄａｔａ ｂａｓｅ）があるが、日本にサバンナが出現するなど誤分類が見受けられる(文献７参照)。これも雲の影響を十分に取り除くことができなかったのが原因である。
【０００７】
さらに、これまでにも、高頻度観測衛星のデータから雲の影響を取り除いてＮＤＶＩ（植生指数）の季節変化を抽出する試みがなされている(文献８参照)。例えば、解像度を落として、複数の画素の最大値を取るものや、１０日〜１ヶ月のデータの最大値を取るものなどがある。しかし、これらは、解像度が低かったり、データの取得間隔が一ヶ月であったりと即時性の要求される防災データとしては使いにくい。データの合成期間が長くなれば、時間的な分解能が低下する。また、これらは、植生の変動を一週間から１０日の間隔で把握する必要のある農林業分野では応用できないという課題がある。
【０００８】
また、Ｆｏｕｒｉｅｒ展開類似の方法を用いて、植生指数の季節変化を三角関数で表すこともなども行われている(文献９−１４参照)。この方法にはどの関数を使うかについて恣意性があり、様々な土地被覆を含む広域データの処理には適していない。また、この方法では、関数の当てはめにおいて、処理に用いる三角関数群はあらかじめ決められており、必ずしも処理に適した関数群が選ばれていないという課題がある。
【０００９】
また、１６×１６などのウインドウを設定し、その中での最大値をウインドウ内のすべての画素値に置き換える方法もある。しかし、このような方法では、空間分解能が低下するという課題がある。
【００１０】
本発明は、以上の点に鑑み、空間分解能を保持したまま、任意の時間分解能を実現させながら、雲やヘイズの影響を取り除くために、時系列モデルを用いた衛星データ処理手法ＬＭＦ（Ｌｏｃａｌ Ｍａｘｉｍｕｍ Ｆｉｔｔｉｎｇ：局所最大値フィッティング）を提案することを目的とする。すなわち、本発明は、画素毎の時系列的な特徴を保存し、原データの時間分解能を完全に保持したまま、雲などの影響およびシステムノイズの除去を目的とするものである。
【００１１】
また、本発明は、時間分解能を保持したまま植生などの変動の様子が捉えられるようにすることで、防災、環境、農業分野での広域現状把握、および予測が容易とすることを目的とする。さらに、本発明は、適切な関数を自動的に選択し、ノイズ除去処理を自動で行えるようにすることにより、定常業務の一部分とすることを可能とすることを目的とする。さらに本発明は、水分指数など他の地球高頻度観測衛星のデータにもそのまま転用および応用を可能とすることを目的とする。
【００１２】
【課題を解決するための手段】
本発明の特徴のひとつとしては、例えば、処理対象画素の時系列データについて、前後ｎ時期について、局所最大値を抽出するフィルターを作用させる。そこで得られた局所最大値群に関数当てはめを行い、画素値を推定する。ある時期に、推定値と実際の画素値のずれが、設定したしきい値よりも大きい場合には雲やシステムノイズの影響とみなしてその時期のデータは採用しない。採用されたデータだけを用いて再度、局所最大値の算出と関数当てはめを行った値で画素値を置き換える。必要であればこの処理を繰り返し、結果が収束するまで行う。
【００１３】
本発明の他の特徴としては、例えば、画素毎の時系列的な局所最大値を求め、それらに対して関数の当てはめをくりかえし行うことで、有効なデータだけを用いた推定を行うものである。上記の処理を画像全体に対して行うため、雲やシステムノイズの影響が除去されるという作用を有する。各画素毎の処理は独立しているため、原データの空間分解能は完全に保持される。
【００１４】
また、本発明は、時期ｋ近傍の局所的最大値を、時間遅れを少なくするように抽出する作用を有する。これにより、雲などの影響を受けた時期の画素値が除去されるという作用を有する。
【００１５】
本発明の第１の解決手段によると、
処理部は、観測対象エリア内の各画素に対応して、衛星による植生の季節変動についての時系列の第１の観測データを記憶した第１の観測データファイルから、画素毎の第１の観測データを読み込むステップと、
【００１６】
処理部は、各画素毎に、ある時期又は時点に対応して、該時期又は時点の前後所定期間の局所最大値を抽出することにより時系列の第２の観測データを求め、それを各画素に対応して第２の観測データファイルに記憶するステップと、
【００１７】
処理部は、第２の観測データファイルに記憶された第２の観測データを、長期的非周期項及び季節変化項を含む次式の関数を用い、各第１の係数ｃｉを最小二乗法又は他の推定法により決定することでフィッティングし、各第１の係数ｃｉを第１のパラメータファイルに記憶するステップと、
【００１８】
【数１０】
ｃｉ：求める第１の係数
ｔ：時期又は時間
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列（例：１，２，３，４，６，１２）
Ｍ：１年間のデータ数（シーン数） ］
【００１９】
処理部は、第１の係数ｃｉで定められた前記関数に従い、第２の観測データから雑音データを除去した第１の推定観測データを求め、それを第１の推定観測データファイルに記憶する、又は出力部に出力、若しくは表示するステップと
を含む地球観測衛星データのノイズ除去方法、これらのステップをコンピュータに実行させるためのノイズ除去処理プログラム、及びノイズ除去処理プログラムを記録した記録媒体が提供される。
【００２０】
【発明の実施の形態】
本実施の形態では、植生の季節変動観測のために、時系列モデルを用いた新しい衛星データ処理法を提案する。また、雲の影響を取り除いた植生の季節変動に関する観測データから、火災危険度を算出した。以下、火災危険度評価の基礎となる植生の季節変動観測のための衛星データ処理法について述べる。なお、植生の季節変動観測データとしては、例えば、植生指数データ、水分指数、中間赤外線バンド、熱赤外線輝度バンド、温度データ等である。また、観測データに含まれる可能性がある雑音データは、例えば、雲、霧、雨、ヘイズ等の大気状態に起因するもの、センサ感度等のシステムに起因するもの等である。
【００２１】
１． ＬＭＦの理論
ＬＭＦ処理は、時系列フィルタリングと関数によるフィッティング手法を組み合わせたアルゴリズムで、対象画素の時系列データから、雲やヘイズなどの影響を除去し、季節変化パターンを抽出するものである。フィッティングに用いる関数を、（１）式に示す。
【００２２】
【数１１】
ｃｉ：求めるパラメータ：ｃ０は平均、ｃ１はトレンド係数、それ以外は周期関数のパラメータとなる。
ｔ：利用するデータの時間間隔に基づく時間単位数（１０日合成のＮＯＡＡデータを用いる場合なら、１０日が１単位）
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：このソフトウェアの構造は周期関数の種類には依存しないため、任意の周期関数と置き換えて利用することが可能。
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列：１，２，３，４，６，１２。これによって、１ヶ月周期、２ヶ月周期、３ヶ月周期、４ヶ月周期、６ヶ月周期および１２ヶ月周期の周期関数（ｓｉｎ、ｃｏｓ）が当てはめられる。
Ｍ：１年間のデータ数：１０日間最大値合成のＮＯＡＡデータを用いる場合、１年間３６シーンになるため、Ｍ＝３６である。
なお、この場合、シーン数Ｎと１年間のデータ数Ｍが一致しないので、一般的な離散 Ｆｏｕｒｉｅｒ級数は計算できない。
【００２３】
１．１． シミュレーションによる処理特性の評価(文献１５参照)
植生指数や熱赤外線輝度温度データでは、観測データには雲などの影響がすでに含まれており、真値は不明である。そこで、本実施の形態では、シミュレーション的手法を用いてフィルタ特性の評価を行った。
まず、図１に、シミュレーションによるフィルタ特性の評価図を示す。例えば、（１）式の係数の組ｃに、発生した乱数（Ｓ１１）を代入して、時系列データを１００００組発生させる（Ｓ１２）。解析のパラメータは、一例として、Ｎ＝５４、Ｍ＝３６とした。次に、雲やヘイズの影響を受けた状況をシミュレートするため、ノイズを発生し（Ｓ１３）、そのノイズを元の値の１０％〜７０％のものに置き換えて、ノイズが加わった状況をシミュレートした（Ｓ１４）。ノイズが加算される確率は各時期について０．５とした。そして、ノイズ加算後の時系列データに対して、ノイズ除去処理し（Ｓ１５）、その結果を調べた（Ｓ１６）。結果の解析は発生させた係数の組ｃと、処理後に再生された係数の組ｃ’とについて、単回帰分析を行った。
【００２４】
１．２． 最小二乗法によるフィッティング
一般に、ノイズ加算後のデータをそのまま（１）式でフィッティングした場合、どの項も係数の約７割が保存される。これに、（２）式で表されるような局所最大値を抽出するような時系列フィルタを作用させる。なお、この局所最大値の計算式は一例であり、適宜の式を用いることもできる。
【００２５】
【数１２】
ｄｔ’ は時間単位ｔの時の局所最大値 （第２の観測データ）
ｄｔ は時間単位ｔの時のデータ値 （第１の観測データ）
ｗ は局所最大値を求める片側時間単位数（ウインドウサイズ）
である。
【００２６】
（２）式では、ｄｔにおける左右ｗ間のデータの最大値をもとめ、その小さい方を局所最大値としている。この考え方が本発明のひとつの特徴である。本実施の形態では高頻度観測衛星からの観測データｄｔ（第１の観測データ）に対して、（２）式のようなフィルタリング処理を実行して局所最大値ｄｔ’（第２の観測データ）を求め、さらに、（１）式の関数で最小二乗法によりフィッティングして各係数ｃｉを求める。これによって、データに影響を与える雲やシステムノイズを軽減して、季節変化する地表の情報として時間単位ｔに近い時期のデータを自動抽出することができる。
【００２７】
１．３． フィッティング関数の自動選択
上述のような式（２）によるフィルタと最小二乗法によるフィッティングを併用したものでは、定数項、一年周期、半年周期項はほぼ回復できるが、それら以外の短周期関数項の係数については保存されない場合がある。これを改善するため、本発明は、さらに以下のように、（１）式の関数のなかで、必要なもの以外はフィッティングに用いないようにしてもよい。なお、最小二乗近似に対しても、適宜のフィッティング方法を用いることができる。
【００２８】
いま、最小二乗近似を用いて時系列データを関数でフィッティングすることを考える。フィッティングに用いる関数の個数（種類）が多ければ残差は減少するが、フィッティングした結果は不安定（結果がデータに大きく依存する）となる。赤池らによると、赤池の情報量基準（ＡＩＣ）が最小となるようなモデルが最適であるとされている(文献１６参照)。残差が正規分布に従うと仮定されているから、残差の分散をσとするとＡＩＣは、対数尤度とパラメータ数の差で表される（３）式のようになる。
【００２９】
【数１３】
N は所定期間のデータ数
j は用いる関数の個数（周期関数のペア数）
σはオリジナルデータとフィッティングによるデータの残差の分散
本実施の形態ではAICの値が最小になるように、周期関数の組み合わせを自動的に選択するようにした。即ち、この値が最小になるようにjを選ぶようにした。なお、正の定数倍は結果を変えないため問題ない。この例においては、シミュレーションの結果、三ヶ月周期より短い周期の項については改善がみられたが、逆に四ヶ月周期の項についての相関係数は低くなる。
【００３０】
図２は、シミュレーション結果の例を示す図である。
１．４． フィッティング処理の繰り返しによる高精度化
本実施の形態ではさらに、処理の高精度化を実現するため、フィッティングを繰り返すアルゴリズムを採用した。最初に、原データに局所最大値抽出フィルタを作用させた後、フィッティングする。その結果（計算値）と、原データとの差をとり、しきい値よりも大きなものを取り除く。データレンジの１５〜２０％をしきい値とすれば、１年および半年に相当する周期を持つ関数の係数は９５％以上、それ以外の短周期成分も８０％が回復できることが明らかになった(文献１５参照)。
【００３１】
１．５． 熱赤外バンドデータの処理
ＮＯＡＡ／ＡＶＨＲＲの熱赤外バンドでは、雲のあるピクセルは温度が低くなる性質を用いれば、ＮＤＶＩと同じ計算方法を用いればよいと考えられる。したがって、本発明は輝度温度についても処理が可能である。
【００３２】
２． カルマンフィルタを用いたＬＭＦアルゴリズムの高精度化
２．１． ＬＭＦ処理の問題点
ＬＭＦは高頻度観測衛星からの植生指数や輝度温度データから雲の影響を軽減し、季節変化パターンを抽出するのに有効な手法である。しかし、処理期間すべてにわたってフィッティング処理を行うため、どうしても平均化された結果が得られてしまう。また、データ収集期間が短い（１年程度）場合には、フィッティングの結果が不安定になりやすく、年々変動の抽出は困難な場合があった。そこで、本実施の形態ではさらに、データ処理に状態空間モデルを採用し、データ処理を行う手法を開発した（文献２２参照）。
【００３３】
また、データを防災に役立てようとするならば、単に現況を把握するだけでなく、短期的にせよ状況の予測をしながら対策を講じなければならない。したがって、時系列データの解析についても、予測が可能な手法を採用する必要がある。状態空間モデルであれば、この点も問題なく適用することができる。本実施の形態では、一例として、状態空間モデルのパラメータの推定には、Ｋａｌｍａｎの提唱した手法であるかルマンフィルタ（Ｋａｌｍａｎ Ｆｉｌｔｅｒｉｎｇ）を応用するが、それに限らず、適宜の推定法を用いることができる。
【００３４】
カルマンフィルタは、人工衛星の軌道・姿勢の予測と制御のように物理的に定式化が可能なものだけでなく、降水量からの洪水予測や経済状況の時間的変化の予想のように、直接に原因と結果の関係が定式化できなくても統計量さえあればモデル構築が可能な点が大きな特徴である(文献２２参照)。さらに、
○ 少数の観測値しか得られなくても複雑な数値モデルを構築できる。
○ 最新の観測値によってパラメータ群を更新・最適化するので、周期的でない要因も考慮できる。
などの優れた特徴を有しているので、エル・ニーニョのように非定常的な擾乱要素の影響も含めたモデルの構築が可能である(文献２２参照)。
【００３５】
２．２． 時系列モデルの設定
高頻度観測衛星データを用いて得られるＮＤＶＩの季節変化プロファイルは、三角関数の組み合わせで表現できることがＬＭＦ処理の結果から明らかとなった。本実施の形態では、この知見を生かし、離散時変スペクトルの推定技術を応用した(文献２３参照)。すなわち、（１）式の係数ｃｋの組を状態空間モデルの状態ベクトルと考え、その値の推定（状態推定）をカルマンフィルタによって行う。
【００３６】
ＬＭＦでは、雲やヘイズの影響のない季節変化プロファイルを、（１）式の形で記述できるとした。ここで、右辺の第二項が非周期的変動を考慮した項である。この式での各係数ｃｋは、（１）式では固定であるが、これを時変として、季節変化プロファイルの変動に対応できるようにした。右辺第一項が時変になったことで、トレンドも表現できることから、第二項は削除できる。したがって、（４）式となる。
【００３７】
【数１４】
この式は先のLMFに対して、パラメータｃが全て時変変数としたものである。これによって、衛星が代わり、観測センサの感度が変わるなどのシステム変動（システム雑音の一部）があっても、季節変化を追うことができる。この点が、特に、このソフトウェアの画期的なところである。
（４）式のｃｋ（ｔ）の時間的変化は自己回帰モデルで現される。
【００３８】
【数１５】
ここでは、Ｋ＝０のときのモデル、
【００３９】
【数１６】
【００４０】
【数１７】
【００４１】
【数１８】
酔歩散乱モデルだけだとデータの追従性が低下し、年々変動を表すことができない場合がある。また、（７−２）式のモデルだけを採用すると、フィルタリング結果の値が大きく変動しすぎる場合がある。（６−１）式を採用することによって、センサ感度の違いなどのノイズによるデータの急激な変化を押さえられる。
K＝０以外の時のモデルは（６−２）式である。
【００４２】
【数１９】
【００４３】
【数２０】
ｃは求めるパラメータ
Ｆは予め定められた行列
Λはシステム雑音が、一度しか作用しないようにするための定数行列。
Ｈは周期関数を表す行列
ｖ（ｔ）は平均０、分散がσvの白色雑音（雲、ヘイズ等の大気の状態な等に起因する観測雑音）
【００４４】
【数２１】
また、各状態変数の雑音の分散Ｑ
【００４５】
【数２２】
三角関数による時変係数推定問題では、複素形式での表記が一般的であり、計算量も半減できるが(文献２３参照)、この例では、使用する関数を限定しないアルゴリズムを開発するため、敢えて複素形式にはしなかった。
【００４６】
実際にはＮＤＶＩデータや熱赤外領域の輝度温度データでは、雲の影響があれば、値が必ず低くなる性質がある。これは、観測雑音ｖ（ｔ）が有色となる可能性を含んでいる。したがって、この状態空間モデルで求められる状態変数の値は、必ずしも最適な量とは言えずに、準最適なものである可能性があるが、データ処理には時系列局所最大値フィルタを併用しているので、観測雑音の有色性は軽減されているものと考える。
【００４７】
２．３． カルマンフィルタ
次に、時変係数ｃ（ｔ）を、カルマンフィルタを用いて推定することについて説明する。ここで、ｃ（ｔ｜ｔ）を時刻ｔまでの観測値がわかっている場合の値とすると、一期先の予測はｃ（ｔ＋１｜ｔ）となる。
カルマンフィルタは（９）式の通りである。観測値をｙ（ｔ）とする。
【００４８】
【数２３】
これらの式を用いれば、Ｈ、Ｆ、Ｑのすべての要素があらかじめ定められているとすれば、ｃ（ｔ｜ｔ−１）やＶ（ｔ｜ｔ−１）が既知のとき、観測値ｙ（ｔ）から他の量を計算することができる。また、ｙ（ｔ）が何らかの理由で欠損になった場合には、カルマンゲインが計算できないので、（１０）式を用いて状態更新を行う。
【００４９】
【数２４】
【００５０】
【数２５】
ε０は正の数（０〜１）、Iは単位行列である。
を用いて最初の計算を行うことがしばしば行われる。
ところで、カルマンフィルタのアルゴリズムを用いて、時間をｔ＝Ｔからｔ＝１に時間を遡って処理すればデータの平滑化（ｓｍｏｏｔｈｉｎｇ）も可能である。それらの式は、
【００５１】
【数２６】
ここで述べた手法をＬＭＦ−ＫＦ（Ｌｏｃａｌ Ｍａｘｉｍｕｍ Ｆｉｔｔｉｎｇ ｗｉｔｈ Ｋａｌｍａｎ Ｆｉｌｔｅｒ ａｌｇｏｒｉｔｈｍ）と以後略記する。
この方法で、ベクトルＱの要素をすべて１０-4程度にすれば、事実上、各状態変数ｃｋは変化しなくなるので、従来のＬＭＦの結果を再現することも可能である。計算量は増大するが、ＬＭＦ−ＫＦはＬＭＦを完全に包含した手法であることがわかる。
また、このＬＭＦ−ＫＦ法の特徴として、各係数についての自己回帰モデルである（５）式の次数ｍを上げれば、長期間データ欠損があってもそれらを推定することが可能であり、また、そのままのアルゴリズムでデータ収集期間の終わりからの数年間先の予測までが可能である。
【００５２】
２．４． 実際のデータ処理手法
画素毎に行われる、実際のデータ処理の流れを以下に説明する。先にも触れたように、カルマンフィルタルゴリズムの適用にあたって、未知量であるシステム雑音の分散σＫおよび観測雑音の分散σｖはパラメータとして与えることにするが、ｃ（１｜０）、Ｖ（１｜０）も必要である。これらの要素を画素毎に推定するのは大変なので、最初に（１）式を適用してＬＭＦにより係数を求めて初期値とする。
【００５３】
【数２７】
ε0は正の数（０〜１）、Iは単位行列である。
本実施の形態では、これらの初期値を用いてフィルタリングと平滑化を一回行って、次回の計算の初期値とする。
【００５４】
【数２８】
このアルゴリズム（ＬＭＦ−ＫＦ）は、ＬＭＦと同様に、画素毎に処理が行われるため、ＬＭＦの時系列解析のサブルーチンを入れ替えるだけでよい。もちろん、プログラム内でのデータ構造も同じなので、並列化が容易である。
【００５５】
３．ハードウェア
図３に、命令列スケジュール処理を実行するための装置の構成図を示す。この装置は、処理部（ＣＰＵ）１、入力部２、出力部３、記憶部４を備える。処理部１は、ＣＰＵであり、命令列スケジュールプログラムによりスケジュールを実行する。記憶部４は、計算パラメータファイル４１、第１の観測データファイル４２、第２の観測データファイルファイル４３、第１の推定観測データァイル４４、第２の推定観測データファイル４５、第１のパラメータファイル４６、第２のパラメータファイル４７を有する。
【００５６】
４．フローチャート
図４に、実際のデータ処理のフローチャートを示す。なお、破線枠内の処理は、画素毎に行われる。
処理部１は、入力部２又は記憶部４の計算パラメータファイル４１から計算パラメータを読み込む（Ｓ１０１）。ここでは例えば、計算パラメータとして、画像数、像ライン数、画像コラム数、関数セット種類数、欠損シーン数、カルマンフィルタ使用時の関数セット数、局所最大値算出時の片側ウインドウサイズ雲検出閾値、カルマンフィルタ処理時の最大繰り返し回数、カルマンフィルタ収束判定値、システムノイズの分散パラメータなどがある。
【００５７】
処理部１は、衛星による画素データ毎の植生の季節変動についての時系列な観測データを記憶した記憶部４の第１の観測データファイル４２から、画素データ毎の第１の観測データを読み込む（Ｓ１０３）。この観測データは、例えば、植生指数データ、水分指数、中間赤外線バンド、熱赤外線輝度バンド、温度データ等がある。処理部１は、記憶部４の第１の観測データファイル４２から、１ライン又は１カラムのデータを読み込み、読み込まれたライン又はカラムのデータを並列演算計算機のパラメータに設定してもよい。
【００５８】
処理部１は、式（２）を用いて各画素毎に、ある時期の前後所定期間の局所最大値を抽出し、第２の観測データを得て、第２の観測データファイル４３に画素毎に記憶する（Ｓ１０５）。
【００５９】
処理部１は、第２の観測データファイル４３から第２の観測データを読み出し、これに基づき式（１）の関数にフィッティング処理を実行して各係数ｃｋを最小二乗法により決定し、その係数を第１パラメータファイル４６に記憶する（Ｓ１０７）。なお、処理部１は、式（３）で示される赤池の情報量基準ＡＩＣが最小となる周期関数の組み合わせを選択するステップをさらに含むことができる。次のような各ステップが実行される。即ち、
・処理部１は、全ての関数の組み合わせを発生させる第１ステップと、
・処理部１は、各組み合わせに対して最小二乗法を実行してフィッティングに用いられる係数の関数を求める第２ステップと、
・処理部１は、最小二乗法により求めた結果について、それぞれＡＩＣを求める第３ステップと、
・処理部１は、ＡＩＣの最小の組み合わせを使用するフィッティングに用いられる関数として選択する第４ステップと
を含む。
【００６０】
また、処理部１は例えば、ステップＳ１０５の前又は後若しくはステップＳ１０７の前又は後等に、欠測値の補間をしてもよい。さらに欠測値の補間後に、補間修正をしてもよい。さらに、処理部１は、原データとフィッティングされた関数による計算値の差を取り、差がしきい値よりも大きなものを取り除いても良い。また、処理部１は採用すべきデータのみを用いてフィッティング、フィルタリング（カルマンフィルタ）等の各処理を繰り返すことができる。
【００６１】
つぎに、処理部１は、定められた係数ｃｋに従い、関数によりノイズが除去された推定値を第１推定観測データファイル４４に記憶する、又は出力部３に出力若しくは表示する。また、処理部１は、第１パラメータファイル４６から係数を読み出し、カルマンフィルタ処理の初期値を式（１３）のように設定する（Ｓ１０９）。
【００６２】
処理部１は、ｙ（ｔ｜Ｔ）ｏｌｄ とＫの値を、ぞれぞれ初期値として０（ｙ（ｔ｜Ｔ）ｏｌｄ＝０，Ｋ＝０）とする（Ｓ１１１）。つぎに、処理部１は、式（９）に基づいて、カルマンフィルタ処理を行う（Ｓ１１３）。さらに、処理部１は、式（１２）を用いて、データの平滑化する（Ｓ１１５）。
【００６３】
処理部１は、Ｋの値を１繰り上げる（Ｋ＝Ｋ＋１） （Ｓ１１７）。さらに、処理部１は、時期Ｋに、推定値と実際の画素値のずれを求め、この値によって次に行う処理を判断する。即ち、カルマンフィルタ処理が収束してなく、このずれが予め設定したしきい値εよりも大きい場合には、次にステップＳ１２１の処理を行い、一方収束して、しきい値εよりも小さい場合には次にステップＳ１２５の処理を行う（Ｓ１１９）。
【００６４】
即ち、収束していない場合、処理部１は、計算パラメータｃ（１｜０）とＶ（１｜０）の初期値を更新し（ｃ（１｜０）＝ｃ（１｜Ｔ）、Ｖ（１｜０）＝Ｖ（１｜Ｔ））（Ｓ１２１）、さらに、ｙ（ｔ｜Ｔ）ｏｌｄの値をｙ（ｔ｜Ｔ）ｏｌｄ＝ｙ（ｔ｜Ｔ）とする（Ｓ１２３）。一方、収束した場合には、処理部１は、計算結果である推定された係数ｃｉ（ｔ）を第２のパラメータファイル４７に記憶し、ノイズを除去した推定された観測値を第２の推定観測データファイル４５に記憶する、又は、出力部３に出力若しくは表示する（Ｓ１２５）。
【００６５】
５． 実際の衛星データ処理例
ＬＭＦ処理を行うためには、地理補正、センサー補正済みであることが必要である。すべてのシーンで、それぞれの画素の地理座標が揃っている必要がある。そのようなデータは、農林水産研究計算センター衛星データアーカイブシステムであるＳＩＤａＢ(文献１７，１８参照)、Ｐａｔｈｆｉｎｄｅｒ プロジェクト(文献１９参照)、ＵＳＧＳ(文献２０参照)などのｗｅｂサイトからネットワークを通じて入手できる。また、ＳＰＯＴ／ｖｅｇｅｔａｔｉｏｎであれば、切り出し範囲と投影法が揃っていれば、ほとんど前処理なしでＬＭＦ処理することが可能である。
【００６６】
計算上、ＬＭＦ処理を行うには、最低１年分（１０日間合成で３６シーン）が必要である。特に雲の影響が強い熱帯では、１年半（５４シーン）以下では、結果が期待できない。３年以上のデータの蓄積があれば、一ヶ月程度のシーン欠損があってもほぼ問題なく処理することが可能である。計算はかなり膨大なものになるので、並列機などの使用が望ましい。（並列機を用いての処理については付記を参照のこと）
【００６７】
５．１． Ｐａｔｈｆｉｎｄｅｒ データ
図５に、ＬＭＦ処理前後のＮＤＶＩデータ図を示す。この図は、Ｐａｔｈｆｉｎｄｅｒデータ（８ｋｍ ＮＯＡＡ／ＡＶＨＲＲ １０日間合成）を用いてＮＤＶＩのＬＭＦ処理を行った例である。対象画素としては、北緯３６．３１度、東経１３６．９２度（岐阜県大野郡白川村）を選んだ。データ収集期間は９５年１月上旬〜９９年１２月下旬の１８０シーンを用いた。
【００６８】
Ｐａｔｈｆｉｎｄｅｒ データには、画素毎に、データ取得時の雲の有無についてのフラグ（ＣＬＡＶＲ データ）が添付されている。ＣＬＡＶＲの値が、１〜１１ は“ｃｌｏｕｄｙ”、１２〜２２は”ｍｉｘｅｄ”、２３〜３１は”ｃｌｅａｒ”にカテゴリが分かれている(文献２１参照)。雪の影響の無い、９５年５月〜９月のＮＤＶＩデータについて、日本付近を例に取り、ＬＭＦ処理の精度について検証した。
【００６９】
図６に、ＬＭＦ処理前後のデータの相関図を示す。この図は、処理データのＮＤＶＩを横軸にとり、ＬＭＦ処理済みデータを縦軸としたものである。処理前のＮＤＶＩについて０．０１毎に区間を区切り、その区間に属する同じ時期、同一画素のＬＭＦ処理後のＮＤＶＩの平均値と標準偏差をプロットした。雲の影響のない“Ｃｌｅａｒ”カテゴリに属する画素では、ＬＭＦ処理前後の値に高い相関が見られるのに対し、部分的に雲が混ざった“Ｍｉｘｅｄ”や雲に覆われた”Ｃｌｏｕｄｙ”カテゴリに属する画素ではＬＭＦ処理後の値が高くなっている。
【００７０】
雲に覆われた画素のＮＤＶＩの真値を知ることは不可能であるが、少なくとも、雲の影響の無い画素については、ＬＭＦ処理の前後でほぼ値が保存されている。
【００７１】
熱赤外バンドデータの処理：
図７に、ＡＶＨＲＲ／Ｃｈ４のデータをＬＭＦ処理したデータ図を示す。この図によると、ＮＤＶＩの時と同様、雲などの影響が軽減されていることがわかる。
最小値画像：
これまで、ＮＤＶＩ最大値画像は作られてきたが、ＮＤＶＩ最小値画像は雲の影響を受けた画素ばかりを集めてしまうために、意味づけをすることはできなかった。ＬＭＦ処理により雲の影響を除去することではじめて、意味のある最小値画像が得られたことになる。
ここで得られたＮＤＶＩ最小値画像は、年間を通じて植物が多く存在するか否かを示している。また、ＮＤＶＩだけでなく、熱赤外バンドの輝度温度に関しても、ＬＭＦ処理済みデータであれば、最小値の意味づけが可能である。
【００７２】
５．２． ＵＳＧＳ １ｋｍ 全球ＡＶＨＲＲデータ
米国地質調査所（ＵＳＧＳ）にアーカイブされている１ｋｍ解像度の ＡＶＨＲＲ のＮＤＶＩについて、タイ周辺のデータの処理前と処理後のものを比較したところ、タイは１０月〜３月までが乾期、４月〜９月までが雨期であるが、雨期の雲の影響を抑えて植生が繁茂する様子がよく捉えられている。また、データの欠損や原データにおけるパス間の質の違いも吸収できていることがわかる。また、１画素毎に、ＮＤＶＩの最小値と最大値、季節変化パターンが得られるので、植生タイプの分類だけでなく、植生の乾燥度や植生由来の地上可燃物量についての知見を得ることができる。
【００７３】
５．３． 実際のデータ処理例：Ｐａｔｈｆｉｎｄｅｒデータの処理
図８は、ＰａｔｈｆｉｎｄｅｒデータのＮＤＶＩをＬＭＦ−ＫＦ処理した結果を示す図である。ＬＭＦとＬＭＦ−ＫＦの結果を比較すると、ＬＭＦ処理の結果が９５年〜９９年の５年間の平均的な季節変化を表しているのに対し、ＬＭＦ−ＫＦでは、雲の影響を取り除きながら、年々変動にも追従している。注目すべきは、９６年早春は値の立ち上がりが遅く、９８年では逆に値の立ち上がりが早いのがはっきり現れていることである。過去の気象データによると、９６年は北陸地方の積雪量は例年と比較して多かったが、９８年早春はエル・ニーニョ現象によって暖冬となり、積雪量が異常に少なかったことが報告されている。
【００７４】
また、図９は、ＡＶＨＲＲ／Ｃｈ４輝度温度データをＬＭＦ−ＫＦ処理した結果を示す図である。ＮＤＶＩデータと同じ傾向を示しており、ＬＭＦ−ＫＦ では９６年、９８年の特徴がはっきりと現れている。
【００７５】
本実施の形態では、ボルネオ周辺のＮＤＶＩデータについて、９５年〜９９年の５年間のＬＭＦ−ＫＦ処理を行い、１年毎にカラー合成をしたものを求めたところ、９５年〜９９年になるにしたがって、年間を通じてＮＤＶＩの高い場所が減少し、代わって、植生指数の低い赤色の部分が増加している。また、季節変動の大きいことを示す白っぽい部分も増えている。９８年以降、ボルネオ島東部に植生の少ない地域が増えている。これは９８年にエル・ニーニョ現象のため大規模森林火災が発生したが、その延焼地域を表している。
【００７６】
ＬＭＦ処理を用いれば年々変動の解析も可能であるが、データを期間毎に切り出して個別に処理する必要があった。さらには、処理の性質としてデータの収集期間が短くなると結果に雲の影響が含まれやすくなる（モデルが不安定化する）ことや、データ収集期間のつなぎ目での計算値の整合性がとれない場合があった。
【００７７】
しかし、ＬＭＦ−ＫＦでは、長期間のデータを一度に処理しても、それぞれの年の特徴が抽出できることが示された。データの収集期間が長いので、モデルの不安定性が軽減されるためであると考えられる。
【００７８】
６． 時系列データ処理法のまとめ
今回開発した、時系列フィルタと状態空間モデルを用いた、高頻度観測衛星データ処理法では、データ処理期間全体の平均的な季節変化を捉えるだけでなく、雲やヘイズ、システムノイズなどの影響を抑えながら、年々変動の抽出も可能であることがわかった。
【００７９】
この方法により、植生指数だけでなく、熱赤外バンドの輝度温度データについても衛星データから季節変動プロファイルを抽出できるようになったが、これにより、一週間〜１０日毎の時間間隔で全球規模の植生環境の変動を監視することがはじめて可能となった。特に、大規模森林火災や乾燥害などの農業・林業災害は植生の乾燥の程度に影響されるため、植生の状態を迅速かつ的確に把握する手法の開発は防災面で大きな意味を持っている。
【００８０】
植生指数や輝度温度の季節変動データから、森林火災危険度を算出できるが、それについては次節で述べる。
また、植生指数の季節変動パターンの分類によって、広域での土地被覆分類も可能である。
【００８１】
ＬＭＦおよび ＬＭＦ−ＫＦは、画素毎に処理が独立しているため、並列機に適している。したがって、ＭＯＤＩＳなど、従来より解像度の高い高頻度観測衛星データについても、時間をかけずに処理が行える利点がある。
【００８２】
７． 輝度温度データを用いた火災検出
つぎに、ＬＭＦおよびＬＭＦ−ＫＦによって処理されたデータの活用、特に輝度温度データを用いた火災検出の高精度化について言及する。
ＤＭＳＰ／ＯＬＳでの雲域の自動推定：
【００８３】
ＤＭＳＰ／ＯＬＳ での雲域除去は、特定のしきい値を定めるか、ＮＯＡＡ／ＮＧＤＣから配信される全球地表面温度データと比較することで行っている。前者の方法では、しきい値の設定に恣意性があったり、雲除去が適切に行われないなど、精度において問題がある。一方、配信される全球地表面温度データは、配信時刻が衛星データと同時でなく１日程度遅れており、リアルタイム処理には間に合わない。また、地上観測での地表面温度と衛星で観測される熱赤外バンドの
輝度温度は必ずしも同じ物理量ではない。
【００８４】
ＤＭＳＰ／ＯＬＳの熱赤外輝度温度をＬＭＦ処理しておけば、各時期、各画素での雲の影響のない値が得られるので、雲の判定は単純にしきい値を設定すれば足りる。また、観測波長の問題はあるが、次節に示す方法により、熱によっても火災検出が可能になる可能性がある。
ＮＯＡＡ／ＡＶＨＲＲ での火災検出アルゴリズムの単純化：
【００８５】
現在、ＮＯＡＡ／ＡＶＨＲＲ でのホットスポット検出は、単純に輝度温度にしきい値を設定するものや、バンド間演算をしたもので判断するものが多い(文献２４参照)。この方法では、どうしても、裸地や火口、人口構造物（石油および製鉄プラント）などの火災でないものも検出してしまう欠点がある。これらを回避するため、複雑な火災検出アルゴリズムも開発されているが２４）、演算量の点から、リアルタイム処理には適していない。さらには、土地被覆によって火災検出に用いるしきい値だけでなく、アルゴリズムも変更する必要があるとの指摘もなされている。これを実現することを考えた場合、特定の狭い範囲で火災検出システムを構築することは可能であっても、国レベルの広域を対象にすることは困難である。
【００８６】
ここで、ＬＭＦ処理済みの、雲無しの輝度温度データがあれば、それらは、対象となる画素の土地被覆や緯度、標高などの影響がすでに含まれているので、単純にしきい値を設定するだけでも火災精度の向上が期待できる。また、次節以下で述べる方法により、画素内の火災割合を算出することも可能となる。
衛星データによる小規模火災の検出：
【００８７】
Ｐｒｉｎｃｅ らが、静止気象衛星 ＧＯＥＳ を用いてブラジルの森林火災監視を行っている(文献２５参照)。一般に、静止気象衛星は高度が非常に高いため、解像度は低くなってしまうが、画素内の火災割合を算出することで、ＮＯＡＡ／ＡＶＨＲＲ と同程度の火災検知を可能にしている。ここで、熱赤外バンドｉでの、温度Ｔの黒体の分光放射輝度をＬｉ（Ｔ）、画素内の火災割合をｐとすると、
Li(Tobs)＝ｐ・Li(Tfire)＋(1−ｐ)・Lｉ(Tｂ) （１５）
の式が成り立つ。ただし、Ｔｏｂｓは観測されたものの温度、Ｔｆｉｒｅは火災温度、Ｔｂは画素内の火災のない部分の温度である。この式には、未知数がｐ、Ｔｆｉｒｅ、Ｔｂの３個あるので、熱赤外バンドは３つ必要である。しかし、Ｐｒｉｎｃｅらは、Ｔｆｉｒｅ＝４００（Ｋ）、Ｔｂを周囲の比較的低温な画素の平均と仮定して火災割合を算出している。しかし、これでは、火災のない場合、対象となる画素と、その周囲の画素の温度が同じであることを仮定している。換言すれば、対象となる画素と、その周囲の画素の土地被覆や標高が同一であることになるが、この仮定が成り立つ場所はかなり限られてくる。ここで、Ｔｂとして、対象となる画素の同時期のＬＭＦ処理された値を使えば、ＤＭＳＰ／ＯＬＳのように熱赤外バンドが一つしかなくても、火災割合ｐを算出することが可能になる。また、ＮＯＡＡ／ＡＶＨＲＲのように、複数の熱赤外バンドがあれば、火災割合と火災温度の両方を算出することが可能である。
【００８８】
ＮＯＡＡ／ＡＶＨＲＲでは、センサーが飽和してしまう温度が低いため、火災温度の推定は困難であるが、ＭＯＤＩＳ等の新型センサーでは、火災温度による消火優先度の設定なども考えられる。
この方法は、ＧＭＳ−５の後継機（ＭＴＳＡＴ）のデータにも応用が可能であり、毎時の火災情報提供への応用が期待できる。
【００８９】
８．付記及び参考文献
本発明の地球観測衛星データのノイズ除去方法又は地球観測衛星データのノイズ除去装置・システムは、その各手順をコンピュータに実行させるための地球観測衛星データのノイズ除去処理プログラム、地球観測衛星データのノイズ除去処理プログラムを記録したコンピュータ読み取り可能な記録媒体、地球観測衛星データのノイズ除去処理プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
以下に、参考文献を示す。
【００９０】
1) インドネシア森林火災予防プロジェクト、“インドネシア国立公園森林火災調査報告書”、国際協力事業団、1998
2) 沢田治雄、“大規模森林火災へのリモートセンシングの利用”、森林科学、24、22-28、1998
3) 本多嘉明、村井俊治、加藤喜久雄、“世界植生モニタリング」、日本写真測量学会誌「写真測量とリモートセンシング」、31、1、4-14、1992
4) Defries,R. S. and J. R. G. Townshed, “NDVI-derived land cover classifications at a global scale”, Int. J. Remote Sens., 15, 3567-3586, 1994.
5) Myneni, R. B., C. J. Tucker, G. Asrar and C. D. Keeling,“Interannual variations in satellite-sensed vegetation index data from 1981 to 1991”, J. Geophys. Res., 103, D6, 6145-6160, 1998.
6) Los, S. O., C. O. Justice and C. J. Tucker, “A global 1 deg. by 1 deg. NDVI data set for climate studies derived from the GIMMS continental NDVI Data.”, Int. J. Remote Sens., 15, 3493-3518, 1994.
7) GLCC web site : http://edcdaac.usgs.gov/glcc/globdoc1_2.html
8) Viovy, N., O. Arino and A. S. Belward, “The Best Index Slope Extraction (BISE): A method for reducing noise in NDVI time-series”, , Int. J. Remote Sens., 13, 1585-1590, 1992.
9) 李 雲慶、大沼一彦、安田嘉純、“時系列植生指数データのフーリエ解析”、日本写真測量学会誌「写真測量とリモートセンシング」、31、4、4-14、1992
10) Sellers, P. J., C. J. Tucker, G. J. Collatz, S. O. Los, C. O. Justice, D. A. Dazlich and D. A. Randall, “A global 1 deg. by 1 deg. NDVI data set for climate studies. Part 2: The generation of global fields of terrestrial biophysical parameters from the NDVI”, Int. J. Remote Sens., 15, 3519-3545, 1994.
11) Verhoef, W., M. Menenti and S. Azzali, “A colour composite of NOAA-AVHRR-NDVI based on time series analysis (1981-1992)”, Int. J. Remote Sens., 17231-225, 1996.
12) Ricotta, C., G. Avena and A. De Palma, “Mapping and monitoring net primary productivity with AVHRR NDVI time-series: statistical equivalence of cumulative vegetation indices”, ISPRS J. Photogramm. Remote Sens., 54, 325-331, 1999.
13) Azzali, S. and M. Menenti, “Mapping vegetation-soil-climate complexes in southern Africa using temporal Fourier analysis of NOAA-AVHRR NDVI data”, , Int. J. Remote Sens., 21, 973-996, 2000.
14) Roerink, G. J., M. Menenti and W. Verhoef, “Reconstructing cloudfree NDVI composites using Fourier analysis of time series”, Int. J. Remote Sens., 21, 1911-1917, 2000.
15) 澤田義人、沢田治雄、斎藤英樹、“高頻度観測衛星による植生の季節変動観測”、日本写真測量学会講演論文集、2000秋季、209-212、2000
16) 北川源四郎、“FORTRAN77 時系列解析プログラミング”、岩波書店、第５章、1993
17) http://www.affrc.go.jp/Agropedia/
18) 児玉正文、宋 献方、“農林水産リモートセンシングデータベースシステムの構築”、日本リモートセンシング学会第２８回学術講演会論文集、259-260、2000
19) Pathfinder web site:
http://daac.gsfc.nasa.gov/CAMPAIGN_DOCS/LAND_BIO/GLBDST_main.html
20) http://edcdaac.usgs.gov/1KM/1kmhomepage.html
21) Agbu, P. A. and M. E. James, “NOAA/NASA Pathfinder AVHRR Land Data Set User's Manual. Goddard Distributed Active Archive Center”, NASA Goddard Space Flight Center, Greenbelt, 1994.
http://daac.gsfc.nasa.gov/CAMPAIGN_DOCS/LAND_BIO/AVHRR_GD.pdf
22) 北川源四郎、“FORTRAN77 時系列解析プログラミング”、岩波書店、第9章、1993
23) 西山清、“周波数が既知である成分のみを含む時変スペクトルの一推定法”、電子情報通信学会論文誌、J74-A、6、916-918、1991
24) Martin, M. P., P. Ceccato, S. Flasse and I. Downey, “Fire detection and fire growth monitoring using satellite data”, Remote Sensing of Large Wildfires in the European Mediterranean Basin, E. Chuvieco (ed.), Splinger-Verlag, chapt. 6, 1999.
25) Prins, E. M. and W. P. Menzel, “Geostationary satellite detection of biomass burning is South America”, Int. J. Remote Sens., 13, 2783-2799, 1992.
【００９１】
【発明の効果】
本発明によると、高頻度観測衛星からのデータの画素毎の時系列的な特徴を保存し、時間分解能を完全に保持したまま、雲などの影響およびシステムノイズの除去をすることができる。
また、本発明によると、時間分解能を保持したまま植生などの変動の様子が捉えられるため、防災、環境、農業分野での広域現状把握、および予測が容易となる。さらに、本発明によると、適切な関数を自動的に選択し、ノイズ除去処理を自動で行えるものであり、定常業務の一部分とすることが可能である。さらに本発明によると、水分指数など他の地球高頻度観測衛星のデータにもそのまま転用および応用ができる。
【図面の簡単な説明】
【図１】シミュレーションによるフィルタ特性の評価図
【図２】シミュレーション結果の例を示す図
【図３】命令列スケジュール処理を実行するための装置の構成図
【図４】データ処理のフローチャート
【図５】ＬＭＦ処理前後のＮＤＶＩデータ図
【図６】ＬＭＦ処理前後のデータの相関図
【図７】ＡＶＨＲＲ／Ｃｈ４のデータをＬＭＦ処理したデータ図
【図８】ＰａｔｈｆｉｎｄｅｒデータのＮＤＶＩをＬＭＦ−ＫＦ処理した結果を示す図
【図９】ＡＶＨＲＲ／Ｃｈ４輝度温度データをＬＭＦ−ＫＦ処理した結果を示す図
【符号の説明】
１ 処理部
２ 入力部
３ 出力部
４ 記憶部
４１ 計算パラメータファイル
４２ 第１の観測データファイル
４３ 第２の観測データファイル
４４ 第１の推定観測データファイル
４５ 第２の推定観測データファイル
４６ 第１のパラメータファイル
４７ 第２のパラメータファイル
Claims (11)
- 処理部は、観測対象エリア内の各画素に対応して、衛星による植生の季節変動についての時系列の第１の観測データを記憶した第１の観測データファイルから、画素毎の第１の観測データを読み込むステップと、
処理部は、各画素毎に、ある時期又は時点に対応して、該時期又は時点の前後所定期間の局所最大値を抽出することにより時系列の第２の観測データを求め、それを各画素に対応して第２の観測データファイルに記憶するステップと、
処理部は、第２の観測データファイルに記憶された第２の観測データを、長期的非周期項及び季節変化項を含む次式の関数を用い、各第１の係数ｃｉを最小二乗法又は他の推定法により決定することでフィッティングし、各第１の係数ｃｉを第１のパラメータファイルに記憶するステップと、
ｃｉ：求める第１の係数
ｔ：時期又は時間
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列（例：１，２，３，４，６，１２）
Ｍ：１年間のデータ数（シーン数） ］
処理部は、第１の係数ｃｉで定められた前記関数に従い、第２の観測データから雑音データを除去した第１の推定観測データを求め、それを第１の推定観測データファイルに記憶する、又は、出力部に出力若しくは表示するステップと
を含む地球観測衛星データのノイズ除去方法。 - 処理部は、第２の観測データについて次式の時系列の状態空間モデルを適用し、該状態空間モデルの各第２の係数ｃｉ（ｔ）を推定する処理を実行し、各第２の係数ｃｉ（ｔ）を第２のパラメータファイルに記憶するステップと、
ｃｉ（ｔ）：求める第２の係数
ｗｋ（ｔ）：システム雑音
ｙ（ｔ）：観測される値
ｖ（ｔ）：平均０、分散がσｖの白色雑音（観測雑音）
処理部は、第２の係数ｃｉ（ｔ）で定められた前記関数に従い、第２の観測データから雑音データを除去した第２の推定観測データを求め、それを第２の推定観測データファイルに記憶する、又は出力部に出力、若しくは表示するステップ
をさらに含む請求項１に記載の地球観測衛星データのノイズ除去方法。 - 時刻ｋに対応して、時期ｋ−ｎから時期ｋまでのデータについて最大値を抽出するステップと、
時刻ｋに対応して、時期ｋから時期ｋ＋ｎのデータについて最大値を抽出するステップと、
時刻ｋに対応して、得られたこれら二つの最大値のうち、値の小さいものを局所的最大値とするステップと、
を含む請求項１又は２に記載の地球観測衛星データのノイズ除去方法。 - 時期ｋに、推定値と実際の画素値のずれが、設定したしきい値よりも大きい場合には雲やシステムノイズの影響とみなしてその時期の観測データは採用しないようにして、採用されたデータだけを用いて再度、係数及び／又は推定観測データを求めることを特徴とする請求項１又は２に記載の地球観測衛星データのノイズ除去方法。
- 前記第１の観測データは、植生指数データ、水分指数、中間赤外線バンド、熱赤外線輝度バンド、温度データのいずれか又は複数であることを特徴とする請求項１又は２に記載の地球観測衛星データのノイズ除去方法。
- 前記雑音データは、雲、霧、雨、ヘイズ等の大気状態に起因するもの、及び／又は、センサ感度等のシステムに起因するものを含むことを特徴とする請求項１又は２に記載の地球観測衛星データのノイズ除去方法。
- 処理部は、観測データおよび推定データの複数を組み合わせて出力部に出力又は表示するステップをさらに含む請求項１又は２に記載の地球観測衛星データのノイズ除去方法。
- 処理部は、観測対象エリア内の各画素に対応して、衛星による植生の季節変動についての時系列の第１の観測データを記憶した第１の観測データファイルから、画素毎の第１の観測データを読み込むステップと、
処理部は、各画素毎に、ある時期又は時点に対応して、該時期又は時点の前後所定期間の局所最大値を抽出することにより時系列の第２の観測データを求め、それを各画素に対応して第２の観測データファイルに記憶するステップと、
処理部は、第２の観測データファイルに記憶された第２の観測データを、長期的非周期項及び季節変化項を含む次式の関数を用い、各第１の係数ｃｉを最小二乗法又は他の推定法により決定することでフィッティングし、各第１の係数ｃｉを第１のパラメータファイルに記憶するステップと、
ｃｉ：求める第１の係数
ｔ：時期又は時間
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列（例：１，２，３，４，６，１２）
Ｍ：１年間のデータ数（シーン数） ］
処理部は、第１の係数ｃｉで定められた前記関数に従い、第２の観測データから雑音データを除去した第１の推定観測データを求め、それを第１の推定観測データファイルに記憶する、又は出力部に出力若しくは表示するステップ
をコンピュータに実行させるための地球観測衛星データのノイズ除去処理プログラム。 - 処理部は、第２の観測データについて次式の時系列の状態空間モデルを適用し、該状態空間モデルの各第２の係数ｃｉ（ｔ）を推定する処理を実行し、各第２の係数ｃｉ（ｔ）を第２のパラメータファイルに記憶するステップと、
ｃｉ（ｔ）：求める第２の係数
ｗｋ（ｔ）：システム雑音
ｙ（ｔ）：観測される値
ｖ（ｔ）：平均０、分散がσｖの白色雑音（観測雑音）
処理部は、第２の係数ｃｉ（ｔ）で定められた前記関数に従い、第２の観測データから雑音データを除去した第２の推定観測データを求め、それを第２の推定観測データファイルに記憶する、又は出力部に出力若しくは表示するステップ
をさらにコンピュータに実行させるための請求項８に記載の地球観測衛星データのノイズ除去処理プログラム。 - 処理部は、観測対象エリア内の各画素に対応して、衛星による植生の季節変動についての時系列の第１の観測データを記憶した第１の観測データファイルから、画素毎の第１の観測データを読み込むステップと、
処理部は、各画素毎に、ある時期又は時点に対応して、該時期又は時点の前後所定期間の局所最大値を抽出することにより時系列の第２の観測データを求め、それを各画素に対応して第２の観測データファイルに記憶するステップと、
処理部は、第２の観測データファイルに記憶された第２の観測データを、長期的非周期項及び季節変化項を含む次式の関数を用い、各第１の係数ｃｉを最小二乗法又は他の推定法により決定することでフィッティングし、各第１の係数ｃｉを第１のパラメータファイルに記憶するステップと、
ｃｉ：求める第１の係数
ｔ：時期又は時間
Ｎ：周期関数のペア数（この場合ｓｉｎとｃｏｓのペア数）：
ｌ：周期関数の仮の番号
ｋｌ：１２ヶ月が回帰の整数倍になる数列（例：１，２，３，４，６，１２）
Ｍ：１年間のデータ数（シーン数） ］
処理部は、第１の係数ｃｉで定められた前記関数に従い、第２の観測データから雑音データを除去した第１の推定観測データを求め、それを第１の推定観測データファイルに記憶する、又は出力部に出力、若しくは表示するステップ
をコンピュータに実行させるための地球観測衛星データのノイズ除去処理プログラムを記録したコンピュータ読み取り可能な記録媒体。 - 処理部は、第２の観測データについて次式の時系列の状態空間モデルを適用し、該状態空間モデルの各第２の係数ｃｉ（ｔ）を推定する処理を実行し、各第２の係数ｃｉ（ｔ）を第２のパラメータファイルに記憶するステップと、
ｃｉ（ｔ）：求める第２の係数
Ｗｋ（ｔ）：システム雑音
ｙ（ｔ）：観測される値
ｖ（ｔ）：平均０、分散がσｖの白色雑音（観測雑音）
処理部は、第２の係数ｃｉ（ｔ）で定められた前記関数に従い、第２の観測データから雑音データを除去した第２の推定観測データを求め、それを第２の推定観測データファイルに記憶する、又は出力部に出力、若しくは表示するステップ
をさらにコンピュータに実行させるための地球観測衛星データのノイズ除去処理プログラムを記録した請求項１０に記載のコンピュータ読み取り可能な記録媒体。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002103966A JP4003869B2 (ja) | 2002-04-05 | 2002-04-05 | 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002103966A JP4003869B2 (ja) | 2002-04-05 | 2002-04-05 | 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2003296702A JP2003296702A (ja) | 2003-10-17 |
| JP2003296702A5 JP2003296702A5 (ja) | 2005-06-16 |
| JP4003869B2 true JP4003869B2 (ja) | 2007-11-07 |
Family
ID=29389479
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002103966A Expired - Fee Related JP4003869B2 (ja) | 2002-04-05 | 2002-04-05 | 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4003869B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11933673B2 (en) | 2018-02-19 | 2024-03-19 | Ihi Corporation | Heat source detection device |
Families Citing this family (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2008281505A (ja) * | 2007-05-14 | 2008-11-20 | Japan Aerospace Exploration Agency | ３次元離散データのリサンプル方法、及び装置 |
| KR101089220B1 (ko) | 2010-02-18 | 2011-12-02 | 공주대학교 산학협력단 | 시공간 연속성을 이용한 식생 지수 보정방법 |
| JP5647814B2 (ja) * | 2010-05-19 | 2015-01-07 | 日本電産エレシス株式会社 | 電子走査型レーダ装置、受信波方向推定方法及び受信波方向推定プログラム |
| JP5600499B2 (ja) * | 2010-07-01 | 2014-10-01 | 日本電産エレシス株式会社 | 電子走査型レーダ装置、受信波方向推定方法及び受信波方向推定プログラム |
| KR101404430B1 (ko) | 2013-06-11 | 2014-06-10 | 서울시립대학교 산학협력단 | 적외선 영상을 이용한 지표온도감률 추정 방법 |
| CN104251662B (zh) * | 2013-06-27 | 2017-10-31 | 杭州中科天维科技有限公司 | 有序点云阈值自适应抑噪技术 |
| CN104502241B (zh) * | 2014-05-14 | 2017-04-05 | 淮海工学院 | 基于gnss的雾霾监测与分析技术 |
| CN107888878B (zh) * | 2017-11-16 | 2018-11-16 | 曹欢欢 | 火灾自动定位式照相平台控制方法 |
| JP7036704B2 (ja) * | 2018-11-16 | 2022-03-15 | ヤフー株式会社 | 情報処理装置、情報処理方法、及び情報処理プログラム |
| CN110378858B (zh) | 2019-07-04 | 2021-04-09 | 浙江大学 | 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 |
| CN111553237B (zh) * | 2020-04-23 | 2023-11-07 | 江西理工大学 | 一种基于多态叠加Gamma分布的LJ1-01夜间灯光数据去噪方法 |
| CN111650102B (zh) * | 2020-05-26 | 2023-08-29 | 北京中科锐景科技有限公司 | 基于卫星数据的霾污染解析方法、装置、介质及设备 |
| WO2022080418A1 (ja) * | 2020-10-14 | 2022-04-21 | 株式会社ブリヂストン | 落葉検出システムおよび落葉検出方法 |
| KR102502154B1 (ko) * | 2022-07-08 | 2023-02-22 | 에스케이임업 주식회사 | 이산화탄소 포집량 예측 시스템 |
| CN117668468B (zh) * | 2024-01-31 | 2024-04-26 | 山东省舜天化工集团有限公司 | 一种化工制备数据智能分析管理系统 |
-
2002
- 2002-04-05 JP JP2002103966A patent/JP4003869B2/ja not_active Expired - Fee Related
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11933673B2 (en) | 2018-02-19 | 2024-03-19 | Ihi Corporation | Heat source detection device |
Also Published As
| Publication number | Publication date |
|---|---|
| JP2003296702A (ja) | 2003-10-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4003869B2 (ja) | 地球観測衛星データのノイズ除去処理方法、ノイズ除去処理プログラム、ノイズ除去処理プログラムを記録した記録媒体 | |
| Rozenstein et al. | Estimating cotton water consumption using a time series of Sentinel-2 imagery | |
| Santos et al. | Savanna and tropical rainforest biomass estimation and spatialization using JERS-1 data | |
| Gao et al. | Toward mapping crop progress at field scales through fusion of Landsat and MODIS imagery | |
| Dixon et al. | Satellite prediction of forest flowering phenology | |
| Ballesteros et al. | Characterization of Vitis vinifera L. canopy using unmanned aerial vehicle-based remote sensing and photogrammetry techniques | |
| Sepulcre-Canto et al. | Assessment of the EUMETSAT LSA-SAF evapotranspiration product for drought monitoring in Europe | |
| Fernandez-Carrillo et al. | Estimating prescribed fire impacts and post-fire tree survival in eucalyptus forests of Western Australia with L-band SAR data | |
| Peng et al. | Relationships between remote-sensing-based agricultural drought indicators and root zone soil moisture: a comparative study of Iowa | |
| Tokura et al. | Long-term variability in vegetation productivity in relation to rainfall, herbivory and fire in Tswalu Kalahari Reserve | |
| Yu et al. | Satellite observations of the seasonal vegetation growth in central Asia: 1982-1990 | |
| Awada et al. | A remote sensing and modeling integrated approach for constructing continuous time series of daily actual evapotranspiration | |
| González-Dugo et al. | Long-term water stress and drought assessment of Mediterranean oak savanna vegetation using thermal remote sensing | |
| Zang et al. | Spatially-explicit mapping annual oil palm heights in peninsular Malaysia combining ICESat-2 and stand age data | |
| Pappalardo et al. | Performance evaluation of a low-cost thermal camera for citrus water status estimation | |
| Reyes-Gonzalez | Using remote sensing to estimate crop water use to improve irrigation water management | |
| Dineshkumar et al. | Rice monitoring using sentinel-1 data in the google earth engine platform | |
| Etchanchu et al. | On the use of high resolution satellite imagery to estimate irrigation volumes and its impact in land surface modeling | |
| Fors | Improved understanding of water balance in the Malwathu Oya river basin using SWAT and remote sensing | |
| Powers et al. | Grand Valley Ecological Forecasting: Assessing Trends in Pinyon-Juniper Habitat Relative to Drought, Beetle Infestation, Wildland Fires, and Treatment to Plan Future Management Strategies | |
| Domingo et al. | Integration of climate time series and MODIS data as an analysis tool for forest drought detection | |
| Khare et al. | Phenology analysis of moist decedous forest using time series Landsat-8 remote sensing data | |
| Elhag et al. | Impact of climate variability on vegetative cover in the Butana area of Sudan | |
| Ning et al. | Predicting hydrological response to forest changes by simple statistical models: the selection of the best indicator of forest changes with a hydrological perspective | |
| Stone et al. | Grand Valley Ecological Forecasting |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A712 Effective date: 20031031 |
|
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20040129 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040922 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040922 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20070711 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20070724 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20070816 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100831 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110831 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110831 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120831 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130831 Year of fee payment: 6 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| LAPS | Cancellation because of no payment of annual fees |