平成 18 年革新的農業技術習得研修

リモートセンシングによる水稲生育ステージの評価法
独立行政法人農業環境技術研究所

生態系計測研究領域

坂本利弘

1. はじめに
近年、稲作の作期早進・夏季気温の高温化傾向等から、登熟期の高温傾向による全国的な米の品質低下が問題になっている。また、地球温暖化・異常気象の多発による水稲被害についてもその関心が高まっている。このような気象条件の変化に対する減収リスクを回避する手段の一つとして、栽培時期の移動が指摘されていが、現地調査データや統計データ等を除けば、現状の水稲栽培時期がどのような空間的分布構造を詳細に知るすべは無い。そのため、気象要因の時空間的な変動特性を考慮した稲作生産の影響評価や作期移動によるリスク回避を評価する上で必要不可欠な水稲生育ステージ・栽培時期の空間分布を広域的に把握する手法の開発が重要であると考える。
地球観測衛星は、様々な種類の光学センサー・マイクロ波センサーによって、広域の地表面状態を均質的かつ継続的に観察しており、農業現場における衛星画像の利用事例も数多く報告されている[1]。本課題で紹介する水稲生育ステージ・栽培時期の把握は、分光反射率から算出される植生指数の時間変化を用いているため、生育期間中における複数枚の光学衛星画像の取得が前提となっている。1972 年からアメリカ NASA と USGS によって運用されてきた Landsat シリーズは、その空間分解能(MSS:80m, TM,ETM+:30m)から、農作物・自然植生の季節変化の把握に利用されてきた、その回帰日数(Landsat-1〜3: 18 日, Landsat-4,5,7: 16 日)による観測頻度の少なさと雲被覆による欠測から、生育ステージを時間的に把握することは困難であった。また、フランスの地球観測衛星 SPOT[2]は、現在３機（SPOT-2,4,5）運用体制であることとそのポインティング機能によって高解像度(10〜20m)でかつ高頻度に農耕地を観測することが可能である[2]。しかしながら、画像取得費用が高額(SPOT5, 20m color, full scene, １枚の場合: アーカイブ：¥294,000, プログラミング: ¥420,000 [2006 年 10 月時点])であることが生育ステージ把握への利用を困難にしていた。そのため、これまでの光学衛星画像を用いた植生生育の時系列モニタリン
グ研究は、日データを無料利用することが可能な NOAA/AVHRR センサー, SPOT/VEGETATION
センサーから得られる植生指数を用いたものが数多かった[3]。しかし、これらの衛星センサーは、空間分解能は約 1km と粗いことが欠点であり、日本のように小規模な水田圃場を観測するには、隣接する土地被覆のミクセル影響が問題であった。そして、NASA の MODIS センサーの運用開始によって、世界各地の中分解能(250, 500m)日々画像を無料で利用できるようになり、農耕地における水稲生育の時間変化をより精緻に把握することが可能となった。本課題では、MODIS 時系列データを用いた水稲生育ステージ評価手法の説明を行なうと共にその適用事例について紹介する。
2. MODISMODIS (MOderate Resolution Imaging Spectroradiometer ) は、NASA の地球観測衛星 EOS
AM-1 (Terra)と PM-1 (Aqua)に搭載された光学センサーであり、地表面、海洋、大気、災害監視を含む多様な観測用途に利用されている。MODIS の諸元と Band1〜7 の観測波長帯を表１,図表１. MODIS の諸元

1 に示す[4]。
表１. MODIS の諸元

軌道

太陽同期準回帰軌道 (高度: 705km)
Terra 午前 10 時 30 分ごろ北から南へ通過
Aqua 午後 1 時 30 分ごろ南から北へ通過

回帰日数

16 日

観測幅

2330km (±55 度)

量子化数12bits

空間解像度 250ｍ

/バンド数

500ｍ

1000ｍ

5(Bands29(Bands2（Bands 1-2) 3-7)8-36)

図１. MODIS(Band1-7), NOAA/AVHRR,
SPOT/VEGETATION の観測波長帯

打上年月日 Terra: 99 年 12 月 18 日
Aqua: 02 年 5 月 4 日
設計寿命6年

MODIS データは、多くの受信機関からインターネットを通して無料配布されているが、地表面反射率プロダクト、雲除去コンポジットデータ等の提供されている高次プロダクトの種類・画像取得期間・解像度・フォーマット等は、配信機関によって異なる[5]。NASA EOS Data Gateway[6]は、2000 年以降の高次プロダクトを配信している点とシステム補正に関する情報が公開されている点で利便性が高く、本課題においても利用しているが、配信されている MODIS 高次プロダク
トの地図投影法が、Integerized Sinusoidal(ISIN)であるため、緯度に応じて経度方向に画像がリサンプリングされ、日本(九州-北海道： 北緯 29〜45.0 度)を対象として場合、本来取得されたデータ量の約 13〜30%が失われている欠点がある(図 2)。したがって、空間分解能が重要になるよ
うな用途については、日本の受信局から配信されているデータを利用することが望ましい。
本課題で用いた MODIS プロダクトは、すでに大気補正が行われている地表面反射率プロダクト(解像度250, 500m)である。500m 解像度の画像については、8日間コンポジットデータも利用している。本コンポジットデータは、観測条件が良く、雲の影響が少ない画素を8 日間ごとに合成した画像であり、あらかじめ雲影響が除かれ、さらにデータ容量が少なくなっている点で、日データよりも扱い易い。

図 ２ . Integerized Sinusoidal(ISIN) 図 法 のMODIS 画像で見た中部および関東地方

3. 植生指数
光学センサー画像で最も良く用いられる植生指数は、NDVI(Normalized Difference VegetationIndex)である。NDVI は、植物体によって反射される近赤外域とクロロフィルが良く吸収する赤色域の反射率の差を和で除したものであり、その値は-1 から 1 にかけて分布する[7]。主に、植物体の量や活性度の高さを示す指標として用いられる。下記に NDVI の計算式を示す。

NDVI =NIR − RED
NIR + RED・ ・・・(1)

NIR: 近赤外領域の反射率(MODIS Band2), RED: 赤色領域の反射率(MODIS Band1)
本課題で用いた植生指数 EVI(Enhanced Vegetation Index)は、赤色域と近赤外域に加えて青色域の反射率を考慮しており、大気や背景土壌の影響を受けにくく、密度の高い植生に対する感度が NDVI よりも良いとされる[8]。下記に本課題で用いた EVI の計算式を示す。

EVI = G *NIR − RED
NIR + C1 * RED − C 2 * BLUE + L

・・・・(2)

NIR: 近赤外領域の反射率(MODIS Band2), RED: 赤色領域の反射率(MODIS Band1) BLUE: 青色領域の反射率(MODIS Band3), G: ゲイン調パラメータ(=2.5), C1: 赤色域大気補正パラメータ(=6),C2:青色域大気補正パラメータ(=7.5), L: 群落背景の輝度補正パラメータ
図 3 に霞ヶ浦南側の水田地帯における 2002 年の NDVI および EVI の時系列プロファイルを示す。全般的に NDVI が EVI よりも高く推移しており、ノイズ状の振幅は、NDVI よりもＥＶＩの方が小さく、また、植生の生育開始前(DOY 150 以前)とピーク時(DOY 200 付近)における植生指数の差は、EVI のほうが NDVI のそれよりも大きいことが観察される。本課題では、大気影響に強いとされる EVI を水稲生育の季節変化を捉える指標として用いていた。

2002 年 6 月 9 日(DOY 160)
R:G:B = Band 6:2:1
図 3. 霞ヶ浦南側における時系列植生指数(EVI, NDVI)プロファイル(500m 解像度)

4. 時系列解析の概要
本課題で用いている MODIS データは、大気補正処理がすでに行われている。しかしながら、厚い雲に遮られてしまえば、地表面からの反射光を観測することができず、大気補正の効果も期待薄で、雲が被覆している日の EVI 値は図Ｘのように低くなる。このような入力データから、雲被覆や観測角度によるミクセル影響のノイズを取り除き、水稲生育に伴う植生指数の季節変化特徴量を抽出するために、青色領域の反射率データによる前処理とウェーブレット変換による時間周波数分解、または、フーリエ変換による周波数分解を適用した。図 4 に 2002 年の滋賀県湖北の水田における 500m 解像度 MODIS データから導いた EVI 時系列データとノイズ除去後の変化曲線を示す(詳しくは、文末引用資料[9]を参考)。図中に示す曲線は、フーリエ変換およびウェーブ
レット変換によってノイズ除去された

0.6EVI の時系列変化を表しており、田植後の0.5EVI の立ち上がりや、出穂日付近のピーク0.4時期をうまく捉えている。
ける極大値を出穂日として、極小値または0.3EVI

ノイズ除去後の EVI 季節変化曲線にお変曲点を田植日、刈取日として、各生育ス0.20.1

テージの期日を推定した[9]。そして、20020年の秋田から鹿児島までの主要な水田地-0.1

帯 30 ヶ所において 500m 解像度 8 日間コ-0.2取日までの期間と定義）と同年の統計デー

60901202703001:1

Error <= ±8days
Error <= ±16days
DOY in Estimation

200
Planting date with Wavelet analysis
Heading season with Wavelet analysis
Harvesting date with Wavelet analysis

160

Planting date with Fourier analysis
Heading date with Fourier analysis

るが、その推定精度（表 2）は、最も結果が良かった場合(ウェーブレット変換, Coiflet

80

Harvesting date with Fourier analysis

80120160

order=4)でも、二乗平均平方誤差で 9.0〜200
Day of year in Statistics

240

280320

図5.統計データにおける田植日・出穂日・刈取日の最盛日と各生育ステージに対応する推定期日との比較

ら１日単位の高い精度が要求される目的に

170

は、すぐに利用可能であるとは言いがたい。

Error <= ±8days

うに、隣接する圃場の生育時期の違いによるミクセルが推定精度に影響している可能性

Error <= ±16days

150
Growing period in Estimation

ものは 2 ヶ月近く)であることからわかるよ

1:1160

検証比較する統計データにおいて田植が地域によって様々(1 週間から 1 ヶ月、長い

140130120110

Using Wavelet analysis （Coiflet 4)
Using Fourier analysis

10090が考えられる。それ以上に、雲被覆による利用可能データが制限されること（8 日間コン
ポジット画像においても、完全に雲が除去さ

360

240

120期・出穂期・刈取期の始期から終期までの間

330

280

期間の違いをうまく捉えることができてい12.1 日程度であり、高温不稔ような数時間か

150 180 210 240
Day of year ( in 2002 )320
平滑化結果（滋賀湖北）

した結果を図 5,6,7 に示す。
地域ごとに特徴的な栽培時期の早晩や生育30

図4. ウェーブレット変換とフーリエ変換による時系列EVIデータのタにおける各時期の最盛期と期間を比較
図 7 に示すように、日本列島の北から南へ、

heading period harvest period
ウェーブレット変換(Coiflet order=4)を用いた時
フーリエ変換を用いた時解析不適と見なされたデータ
時系列解析に用いたデータ

0

ンポジットデータから推定された田植日、出穂日、刈取日、生育期間（田植日から刈planting period

90100110120130140
Growing period in Statistics

150160170

図6. 統計に記録された生育期間（田植日から収穫日までの期間）と推定された期間との比較れることはない。つまり、入力データの時間解像度は 8 日間以上）が、推定精度に強く影響していると考える。いずれにせよ、現行の 10 日間前後の推定精度では、1 週間〜10 日前後の各生育ステージの年変動を把握することが困難で誤差に埋もれてしまう可能性がある。しかしながら、定性的であるにせよ水稲生育ステージを面的に広域把握することができる点では、十分意義深いと考える。今後、Aqua 搭載の MODIS のように上空通過時間の異なる衛星画像の利用や、250m 解像度のデータセットの利用によって、さらなる推定精度の向上が期待される。

図7．ウェーブレット変換、フーリエ変換を利用した場合の推定値と2002年の生育ステージを記録した統計データとの比較 *青：田植期、緑:出穂期、赤：刈取期

表2．各種ウェーブレット変換、フーリエ変換を用いた時の推定誤差（RMSE・日）

5. 作付け体系の空間的把握
前述のように、圃場規模が小さく作期が大幅に年変化しない日本の水田においては、MODIS時系列データによる水稲生育ステージ把握手法のさらなる推定精度の向上が必要とされた。しかしながら、自然環境の変動影響を強く受け、大規模圃場で栽培される海外の農産物の生育モニタリングを対象とした場合では、現行の推定精度でも十分利用可能であると考える。特に、栽培履歴が統計情報として整備されていない発展途上国の栽培時期や作付け体系の空間的把握や我が国の主な輸入農産物を栽培している地域における作物生育の監視等への利用が期待される。ここでは、ベトナム国南部メコンデルタにおける作付け体系・生育ステージ（出穂日）の時空間把握を行った適用事例について紹介する。

5.1 メコンデルタ
ベトナム南部メコンデルタ[Lat. 8.5-11˚N, Long.
104.5-106.8˚E]は、インドシナ半島最南端、メコン川河口部に位置する（図 8）。メコン川は、カンボジアの首都プノンペンを過ぎたあたりからバサック川とメコン川に分岐し、さらなる分流を経て東シナ海に注ぎ込む。また、メコン・バッサク両河川は、毎年、雨期に集中する降雨によってその水位を急上昇させ、本川とカンボジア国境から溢れ出た表層水が洪水となって、メコンデルタの半分以上を浸水させる。メコン川のもたらす豊富な水資源と一年中温暖な気候によって、ベトナム領メコンデルタでは、２期・３期作の水稲栽培が可能である。2003 年におけるメコンデルタの米生産量と単収は、冬春米[849 万ｔ,5.8t/ha],夏秋米[766 万ｔ,4.0t/ha],雨期作米[137 万ｔ,3.6t/ha]であった。また、

図8．研究対象領域 *赤色で囲まれた部分がベトナムメコンデルタ１３州を示す。

ベトナムは、タイに次ぐ世界第２位の米輸出国（380 万ｔsource: FAOSTAT 2003）でもあり、その輸出米の大半は、メコンデルタにおいて生産されたものと言われている。ベトナム南部メコンデルタは、ベトナムにとって米の国際取引を通じた貴重な外貨獲得手段であり、米自給が不安定な国々の食料供給不足を補う役目を果す重要な米生産地帯といえる。
5.2 データと方法
ここで利用したデータは、MODIS/Terra 地表面反射率の 8 日間コンポジットプロダクト（2002年と 2003 年の 2 年分、500ｍ分解能）である。日本の場合と同様に時系列 EVI データにウェーブレット変換を適用し、ノイズ除去を行った時系列 EVI 変化曲線を、水稲生育を表す指標として用いた。また、EVI 季節変化曲線の極大値（EVI≧0.4）を水稲出穂日として推定し、出穂日の出現回数と出現時期の組み合わせによって、水稲の作付体系の分類を行なった。水田分布については、ベトナム南部農業計画設計研究所が作成した 2002 年の土地利用図を用いている(詳しくは、
文末引用資料[10]参考)。
5.3

結果：2002 年・2003 年の推定出穂日の時空間分布・作付体系

2002・2003 年の推定出穂日空間分布を 3 ヶ月ごとに比較した結果、ベトナム南部メコンデルタの上流域・中流域・沿岸域では、その推定出穂日の出現時期と回数が明らかに異なっていることが確認された（図 9）。上流域における出穂日は、1〜6 月にかけて、中流域では、１〜9 月にかけて、沿岸域では、７〜12 月にかけて推定された。3 ヶ月ごとに区分した推定出穂日の出現回数とパターンによって作付体系を分類した結果を図 10 に示す。分類された作付体系は、上流域から沿岸域に広がる扇状に分布しており、この特徴的な分布は、雨期（8〜11 月）に中上流域で発生する洪水と乾期（３〜4 月）に沿岸域で発生する塩水遡上によって説明することができる。また、上流域における 2003 年の洪水後(1〜9 月)の推定出穂日は、2002 年の推定出穂日よりも 20〜30日早く推定されていた。これは、2002 年雨季に発生した洪水のほうが 2001 年よりもその規模が小さく、直播に適度な水位にまで減水した時期が2003 年乾季の方が早かったためと考えられる。
Chau Doc Station におけるメコン川水位の日デー
タからも、2002 年の減水時期が 2001 年よりも早かったことを示していた。
5.4

結果: 2002 年と 2003 年の作付体系の変化

2002 年と 2003 年の作付体系を比較した結果、2003 年の AnGiang 省南部（図１0α）と SocTrang
省東部（図１0β）において、作付回数が２回から３回に増加している地域を確認することができた。
図１0αに示す地域は、堤防を建設し輪中を形成することによって雨期作が可能になったことを意味している。また、沿岸域にある図１0βの作付回数の増加は、乾期作が可能な地域が拡大していることを意味している。図１0γの示す TienGiang 省
GoCong 地区では、堤防と樋門制御により、乾期

図9．2002年・2003年のMODISデータから３ヶ月ごとの灌漑水塩分濃度上昇を抑え、三期作を行なっての推定出穂日の空間分布図いることが知られている。よって、図１0βの示す SocTrang 省東部についても、GoCong 地区と同様に何らかの塩水遡上対策が施されたことによって、乾期作を含む三期作が可能になったと言える。このように、水資源の量的・質的な季節変化をインフラ整備によって克服し、作付け回数を 2 回から 3 回へ増加させることによる大規模な食糧増産を面的に確認することができた。
2002 年

α2003 年

γ

作付体系

β

三期作
乾季二期作
雨季二期作
乾季・雨季二期作
乾季一作
雨季一作

図 10．2002・2003 年における作付体系分類図6. おわりに本研修で紹介する手法は、日本の水稲生育ステージを定性的に把握・評価することを目的としているが、現状の精度では、日単位の精度が必要とされる場面には対応できていない。

今後、推定精度・ミクセル等の数々の問題を克服するべく、さらなる手法の開発・改良を引き続き行っていく必要がある。しかし、日本において空間解像度、時間推定精度点で問題を残した本手法も、ベトナム南部メコンデルタの適用事例に示したように、その適用場所と利用目的を変えれば、有益な情報を導きだすことができる。メコンデルタのように広域の作付け時期の違いを定性的であれ一目瞭然に表すことができるので、限りある水資源の分配方法・利用計画についての流域住民の合意形成に役立てることが期待される。また、海外の主要な食料生産地域の環境変化や作物の生育状態を広域に監視することは、食料の 60%（カロリーベース）を海外輸入に依存している我が国にとっても重要な意義があると考える。最近
では、2004 年に台湾によって打ち上げられた FORMASAT-2[11]は、毎日同じ軌道を通ることから、観測幅は 24km と狭いものの高解像度(マルチスペクトル 8m, パンクロ 2m)で同じ場所を毎日撮影することが可能であり、日本の生産現場における実利用も期待されている。残念ながら、FORMASAT-2 データも冒頭で述べた SPOT(HRV-XS, HRVIR-X, HRG-X)と同様に有償であり、本研修で紹介したような時系列解析を行なうには、莫大な研究予算を準備する必要がある。米国から軍事用に打ち上げられ、後に民間に開放された GPS 衛星のように、今後、各国政府の所有する非民生用の高解像度衛星データが、温暖化研究と平和目的のために開放され、無料で多時期の高解像度画像が解析可能になることを期待したい。
7. 参考文献・Web サイト・出典[1] 財団法人

リモート・センシング技術センター

http://www.restec.or.jp/
農業リモートセンシング推進委員会 http://www.restec.or.jp/eeoc/agri/agri.html[2] 東京スポットイマージュ

http://www.spotimage.co.jp/html/_511_.php[3] Albert van Dijk and Wayne L. Decker, 1987, Smoothing Vegetation Index Profile: An
Alternative Method for Reducing Radiometric Disturbance in NOAA/AVHRR Data.
Photogrammetric Engineering and Remote Sensing. 1059-1067.
原ほか, 2003. 時系列衛星画像のノイズ除去アルゴリズムの開発と評価

写真測量とリモートセンシング Vol.42, 49-59
澤田ほか, 2000. 高頻度観測衛星による植生の季節変動観測 日本写真測量学会秋季学術講演
論文集, 209-212
朴ほか, 1999. AVHRR NDVI 時系列データの高周波ノイズの除去のための TWO 法の提案 写
真測量とリモートセンシング Vol.38, 36-47[4] MODIS Web http://modis.gsfc.nasa.gov/[5] 農林水産研究計算センターSIDaB
http://rms1.agsearch.agropedia.affrc.go.jp/sidab/index-ja.html
東京情報大学フロンティアプロジェクト
http://www.frontier.tuis.ac.jp/modis/frontier/index.html
東京大学 Web MODIS http://webmodis.iis.u-tokyo.ac.jp/index_j.php
JAXA/地球観測利用推進センター
http://kuroshio.eorc.jaxa.jp/ADEOS/mod_nrt/

[6] NASA Earth Observing System data gateway
http://edcimswww.cr.usgs.gov/pub/imswelcome/[7] 改訂版 図解リモートセンシング 日本リモートセンシング研究会編, 発行:社団法人日本測量協会

pp.325[8] Huete, et al., 2002, Overview of the radiometric and biophysical performance of theMODIS vegetation indices. Remote Sensing of Environment, 83, 195-213[9] Sakamoto, et al., 2005, A crop phenology detection method using time-series MODIS data.
Remote Sensing of Environment, 96, 366-374
坂本ほか, 2004, MODIS 画像時系列解析による水稲の生育モニタリング, システム農学

第20 巻別号 2, 47-48[10] Sakamoto, et al., 2006, Spatio-temporal distribution of rice phenology and cropping
systems in the Mekong Delta with special reference to the seasonal water flow of the Mekong
and Bassac rivers. Remote Sensing of Environment, 100, 1-16
坂本ほか, 2006, メコンデルタにおける水稲作付体系の把握, システム農学
1, 39-40[11] 台湾 国家宇宙計画局(NSPO) http://www.nspo.org.tw/2005e/
株式会社イメージワン (販売代理店)
http://www.imageone.co.jp/enterprise/geo/sat/fs2/index.html

第 22 巻別号