(57)【特許請求の範囲】
【請求項１】
農作物の正常株率を推定する正常株率推定装置であって、
農作物の作付を一意に識別する作付ＩＤと、作付が実施された圃場を一意に識別する圃場ＩＤと、作付が開始された日付である作付開始日と、農作物を収穫し作付を終了する予定の日付である作付終了予定日と、が登録されている作付情報と、
圃場ごとに、圃場の位置と形状とを表す情報が登録されている、圃場別ＧＩＳ情報と、
農作物の作付や生育の段階を表す生育ステージと、画像データ解析対象の生育ステージであるかどうかを表す情報と、が登録されている生育ステージ情報と、
第一の解像度の画像を一意に識別する画像ＩＤと、前記第一の解像度の画像を撮影した日付である画像撮影日と、前記第一の解像度の画像データ本体と、が登録されている第一の解像度の画像情報と、
第二の解像度の画像を一意に識別する画像ＩＤと、前記第二の解像度の画像を撮影した日付である画像撮影日と、前記第二の解像度の画像データ本体と、が登録されている第二の解像度の画像情報と、を格納する記憶装置と関連付けされており、
入力装置から予測対象の作付ＩＤを設定する予測対象作付設定部と、
設定された前記作付ＩＤに基づいて、前記作付情報と前記第一の解像度の画像情報とを参照して取得した、前記第一の解像度の画像データから、画像撮影日における予測対象作付の圃場の植生指数代表値を算出する画像データ分析部と、
前記生育ステージ情報を参照して、前記植生指数代表値の時系列データに前記農作物の生育モデルを適合して推定植生指数代表値を計算する生育モデル適合部と、
適合した前記農作物の生育モデルの変化率及び日数の少なくともいずれか一方から、作付開始日からの経過日数ごとに生育ステージを推定する生育ステージ推定部と、
画像データ解析対象の生育ステージごとに、解析対象とする第二の解像度の画像データを選択する第二の解像度に対する画像データ選択部と、前記第二の解像度の画像における予測対象作付の圃場のピクセルごとに、前記解析対象とする前記第二の解像度の画像データから計算したピクセル別植生指数を結合して多次元ベクトルを算出する第二の解像度に対する画像データ分析部と、
前記多次元ベクトルを当該ベクトルの値のみからクラスタリングにより分類し、得られたクラスタから正常株の部分のピクセルが分類されるクラスタを所定の条件により一つ以上特定し、圃場の全ピクセルに対する前記特定したクラスタに属するピクセルの割合から圃場内の正常株率を推定する正常株率推定部と
を有することを特徴とする正常株率推定装置。
【請求項２】
前記第一の解像度は前記第二の解像度よりも低い解像度であることを特徴とする請求項１に記載の正常株率推定装置。
【請求項３】
前記正常株率推定部は、
前記クラスタの中心点に対する条件から、正常株の部分のピクセルが分類されるクラスタを特定する請求項１又は２に記載の正常株率推定装置。
【請求項４】
前記植生指数代表値および前記ピクセル別植生指数はそれぞれNDVI、EVIまたはGRVIのいずれか、であることを特徴とする請求項１から３までのいずれか１項に記載の正常株率推定装置。
【請求項５】
農作物の収穫量を予測する収穫量予測装置であって、請求項１から４までのいずれか１項に記載の正常株率推定装置と、
少なくとも前記正常株率を説明変数に含む収穫量予測モデルで、予測対象作付の収穫量を予測する収穫量予測部と、を有する収穫量予測装置。
【請求項６】
前記生育ステージ情報に、生育度集計対象の生育ステージであるかどうかを表す情報が登録されており、
さらに、
生育度集計対象の生育ステージにおける予測対象作付の前記推定植生指数代表値を合計することで予測対象作付の推定生育度を計算する、生育度推定部を備え、
前記収穫量予測部は、少なくとも前記推定生育度と前記正常株率を説明変数に含む収穫量予測モデルとに基づいて、予測対象作付の収穫量を予測することを特徴とする請求項５に記載の収穫量予測装置。
【請求項７】
農作物の正常株率を推定する正常株率推定方法であって、
農作物の作付を一意に識別する作付ＩＤと、作付が実施された圃場を一意に識別する圃場ＩＤと、作付が開始された日付である作付開始日と、農作物を収穫し作付を終了する予定の日付である作付終了予定日と、が登録されている作付情報と、
圃場ごとに、圃場の位置と形状とを表す情報が登録されている、圃場別ＧＩＳ情報と、
農作物の作付や生育の段階を表す生育ステージと、画像データ解析対象の生育ステージであるかどうかを表す情報と、が登録されている生育ステージ情報と、
第一の解像度の画像を一意に識別する画像ＩＤと、前記第一の解像度の画像を撮影した日付である画像撮影日と、前記第一の解像度の画像データ本体と、が登録されている第一の解像度の画像情報と、
第二の解像度の画像を一意に識別する画像ＩＤと、前記第二の解像度の画像を撮影した日付である画像撮影日と、前記第二の解像度の画像データ本体と、が登録されている第二の解像度の画像情報と、を格納する記憶装置の各情報を参照し、
入力装置から予測対象の作付ＩＤを設定する予測対象作付設定ステップと、
設定された前記作付ＩＤに基づいて、前記作付情報と前記第一の解像度の画像情報とを参照して取得した、前記第一の解像度の画像データから、画像撮影日における予測対象作付の圃場の植生指数代表値を算出する画像データ分析ステップと、
前記生育ステージ情報を参照して、前記植生指数代表値の時系列データに前記農作物の生育モデルを適合して推定植生指数代表値を計算する生育モデル適合ステップと、
適合した前記農作物の生育モデルの変化率及び日数の少なくともいずれか一方から、作付開始日からの経過日数ごとに生育ステージを推定する生育ステージ推定ステップと、
画像データ解析対象の生育ステージごとに、解析対象とする第二の解像度の画像データを選択する第二の解像度に対する画像データ選択ステップと、前記第二の解像度の画像における予測対象作付の圃場のピクセルごとに、前記解析対象とする前記第二の解像度の画像データから計算したピクセル別植生指数を結合して多次元ベクトルを算出する第二の解像度に対する画像データ分析ステップと、
前記多次元ベクトルを当該ベクトルの値のみからクラスタリングにより分類し、得られたクラスタから正常株の部分のピクセルが分類されるクラスタを所定の条件により一つ以上特定し、圃場の全ピクセルに対する前記特定したクラスタに属するピクセルの割合から圃場内の正常株率を推定する正常株率推定ステップと、
を有することを特徴とする正常株率推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、農作物の収穫量予測（正常株率推定）技術に関する。
【背景技術】
【０００２】
農作物の供給量や価格は、天候などの影響により大きく変動することがある。例えば、需要者は、安定した農作物の調達を行うため、出荷量、価格、出荷日などについて事前に生産者と契約を交わす、契約取引を行っている。
【０００３】
しかしながら、契約生産者全体の生産量が予定した調達量に満たない場合は、需要者は海外から農作物を輸入するなどして、予定した数量に達するように農作物を別途調達しなければならない。従って、需要者は契約生産者の農作物の収穫量を精度よく予測し、別途調達する農産物の数量を適切に決定する必要がある。
【０００４】
ところで、根菜類などのように、収穫部分が地中に生育する農作物の場合は、収穫前に圃場調査を行い、圃場内の欠株率や、一株当りの収穫部分の重量を調査することで、収穫量予測を行っていた。
【０００５】
しかしながら、この方法による収穫量予測は、契約生産者の圃場が多数あり、かつ、広域に分散している場合などにおいて、手間がかかり非効率的である。また、圃場調査は圃場内の一部の区画でのみ行うため、欠株率や株当りの重量に圃場内でバラつきがあると、精度よく予測できないという問題がある。
【０００６】
このような問題を解決するため、人工衛星によるリモートセンシングで観測された衛星画像から農作物の収穫量を予測する方法が考えられている。
【０００７】
例えば、下記非特許文献１では、人工衛星によるリモートセンシングで観測された予測対象地域の衛星画像から、NDVI（Normalized Difference Vegetation Index）などの植生指数を算出し、ある期間で累積したNDVIの値から当該地域における農作物の収穫量を予測している。
【０００８】
また、下記特許文献１（例えば請求項６）に記載の植生生長分析システムでは、低解像度の衛星データ及び少数の高解像度の衛星データを用いて、NDVIなどの植生指数や植生生長曲線を補正する。これにより、植生指数に基づく収穫量予測の精度を向上することができる。
【０００９】
下記非特許文献２では、圃場調査の結果を教師データとして、解像度衛星画像から農作物の病害の程度を推定する分類木を学習し、圃場内の推定された病害発生状況と実際の収穫量との関係を分析している。
【００１０】
特許文献１（たとえば請求項２）に記載の植生生長分析システムでは、倒伏、病虫害などの災害に対して、予め格納された過去の災害データを教師データとして、高次元分析モデルデータの各点群から災害ありデータと正常生長データを分ける閉曲面・線を学習する。これにより、衛星画像から対象植物の生長が正常か否かを判断することができ、その結果は収穫量予測に適用できる。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】特開２０１５−１８８３３３号公報
【非特許文献】
【００１２】
【非特許文献１】C. J. Weissteiner, W. Kuhbauch, Regional Yield Forecasts of Malting Barley (Hordeum vulgare L.) by NOAA‐AVHRR Remote Sensing Data and Ancillary Data, Journal of Agronomy and Crop Science, Vol. 191 Issue 4, pp.308-320 (2005).
【非特許文献２】Jonas Franke, Gunter Menz, Multi-temporal wheat disease detection by multi-spectral remote sensing, Precision Agriculture, Vol. 8, Issue 3, pp. 161-172 (2007).
【発明の概要】
【発明が解決しようとする課題】
【００１３】
しかしながら、非特許文献１や特許文献１に記載の技術では、地上部分の植生の情報のみから農作物の収穫量を予測している。そのため、欠株による収穫量の減少を早期に予測することが難しく、特に収穫部分が地中に生育する農作物では予測誤差が大きくなってしまうという問題がある。
【００１４】
また、非特許文献２や特許文献１に記載の方法では、分類モデルの学習用に教師データが必要となる。これらの文献では、欠株の推定の場合は、正常株もしくは欠株がラベリングされた衛星画像データを教師データとして、衛星画像データが正常株のものか欠株のものかを判別する分類モデルを学習する。
【００１５】
ここで、農作物の地上部分の植生は生育の段階に応じて異なる。そのため、同じ正常株であっても、生育の段階によって観測される衛星画像データは異なる。
【００１６】
したがって、生育の段階ごとに別の教師データを用意して分類モデルを学習する必要があるが、そのための教師データを実際の圃場調査や衛星画像の分析から十分な量だけ得るのは難しいという問題がある。
【００１７】
本発明は、欠株による収穫量の減少を早期に予測し、農作物の収穫量を高い精度で予測することができる技術を提供することを目的とする。
【課題を解決するための手段】
【００１８】
本発明の一観点によれば、農作物の正常株率を推定する正常株率推定装置であって、農作物の作付を一意に識別する作付ＩＤと、作付が実施された圃場を一意に識別する圃場ＩＤと、作付が開始された日付である作付開始日と、農作物を収穫し作付を終了する予定の日付である作付終了予定日と、が登録されている作付情報と、圃場ごとに、圃場の位置と形状とを表す情報が登録されている、圃場別ＧＩＳ情報と、農作物の作付や生育の段階を表す生育ステージと、画像データ解析対象の生育ステージであるかどうかを表す情報と、が登録されている生育ステージ情報と、第一の解像度の画像を一意に識別する画像ＩＤと、前記第一の解像度の画像を撮影した日付である画像撮影日と、前記第一の解像度の画像データ本体と、が登録されている第一の解像度の画像情報と、第二の解像度の画像を一意に識別する画像ＩＤと、前記第二の解像度の画像を撮影した日付である画像撮影日と、前記第二の解像度の画像データ本体と、が登録されている第二の解像度の画像情報と、を格納する記憶装置と関連付けされており、入力装置から予測対象の作付ＩＤを設定する予測対象作付設定部と、設定された前記作付ＩＤに基づいて、前記作付情報と前記第一の解像度の画像情報とを参照して取得した、前記第一の解像度の画像データから、画像撮影日における予測対象作付の圃場の植生指数代表値を算出する画像データ分析部と、前記生育ステージ情報を参照して、前記植生指数代表値の時系列データに前記農作物の生育モデルを適合して推定植生指数代表値を計算する生育モデル適合部と、適合した前記農作物の生育モデルの変化率及び日数の少なくともいずれか一方から、作付開始日からの経過日数ごとに生育ステージを推定する生育ステージ推定部と、画像データ解析対象の生育ステージごとに、解析対象とする第二の解像度の画像データを選択する第二の解像度に対する画像データ選択部と、前記第二の解像度の画像における予測対象作付の圃場のピクセルごとに、前記解析対象とする前記第二の解像度の画像データから計算したピクセル別植生指数を結合して多次元ベクトルを算出する第二の解像度に対する画像データ分析部と、前記多次元ベクトルを当該ベクトルの値のみからクラスタリングにより分類し、得られたクラスタから正常株の部分のピクセルが分類されるクラスタを所定の条件により一つ以上特定し、圃場の全ピクセルに対する前記特定したクラスタに属するピクセルの割合から圃場内の正常株率を推定する正常株率推定部と、を有することを特徴とする正常株率推定装置が提供される。
【００１９】
上記の正常株率推定装置によれば、前記多次元ベクトルを当該ベクトルの値のみからクラスタリングにより分類し、正常株の部分のピクセルが分類されるクラスタを所定の条件により特定するため教師なし学習をすることができる。
【００２０】
また、前記第一の解像度は前記第二の解像度よりも低い解像度であることが好ましい。  前記正常株率推定部は、前記クラスタの中心点に対する条件から、正常株の部分のピクセルが分類されるクラスタを特定するようにすると良い。一例として、「植生停止」の生育ステージにおける植生指数の値が最も大きいｃｋを選択することで、正常株のピクセルのクラスタを特定する方法を用いることができる。
【００２１】
正常な株では、「植生停止」の生育ステージまで植生指数が増加していくのに対して、欠株部分では「植生増加」の生育ステージの途中から植生指数が増加しなくなる。また、株が作付されていない土壌部分は「定植直後」の生育ステージから継続して植生指数が低いままとなる。
【００２２】
前記植生指数代表値および前記ピクセル別植生指数はそれぞれNDVI、EVIまたはGRVIのいずれか、であることを特徴とする。
【００２３】
前記画像データ分析部は、前記第二の解像度の画像データから当該圃場の画像を抽出し、当該圃場の位置や向きが一致するように座標変換を行い、座標変換処理後の当該圃場の画像の各ピクセルＰに対して、ピクセルＰにおける植生指数ＶＩ（Ｓ，Ｐ）を計算することが好ましい。
【００２４】
本発明は、農作物の収穫量を予測する収穫量予測装置であって、上記に記載の正常株率推定装置と、少なくとも前記正常株率を説明変数に含む収穫量予測モデルで、予測対象作付の収穫量を予測する収穫量予測部と、を有する収穫量予測装置であっても良い。ここで、推定率は、欠株率でも良い。
【００２５】
前記生育ステージ情報に、生育度集計対象の生育ステージであるかどうかを表す情報が登録されており、さらに、生育度集計対象の生育ステージにおける予測対象作付の前記推定植生指数代表値を合計することで予測対象作付の推定生育度を計算する、生育度推定部を備え、前記収穫量予測部は、少なくとも前記推定生育度と前記正常株率を説明変数に含む収穫量予測モデルとに基づいて、予測対象作付の収穫量を予測することを特徴とする。ここで、前記推定植生指数代表値は、例えば、植生指数の中央値や最頻値、平均値などである。
【００２６】
本発明の他の観点によれば、農作物の正常株率を推定する正常株率推定方法であって、農作物の作付を一意に識別する作付ＩＤと、作付が実施された圃場を一意に識別する圃場ＩＤと、作付が開始された日付である作付開始日と、農作物を収穫し作付を終了する予定の日付である作付終了予定日と、が登録されている作付情報と、圃場ごとに、圃場の位置と形状とを表す情報が登録されている、圃場別ＧＩＳ情報と、農作物の作付や生育の段階を表す生育ステージと、画像データ解析対象の生育ステージであるかどうかを表す情報と、が登録されている生育ステージ情報と、第一の解像度の画像を一意に識別する画像ＩＤと、前記第一の解像度の画像を撮影した日付である画像撮影日と、前記第一の解像度の画像データ本体と、が登録されている第一の解像度の画像情報と、第二の解像度の画像を一意に識別する画像ＩＤと、前記第二の解像度の画像を撮影した日付である画像撮影日と、前記第二の解像度の画像データ本体と、が登録されている第二の解像度の画像情報と、を格納する記憶装置の各情報を参照し、入力装置から予測対象の作付ＩＤを設定する予測対象作付設定ステップと、設定された前記作付ＩＤに基づいて、前記作付情報と前記第一の解像度の画像情報とを参照して取得した、前記第一の解像度の画像データから、画像撮影日における予測対象作付の圃場の植生指数代表値を算出する画像データ分析ステップと、前記生育ステージ情報を参照して、前記植生指数代表値の時系列データに前記農作物の生育モデルを適合して推定植生指数代表値を計算する生育モデル適合ステップと、適合した前記農作物の生育モデルの変化率及び日数の少なくともいずれか一方から、作付開始日からの経過日数ごとに生育ステージを推定する生育ステージ推定ステップと、画像データ解析対象の生育ステージごとに、解析対象とする第二の解像度の画像データを選択する第二の解像度に対する画像データ選択ステップと、前記第二の解像度の画像における予測対象作付の圃場のピクセルごとに、前記解析対象とする前記第二の解像度の画像データから計算したピクセル別植生指数を結合して多次元ベクトルを算出する第二の解像度に対する画像データ分析ステップと、前記多次元ベクトルを当該ベクトルの値のみからクラスタリングにより分類し、得られたクラスタから正常株の部分のピクセルが分類されるクラスタを所定の条件により一つ以上特定し、圃場の全ピクセルに対する前記特定したクラスタに属するピクセルの割合から圃場内の正常株率を推定する正常株率推定ステップと、を有することを特徴とする正常株率推定方法が提供される。
【００２７】
本発明は、上記に記載の正常株率推定方法をコンピュータに実行させるためのプログラムであっても良く、当該プログラムを記録するコンピュータ読み取り可能な記憶媒体であっても良い。
【発明の効果】
【００２８】
本発明によれば、欠株による収穫量の減少を早期に予測し、農作物の収穫量を高い精度で予測することができる技術を提供することができる。
【図面の簡単な説明】
【００２９】
【図１Ａ】本発明の実施の形態による収穫量予測装置の概略構成例を示す機能ブロック図である。
【図１Ｂ】図１Ａの構成のより詳細な構成例を示す機能ブロック図である。
【図２】作付情報テーブルの一構成例を示す図である。
【図３】圃場別GIS情報テーブルの一構成例を示す図である。
【図４】生育ステージテーブルの一構成例を示す図である。
【図５】低解像度画像情報テーブルの一構成例を示す図である。
【図６】高解像度画像情報テーブルの一構成例を示す図である。
【図７】作付別植生指数テーブルの一構成例を示す図である。
【図８】推定生育度テーブルの一構成例を示す図である。
【図９】ピクセル別分析結果テーブルの一構成例を示す図である。
【図１０】推定正常株率テーブルの一構成例を示す図である。
【図１１】予測収穫量テーブルの一構成例を示す図である。
【図１２】実績収穫量テーブルの一構成例を示す図である。
【図１３】本実施の形態による農作物の収穫量予測処理の流れを示すフローチャート図である。
【図１４】低解像度画像データ分析部における解像度画像データ分析処理の流れを示すフローチャート図である。
【図１５】生育モデル適合部における生育モデル適合処理の流れを示すフローチャート図である。
【図１６】生育度推定部における生育度推定処理の流れを示すフローチャート図である。
【図１７Ａ】高解像度画像データ選択部における処理の流れを示すフローチャート図である。
【図１７Ｂ】図１７Ａに続く図である。
【図１８】高解像度画像データ分析部における高解像度画像データ分析処理の流れを示すフローチャート図である。
【図１９】正常株率推定部における正常株率推定処理の流れを示すフローチャート図である。
【図２０】生育モデル適合部における生育モデル適合処理のイメージを示す図である。
【図２１】生育ステージ推定部における生育ステージ推定処理のイメージを示す図である。
【図２２】生育度推定部における生育度推定処理のイメージを示す図である。
【図２３】高解像度画像データ選択部における高解像度画像データ選択処理のイメージを示す図である。
【図２４】高解像度画像データ選択部の経過日数選択処理のイメージを示す図である。
【図２５】座標変換処理のイメージを示す図である。
【図２６】ピクセルＰにおけるピクセル別植生指数ＶＩ（Ｓ，Ｐ）を計算する処理のイメージ図である。
【図２７Ａ】第１のクラスタ処理例のイメージを示す図である。
【図２７Ｂ】第２のクラスタ処理例のイメージを示す図である。
【発明を実施するための形態】
【００３０】
以下、本発明の一実施の形態による農作物の正常株率推定技術に基づく収穫量予測技術について、農作物として根菜類を例にし、図面を参照しながら詳細に説明を行う。但し、農作物は根菜類に限定されるものではない。
【００３１】
まず、本実施の形態による農作物の収穫量予測装置のシステム構成例について説明する。図１Ａは、本実施の形態による収穫量予測装置Ａの概略構成例を示す機能ブロック図である。図１Ａに示すように、収穫量予測装置Ａは、収穫量予測に利用する情報等を格納するデータベース１ａ、１ｂ、１ｃ、…を有するデータベース群１と、収穫量予測演算を行うための予測演算部Ｂと、これらを関連付けする例えばインターネットなどのネットワークＮＴとを有している。予測演算部Ｂは、各種インタフェース部Ｃと、制御部７と、演算処理を行うためのプログラムなどを記憶するメモリＤとを有している。収穫量予測装置Ａと予測演算部Ｂとは、例えばパーソナルコンピュータＰＣなどからなり、インターネットなどを介さない一体化した構成であっても良い。
【００３２】
図１Ｂは、図１Ａの構成例の、より詳細な構成例を示す機能ブロック図である。図１Ｂに示すように、本実施の形態による収穫量予測装置Ａは、データベース群１に対応する補助記憶装置１と、各種インタフェース部Ｃに対応する、例えば、キーボードなどの入力装置３、ディスプレイなどの出力装置５と、中央演算装置（ＣＰＵ）７を、含み、さらに、各種処理部をソフトウェア構成により実現するためのプログラムを記録する主記憶装置２１を有する。各種処理部は、ハードウェア構成により構成されていても良い。
【００３３】
より詳細には、主記憶装置２１とＣＰＵ７とにより、以下のような各種処理部における処理を実現することができる。
【００３４】
尚、気象情報取得手段１１は、インターネット等を介して、外部のシステムから気象予報を取得する。気象予報には、１か月予報や２か月予報などの長期予報を含む。また、気象情報取得手段１１は、過去の気象観測結果（実績データ）も取得する。その詳細については後述する。
【００３５】
≪各種情報≫
各種情報は、例えば、以下に説明するようなテーブル形式でメモリに格納されている。図２から図１２までは、各種テーブルのデータ構成例を示す図である。
【００３６】
＜作付情報テーブル１１−１＞
図２に示すように、作付情報テーブル１１−１は、少なくとも収穫量の予測対象となる作付の情報が登録されているテーブルである。
【００３７】
例えば、以下のデータを格納する。
農作物の作付を一意に識別する作付ＩＤ、作付が実施された圃場を一意に識別する圃場ＩＤ、作付が開始された日付である作付開始日、農作物を収穫し作付を終了する予定の日付である作付終了予定日、当該作付において作付した株数である作付株数。
【００３８】
＜圃場別GIS情報テーブル１１−２＞
図３に示すように、圃場別GIS情報テーブル１１−２は、作付を実施している圃場の地図上の位置と形状とを表す情報が登録されているテーブルである。
【００３９】
例えば、以下のデータを格納する。
圃場ID、圃場の地図上の位置と形状を表す情報を格納したシェープファイル、当該圃場の面積。ここで、シェープファイルとは、圃場の形状をポリゴンで表現することができるデータ形式であり、座標系に経緯度を使用することで地図上での圃場の位置と形状とを表すことができる。圃場の面積は、単位面積当たりの収穫量を取得したいときに利用することができる。必要に応じて、事前にシェープファイルから面積を求めて登録しておいてもよい。
【００４０】
＜生育ステージテーブル１１−３＞
図４に示すように、生育ステージテーブル１１−３は、農作物の各生育ステージとそれに付随する情報が登録されているテーブルである。
【００４１】
例えば、以下のデータを格納する。
作付や生育の段階を表す、「定植直後」、「植生増加」、「植生停止」、「植生減少」などの生育ステージ、後述する生育度推定部において生育度を推定する際に、推定植生指数平均値の集計対象となる生育ステージかどうかを表す生育度集計対象であるか否かに関する情報、後述する高解像度画像データ選択部において高解像度画像を選択する際に、対象となる生育ステージかどうかを表す高解像度画像データ解析対象であるか否かに関する情報、後述する高解像度画像データ選択部において高解像度画像を選択する際に、作付開始日からの日数である経過日数をもとに、同一の生育ステージ内の複数の高解像度画像から解析対象の画像を選択するための、経過日数選択式等。
【００４２】
尚、高解像度画像データ解析対象のカラムに“No”が登録されているレコードでは、経過日数選択式は登録しなくてよい。
【００４３】
経過日数選択式の例として、以下の式を用いることができる。ここで、Ｓは生育ステージ、DSTART(S) は生育ステージＳの最初の経過日数、DEND(S) は生育ステージSの最後の経過日数、{ Dn(S) }n=1,…,N は高解像度画像が取得可能な生育ステージＳの経過日数であり、DSTART(S) ≦ Dn(S) ≦ DEND(S) を満たすものとする。
【００４４】
【数１】
0006639300.tif 000002
【００４５】
【数２】
0006639300.tif 000003
【００４６】
【数３】
0006639300.tif 000004
【００４７】
ここで、式（１）は、高解像度画像が取得可能な生育ステージSの経過日数から、最も時期が早いものを選択する式である。
【００４８】
式（２）は、高解像度画像が取得可能な生育ステージＳの経過日数から、生育ステージＳ におけるちょうど半分の経過日数である (DEND(S) - DSTART(S)) / 2 に最も近いものを選択する式である。式（３）は、高解像度画像が取得可能な生育ステージSの経過日数から、最も時期が遅いものを選択する式である。
【００４９】
＜低解像度（第１の解像度）画像情報テーブル１１−４＞
図５に示すように、低解像度画像情報テーブル１１−４は、収穫量予測に使用する低解像度の衛星画像の情報が登録されているテーブルである。
【００５０】
例えば、以下のデータを格納する。
低解像度画像を一意に識別する低解像度画像ＩＤ、低解像度画像を撮影した日付である画像撮影日、低解像度画像データ本体のファイルである画像ファイル。
【００５１】
尚、衛星画像には、Landsatなど、農作物、特に地表に現れる葉を検出可能な可視光と近赤外線とのバンドを含むマルチスペクトルの観測衛星センサーで取得されたものであり、かつ、圃場が識別できる程度の解像度で取得されたものを使用する。但し、衛星画像の種類は、実際には、解像度が低いものに限定しなくても良い。
【００５２】
＜高解像度（第２の解像度）画像情報テーブル１１−５＞
図６に示すように、高解像度画像情報テーブル１１−５は、収穫量予測に使用する高解像度の衛星画像の情報が登録されているテーブルである。
【００５３】
例えば、以下のデータを格納する。
高解像度画像を一意に識別する高解像度画像ＩＤ、高解像度画像を撮影した日付である画像撮影日、高解像画像データ本体のファイルである画像ファイル。
【００５４】
衛星画像は、GeoEyeなど、可視光と近赤外線のバンドを含むマルチスペクトルの観測衛星センサーで取得されたもので、かつ圃場内の正常株と欠株の状況が推定できるよう、一つのピクセル内で観測される株数がなるべく少ないものを使用する。ただし、衛星画像の種類は同様に限定しない。
【００５５】
ここで、低解像度と高解像度との区別は、相対的な解像度の違いによるものであり、絶対的な解像度により利用可能な衛星画像を限定するものではない。さらに、低解像度衛星画像には高解像度衛星画像と同一の観測衛星センサーで取得された画像を用いてもよい。
【００５６】
＜作付別植生指数テーブル１１−６＞
図７に示すように、作付別植生指数テーブル１１−６は、作付ごと、作付開始日からの経過日数ごとに、圃場の植生指数平均値の情報を登録するテーブルである。
【００５７】
例えば、以下のデータを格納する。
作付ＩＤ、当該作付の作付開始日からの経過日数、植生指数平均値、後述する低解像度画像データ分析部で計算された、当該経過日数における当該圃場の推定植生指数平均値、後述する生育ステージ推定部により推定された、当該作付の当該経過日数における生育ステージ、後述する高解像度画像データ選択部により登録された、当該経過日数における当該圃場を含む高解像度画像の高解像度画像ＩＤ、高解像度画像データ選択部により登録された、当該高解像度画像を後述する高解像度画像データ分析部で使用するかどうかを表す、高解像度画像選択フラグ。
【００５８】
当該経過日数における当該圃場を含む低解像度画像が取得できない場合は、当該レコードの植生指数平均値は登録されない。また、当該経過日数における当該圃場を含む高解像度画像が取得できない場合は、当該レコードの高解像度画像ＩＤは登録されない。さらにその場合、高解像度画像選択フラグも登録しない。
【００５９】
＜推定生育度テーブル１１−７＞
図８に示すように、推定生育度テーブル１１−７は、収穫部分の推定生育度の情報が登録されるテーブルである。
【００６０】
例えば、以下のデータを格納する。
作付ＩＤ、後述する生育度推定部により推定された、当該作付の推定生育度。
【００６１】
＜ピクセル別分析結果テーブル１１−８＞
図９に示すピクセル別分析結果テーブル１１−８は、作付ごと、ピクセルごと、高解像度画像データ解析対象の生育ステージごとに、ピクセル別植生指数の情報が登録されるテーブルである。
【００６２】
例えば、以下のデータを格納する。
作付ＩＤ、当該圃場の高解像度画像におけるピクセルを一意に識別するピクセルＩＤ、高解像度画像データ解析対象の生育ステージごとに、高解像度画像データ分析部により分析されたピクセル別植生指数。尚、高解像度画像データ解析対象の生育ステージが複数ある場合は、その数だけ植生指数を登録する列が用意される。
【００６３】
＜推定正常株率テーブル１１−９＞
図１０に示すように、推定正常株率テーブル１１−９は、作付ごとに、圃場内の正常株率の推定値の情報が登録されるテーブルである。
【００６４】
例えば、以下のデータを格納する。
作付ＩＤ、正常株率推定部により推定された推定正常株率。
【００６５】
＜予測収穫量テーブル１１−１０＞
図１１に示すように、予測収穫量テーブル１１−１０は、予測対象となる作付ごとに、予測収穫量の情報が登録されるテーブルである。
【００６６】
例えば、以下のデータを格納する。
作付ＩＤ、圃場ＩＤ、作付開始日、当該作付の収穫が予定されている日付である収穫予定日、後述する収穫量予測部により予測された予測収穫量。
【００６７】
＜実績収穫量テーブル１１−１１＞
図１２に示すように、実績収穫量テーブル１１−１１は、過去の作付における実績収穫量の情報が登録されているテーブルである。
【００６８】
例えば、以下のデータを格納する。
作付ＩＤ、圃場ＩＤ、作付開始日、当該作付の収穫が実施された日付である収穫実施日、当該作付の収穫量の実績値である実績収穫量。
【００６９】
次に、全体処理部と各処理部の処理等について説明する。
≪処理の説明≫
＜全体処理＞
図１３は、本実施の形態による農作物の収穫量予測処理の流れを示すフローチャート図である。図１Ｂ等も参照しながら説明を行う。
【００７０】
（ステップＳ１）  予測対象作付設定部２１−１は、入力装置３を介して操作された予測対象作付の作付ＩＤを設定する。
【００７１】
（ステップＳ２）  低解像度画像データ分析部２１−２は、低解像度画像情報テーブル１１−４に登録された低解像度画像データから、当該画像データの画像撮影日における予測対象作付の圃場の植生指数平均値を算出し、作付別植生指数テーブル１１−６に登録する。尚、平均値を用いるのは一例であり、その他の統計処理に基づく代表値を用いても良い。以下の処理においても同様である。
【００７２】
（ステップＳ３）  生育モデル適合部２１−３は、植生指数平均値の時系列データに、植生指数の日々の変化を表す曲線である生育モデルを適合し、適合した値である推定植生指数平均値を計算して、作付別植生指数テーブル１１−６に登録する。
【００７３】
（ステップＳ４）  生育ステージ推定部２１−４は、適合した生育モデルの変化率及び／または日数から、作付開始日からの経過日数ごとに、予測対象作付の生育ステージを推定し、作付別植生指数テーブル１１−６に登録する。
【００７４】
（ステップＳ５）  生育度推定部２１−５は、生育ステージテーブル１１−３から、生育度集計対象の生育ステージを取得し、作付別植生指数テーブル１１−６から予測対象作付の推定植生指数平均値と生育ステージとを取得し、生育度集計対象の生育ステージに該当する予測対象作付の推定植生指数平均値を合計することで、予測対象作付の推定生育度を計算し、推定生育度テーブル１１−７に登録する。
【００７５】
（ステップＳ６）  高解像度画像データ選択部２１−６は、高解像度画像情報テーブル１１−５に登録された高解像度画像から、予測対象作付の圃場を含む高解像度画像を抽出し、生育ステージテーブル１１−３に予め登録された経過日数選択式に基づき、予測対象作付ごと、高解像度画像データ解析対象の生育ステージごとに、解析対象とする高解像度画像データを選択し、作付別植生指数テーブル１１−６に登録する。
【００７６】
（ステップＳ７）  高解像度画像データ分析部２１−７は、解析対象とする高解像度画像データから当該圃場の画像を抽出し、当該圃場の位置が一致するように座標変換を実行し、当該圃場のピクセルごとに、解析対象とする高解像度画像データから計算したピクセル別植生指数を結合して多次元ベクトルを算出し、ピクセル別分析結果テーブル１１−８に登録する。
【００７７】
（ステップＳ８）  正常株率推定部２１−８は、多次元ベクトルをクラスタリングによる教師なし学習で分類することで正常株の部分のピクセルを抽出し、圃場の全ピクセルに対する抽出したピクセルの割合から圃場内の正常株率を推定し、推定正常株率テーブル１１−９に登録する。
ここで、正常株率に加えて或いは正常株率に代えて、欠株率を求めることもできる。
【００７８】
（ステップＳ９）  収穫量予測部２１−９は、少なくとも推定生育度と正常株率を説明変数に含む収穫量予測モデルで予測対象作付の収穫量を予測し、予測収穫量を予測収穫量テーブル１１−１０に登録する。
【００７９】
（ステップＳ１０）  予測収穫量出力部２１−１０は、予測収穫量テーブル１１−１０に登録された予測収穫量を取得し、出力装置５に表示させる。
【００８０】
次に、各処理部における処理例について詳細に説明する。
＜低解像度画像データ分析部２１−２における処理（図１４）＞
ステップＳ２は、詳細には、以下の処理を含む。
【００８１】
（ステップＳ２−１）  低解像度画像データ分析部２１−２は、作付情報テーブル１１−１から、予測対象作付の作付ＩＤをキーとして、予測対象作付の圃場ＩＤと作付開始日ＤＳＴＡＲＴを取得する。
【００８２】
（ステップＳ２−２）圃場別ＧＩＳ情報テーブル１１−２から、予測対象作付の圃場ＩＤをキーとして、当該圃場のＧＩＳ情報を格納したシェープファイルを取得する。ここで、前述の通り、シェープファイルとは、圃場の形状をポリゴンで表現することができるデータ形式であり、座標系に経緯度を使用することで地図上での圃場の位置を表すことができる。シェープファイルを用いることで、衛星画像に当該圃場の画像が含まれているか否かを判定したり、衛星画像から当該圃場の画像を抽出したりすることができる。
【００８３】
（ステップＳ２−３）  低解像度画像情報テーブル１１−４に登録されている全ての低解像度画像について低解像度画像データ分析処理が完了しているか否かを判定する。完了していない場合、ステップＳ２−４に進む。完了している場合、低解像度画像データ分析部２１−２の処理を終了する。
【００８４】
（ステップＳ２−４）  低解像度画像情報テーブル１１−４に登録されている次の低解像度画像データの低解像度画像ＩＤを取得する。
【００８５】
（ステップＳ２−５）  低解像度画像ＩＤをキーとして、低解像度画像情報テーブル１１−４から、当該画像の画像撮影日DIMAGEを取得し、Ｄ＝ＤＩＭＡＧＥ −ＤＳＴＡＲＴ＋１により作付開始日からの経過日数Ｄを算出する。
【００８６】
（ステップＳ２−６）  低解像度画像ＩＤをキーとして、低解像度画像情報テーブル１１−４から低解像度画像ファイルを取得し、シェープファイルを用いて、当該低解像度画像ファイルに当該圃場の画像が含まれているか否かを判定する。
作付開始日からの経過日数Ｄ≧１、かつ、当該低解像度画像ファイルに当該圃場の画像が含まれている場合（ｙｅｓ）、ステップＳ２−７に進む。そうでない場合（ｎｏ）、ステップＳ２−３に戻る。
【００８７】
（ステップＳ２−７）  当該低解像度画像ファイルから、前記シェープファイルを用いて、当該圃場の画像を抽出する。
【００８８】
（ステップＳ２−８）  前記当該圃場の画像に対して、ピクセルごとに植生指数を算出して平均し、植生指数平均値ＶＩを計算する。
【００８９】
植生指数としては、一例として、ＮＤＶＩ（Normalized Difference Vegetation Index）を使用することができる。ピクセルＰにおける近赤外線と可視光赤の反射率を、それぞれ ρNIR(P)およびρRED(P)としたとき、ピクセルＰにおけるＮＤＶＩは、以下の数式により算出する。
【００９０】
【数４】
0006639300.tif 000005
【００９１】
さらに、当該圃場の画像がピクセルP1, … , PN から構成される場合、植生指数平均値VIは以下の数式で算出する。
【００９２】
【数５】
0006639300.tif 000006
【００９３】
その他の植生指数の例として、以下の数式で定義されるEVI （Enhanced Vegetation Index）およびGRVI（Green-Red Vegetation Index）がある。ここで、ρGREEN(P) と ρBLUE(P) はそれぞれピクセルPにおける可視光緑と可視光青の反射率、C1、C2、Lはそれぞれ定数である。
【００９４】
【数６】
0006639300.tif 000007
【００９５】
【数７】
0006639300.tif 000008
【００９６】
EVIは大気や背景土壌の影響を補正した指標である。NDVIに比べて密度の高い植生に対する感度が高く、飽和しにくいとされる。GRVIは可視光緑と可視光赤から計算される指標で、植生の緑色の濃さに反応しやすい。
【００９７】
（ステップＳ２−９）  作付別植生指数テーブル１１−６に、予測対象作付の作付ＩＤと経過日数Ｄを付加して、植生指数平均値ＶＩを登録する。
【００９８】
＜生育モデル適合部２１−３における処理（図１５）、（図２０）＞
ステップＳ３は、詳細には、以下の処理を含む。
【００９９】
（ステップＳ３−１）  生育モデル適合部２１−３は、入力装置３を介して、植生指数の日々の変化を表す曲線である生育モデルｆ（Ｘ）を設定する。低解像度画像を入手できない日の圃場全体での推定植生指数平均値を計算するため、生育モデルを使用する。生育モデル ｆ（Ｘ）の経過日数Ｄにおける値ｆ（Ｄ）を、予測対象作付の経過日数Ｄにおける圃場全体での推定植生指数平均値とする。生育モデルには、例として、以下の式（８）のような多項式関数や、式（９）のガウス関数を用いることができる。ここで、ａ，ｂ，ｃは、生育モデルｆ（ｘ）のパラメータである。
【０１００】
【数８】
0006639300.tif 000009
【０１０１】
【数９】
0006639300.tif 000010
【０１０２】
（ステップＳ３−２）  入力装置３を介して、外れ値の閾値θを設定する。ここで、θは0以上の定数である。尚、外れ値とは、統計において他の値から大きく外れた値として定義できる。
【０１０３】
（ステップＳ３−３）  作付別植生指数テーブル１１−６から、予測対象作付の作付ＩＤをキーとして、予測対象作付に対する経過日数Ｄｎと植生指数平均値ＶＩｎとの組の時系列データ Ｔ＝{（Ｄｎ，ＶＩｎ）}ｎ＝１，２，…，Ｎを取得する。
ここで、Nは予測対象作付の作付ＩＤを含む作付別植生指数テーブル１１−６のレコードで、かつ植生指数平均値が登録されているレコードのレコード数である。
【０１０４】
（ステップＳ３−４）  生育モデルｆ（Ｘ）と時系列データＴとの誤差の二乗和が最小になるように生育モデルｆ（Ｘ）のパラメータを学習する。
【０１０５】
【数１０】
0006639300.tif 000011
【０１０６】
（ステップＳ３−５）  ｆ（Ｄｎ）−ＶＩｎ≧θとなる１≦ｎ≦Ｎが存在するか否かを判定する。存在する場合（ｙｅｓ）、ステップＳ３−６に進む。存在しない場合（ｎｏ）、ステップＳ３−７に進む。
【０１０７】
（ステップＳ３−６）  時系列データＴから点（Ｄｎ，ＶＩｎ）を除外し、Ｎから１を減算して、ステップＳ３−４に戻る。
【０１０８】
一般的に、衛星画像では雲など大気の影響や、土壌水分、太陽光の入射・反射角の影響により、本来観測されるべき植生指数平均値よりも小さな値が植生指数平均値ＶＩｎとして観測されることがある。そのようなＶＩｎを使用して生育モデルｆ（Ｘ）のパラメータを学習すると、推定植生指数平均値が過小評価されるおそれがある。例えば、図２０上図に示すように、外れ値があると、推定植生指数平均値が過小評価される。そこで、図２０下図に示すように、外れ値を除外すると、本来観測されるべき植生指数に近い値になる。
【０１０９】
そこで、推定植生指数平均値ｆ（Ｄｎ）より植生指数平均値ＶＩｎが閾値θ以上下回るときは、点（Ｄｎ，ＶＩｎ）を外れ値として学習用データから除外し、再度パラメータを学習する。これにより、生育モデルｆ（Ｘ）による推定植生指数平均値が、雲などの影響を受けない本来観測されるべき植生指数平均値により近い値になる。
【０１１０】
（ステップＳ３−７）  作付情報テーブル１１−１から、予測対象作付の作付ＩＤをキーとして、予測対象作付の作付開始日ＤＳＴＡＲＴと作付終了予定日ＤＥＮＤを取得する。また、経過日数Ｄを１で初期化する。
【０１１１】
（ステップＳ３−８）  生育モデルｆ（Ｘ）の経過日数Ｄにおける値 ｆ（Ｄ）を予測対象作付の経過日数Ｄにおける推定植生指数平均値として算出し、作付別植生指数テーブル１１−６に、予測対象作付の作付ＩＤと経過日数Ｄとを付加して、推定植生指数平均値ｆ（Ｄ）の値を登録する。
【０１１２】
（ステップＳ３−９）  Ｄに１を加算する。
【０１１３】
（ステップＳ３−１０）  Ｄ＞ＤＥＮＤ−ＤＳＴＡＲＴ＋１か否かを判定する。そうであれば（ｙｅｓ）、生育モデル適合部２１−３の処理を終了する。そうでなければ（ｎｏ）、ステップＳ３−８に戻る。
【０１１４】
＜生育ステージ推定部２１−４における処理、図２１＞
ステップＳ４に示すように、生育ステージ推定部２１−４では、生育モデル適合部２１−３により適合した生育モデルの変化率及び／または日数から、図２１に示すように、作付開始日からの経過日数ごとに、作付や生育の段階を表す「定植直後」、「植生増加」、「植生停止」、「植生減少」などの生育ステージを推定する。生育ステージの分割には、例えば、以下のような基準を用いることができる。
【０１１５】
生育モデルｆ（Ｘ）の変化率を表す微分ｄｆ（Ｘ）／ｄＸの値が正から負に変わり、ｆ（Ｘ）の値が増加から減少に転じる日（以下「生育モデル極大の日」と称する）ＤＭＡＸを求めることができる。ＤＭＡＸは、作付開始（図２１の原点（横軸））から数えた日数である。以下、作付開始日を１日目とし、生育モデル極大の日をＤＭＡＸ日目、作付終了予定日をＤＥＮＤ日目とする。生育ステージ分割の一例として、以下のように定義する方法が考えられる。
【０１１６】
生育ステージは、図２１に示すように、以下のように１）〜４）までのステージが設定される。
１）「定植直後」は、１日目〜１０日目、
２）「植生増加」は、１１日目〜(ＤＭＡＸ−３１)日目、
３）「植生停止」は、(ＤＭＡＸ−３０)日目〜(ＤＭＡＸ−１)日目、
４）「植生減少」は、ＤＭＡＸ日目〜ＤＥＮＤ日目。
【０１１７】
生育ステージテーブル１１−３において生育度集計対象とされる生育ステージには、収穫部分が成長する時期の生育ステージを指定することが望ましい。詳細は、以下の＜生育度推定部＞の説明欄に記載する。
【０１１８】
また、生育ステージテーブル１１−３において高解像度画像データ解析対象とされる生育ステージには、当該生育ステージにおいて高解像度画像の撮影日が最低１回は存在している必要がある。さらに、雲など天候の影響を考慮すると、同一生育ステージ内に高解像度画像の撮影日が複数回含まれていることが望ましい。そのため、衛星の撮影周期に応じて、最低１０日間程度の十分な期間を設けるとよい。
【０１１９】
＜生育度推定部２１−５における処理、図１６、図２２＞
ステップＳ５は、詳細には、以下の処理を含む。
【０１２０】
生育度推定部２１−５では、収穫量予測に用いるため、予測対象作付における収穫時点での収穫部分の生育度を推定する。
【０１２１】
事前の準備として、生育ステージテーブル１１−３には、収穫部分が成長する時期の生育ステージが、生育度集計対象の生育ステージとして登録されている。この時期の地上部分の植生は、収穫部分が成長するための光合成活動に用いられる。そのため、この時期の植生指数平均値を累計することで、収穫時点での収穫部分の成長度合いを推定することができる。そこで、この合計値を予測対象作付の推定生育度として用いる。
【０１２２】
ステップＳ５−１において、生育度推定部２１−５は、生育ステージテーブル１１−３から、生育度集計対象のカラムが“Ｙｅｓ”となっているレコード抽出することで、生育度集計対象の生育ステージＳ１，…，ＳＮを取得する。ここで、Ｎは抽出したレコード数である。一例として、図２２に示すように、地上の植生部分がほぼ最大に達して生育が停止し、光合成活動が最も盛んになる「植生停止」から、収穫時期を含む「植生減少」までの生育ステージを生育度集計対象とする。
【０１２３】
ステップＳ５−２において、作付情報テーブル１１−１から、予測対象作付の作付ＩＤをキーとして、予測対象作付の作付開始日ＤＳＴＡＲＴと作付終了予定日ＤＥＮＤを取得する。また、経過日数Ｄを１で初期化する。
【０１２４】
ステップＳ５−３で、推定生育度Ｖを０で初期化する。
【０１２５】
ステップＳ５−４で、作付別植生指数テーブル１１−６から、予測対象作付の作付ＩＤと経過日数Ｄをキーとして、予測対象作付の経過日数Ｄにおける推定植生指数平均値ｆ（Ｄ）と、生育ステージＳ（Ｄ）とを取得する。
【０１２６】
ステップＳ５−５において、
【数１１】
0006639300.tif 000012
か否か、すなわち生育ステージＳ（Ｄ）が生育度集計対象の生育ステージＳ１，…，ＳＮのいずれかであるか否かを判定する。そうであれば（ｙｅｓ）、ステップＳ５−６に進む。そうでなければ（ｎｏ）、ステップＳ５−７に進む。
【０１２７】
ステップＳ５−６において、推定生育度Ｖに ｆ（Ｄ）を加算する。
【０１２８】
ステップＳ５−７において、Ｄに１を加算する。
【０１２９】
ステップＳ５−８において、Ｄ＞ＤＥＮＤ-ＤＳＴＡＲＴ+１か否かを判定する（集計終了日に到達）。そうであれば（ｙｅｓ）、ステップＳ５−９に進む。そうでなければ（ｎｏ）、ステップＳ５−４に戻る。
【０１３０】
ステップＳ５−９において、推定生育度テーブル１１−７に、予測対象作付の作付ＩＤを付加して、推定生育度Ｖを登録する。
【０１３１】
＜高解像度画像データ選択部２１−６における処理、図１７Ａ、図１７Ｂ、図２３，図２４＞
ステップＳ６は、詳細には、以下の処理を含む。
【０１３２】
ステップＳ６−１において、高解像度画像データ選択部２１−６は、作付情報テーブル１１−１から、予測対象作付の作付ＩＤをキーとして、予測対象作付の圃場ＩＤと作付開始日ＤＳＴＡＲＴを取得する。
【０１３３】
ステップＳ６−２において、圃場別ＧＩＳ情報テーブル１１−２から、予測対象作付の圃場ＩＤをキーとして、当該圃場のＧＩＳ情報を格納したシェープファイルを取得する。シェープファイルについては、前述の通りである。
【０１３４】
ステップＳ６−３において、高解像度画像情報テーブル１１−５に登録されている全ての高解像度画像について高解像度画像データ選択処理が完了しているか否かを判定する。完了していない場合（ｎｏ）、ステップＳ６−４に進む。完了している場合（ｙｅｓ）、ステップＳ６−８に進む。
【０１３５】
ステップＳ６−４において、高解像度画像情報テーブル１１−５に登録されている次の高解像度画像データの高解像度画像ＩＤを取得する。
【０１３６】
ステップＳ６−５において、この高低解像度画像ＩＤをキーとして、高解像度画像情報テーブル１１−５から、当該画像の画像撮影日ＤＩＭＡＧＥを取得し、Ｄ＝ＤＩＭＡＧＥ−ＤＳＴＡＲＴ ＋１により作付開始日からの経過日数Ｄを算出する。
【０１３７】
ステップＳ６−６において、ステップＳ６−４で取得した高解像度画像ＩＤをキーとして、高解像度画像情報テーブル１１−５から高解像度画像ファイルを取得し、ステップＳ６−２で収録したシェープファイルを用いて、当該高解像度画像ファイルに当該圃場の画像が含まれているか否かを判定する。Ｄ≧１、かつ、当該高解像度像ファイルに当該圃場の画像が含まれている場合（ｙｅｓ）、ステップＳ６−７に進む。そうでない場合（ｎｏ）、ステップＳ６−３に戻る。
【０１３８】
ステップＳ６−７において、作付別植生指数テーブル１１−６に、予測対象作付の作付ＩＤと経過日数Ｄとを付加して、ステップＳ６−４において取得した高解像度画像ＩＤを登録する。その後、ステップＳ６−３に戻る。
【０１３９】
ステップＳ６−８において、生育ステージテーブル１１−３に登録されている全ての生育ステージに対して、以下の処理が完了しているか否かを判定する。完了していない場合（ｎｏ）、ステップＳ６−９に進む。完了している場合（ｙｅｓ）、高解像度画像データ選択部２１−６の処理を終了する。
【０１４０】
ステップＳ６−９において、生育ステージテーブル１１−３から、次の生育ステージＳを取得する。
【０１４１】
ステップＳ６−１０において、当該レコードの高解像度画像データ解析対象のカラムが“ｙｅｓ”であれば、生育ステージＳは高解像度画像データ解析対象であると判定し、ステップステップＳ６−１１に進む。“ｎｏ”であれば、ステップＳ６−８に戻る。
【０１４２】
ステップＳ６−１１において、生育ステージテーブル１１−３から、生育ステージＳをキーとして、経過日数選択式Ｑを取得する。
【０１４３】
図４に記載の経過日数選択式の例に従えば、図２３、図２４に示すように、高解像度画像データ選択処理は、以下のようになる。
ａ）「定植直後」は経過日数の早いものが選択されるため、高解像度画像撮影日１）となる（Ｓ１）。
ｂ）「植生増加」は当該生育ステージにおけるちょうど半分の経過日数（図２４の点線の箇所）に最も近いものが選択されるため、高解像度画像撮影日４）となる（Ｓ２）。
ｃ）「植生停止」は、経過日数の最も遅いものが選択されるため、高解像度画像撮影日９）となる（Ｓ３）。
尚、経過日数選択式の詳細は＜生育ステージテーブル＞の欄で説明する。
【０１４４】
ステップＳ６−１２において、作付別植生指数テーブル１１−６から、予測対象作付の作付ＩＤと生育ステージＳをキーとしてレコードを抽出し、抽出したレコードの経過日数から、最初の経過日数ＤＳＴＡＲＴ(Ｓ)、最後の経過日数 ＤＥＮＤ(Ｓ)、高解像度画像ＩＤが登録されている経過日数{Ｄｎ（Ｓ）}ｎ＝１，…，Ｎを取得する。ここで、Ｎは抽出したレコードのうち高解像度画像ＩＤが登録されているレコードのレコード数である。尚、図には記載していないが、Ｎ＝０のときは以下の処理を省略し、ステップＳ６−８に戻る。
【０１４５】
ステップＳ６−１３において、経過日数の添え字ｎを、ｎ＝Ｑ（ＤＳＴＡＲＴ（Ｓ），ＤＥＮＤ（Ｓ）， {Ｄｎ（Ｓ）}ｎ＝１，…，Ｎ）により計算する。
【０１４６】
ステップＳ６−１４において、作付別植生指数テーブル１１−６に、予測対象作付の作付ＩＤと経過日数Ｄｎ（Ｓ）を付加して、高解像度画像選択フラグのカラムに“Ｙｅｓ”を登録する。ここでｎは、ステップＳ６−１３において計算された経過日数の添え字ｎである。
【０１４７】
＜高解像度画像データ分析部２１−７における処理、図１８、図２５、図２６＞
ステップＳ７は、詳細には、以下の処理を含む。
【０１４８】
ステップＳ７−１において、高解像度画像データ分析部２１−７は、作付情報テーブル１１−１から、予測対象作付の作付ＩＤをキーとして、予測対象作付の圃場ＩＤを取得する。
【０１４９】
ステップＳ７−２において、圃場別ＧＩＳ情報テーブル１１−２から、予測対象作付の圃場ＩＤをキーとして、当該圃場のＧＩＳ情報を格納したシェープファイルを取得する。シェープファイルについては、前述の通りである。
【０１５０】
ステップＳ７−３において、生育ステージテーブル１１−３に登録されている全ての生育ステージに対して、以下の処理が完了しているか否かを判定する。完了していない場合（ｎｏ）、ステップＳ７−４に進む。完了している場合（ｙｅｓ）、高解像度画像データ分析部２１−７の処理を終了する。
【０１５１】
ステップＳ７−４においては、生育ステージテーブル１１−３から、次の生育ステージＳを取得する。
【０１５２】
ステップＳ７−５において、生育ステージテーブル１１−３から、生育ステージＳをキーとしてレコードを取得する。当該レコードの高解像度画像データ解析対象のカラムが“ｙｅｓ”であれば、生育ステージＳは高解像度画像データ解析対象であると判定し、ステップＳ７−６に進む。“ｎｏ”であれば、ステップＳ７−３に戻る。
【０１５３】
ステップＳ７−６において、作付別植生指数テーブル１１−６から、予測対象作付の作付ＩＤと生育ステージＳをキーとしてレコードを抽出し、さらに高解像度画像選択フラグのカラムが“Ｙｅｓ”であるレコードを抽出して、当該レコードの高解像度画像ＩＤを取得する。なお、高解像度画像データ選択部２１−６の処理より、当該レコードは存在すれば一意に定まる。存在しない場合は、予測対象作付の収穫量予測はできないものとし、これ以降の全ての処理を中断する。
【０１５４】
ステップＳ７−７において、当該高解像度画像データから当該圃場の画像を抽出し、当該圃場の位置や向きが一致するように座標変換を実行する。図２５は、座標変換処理の概要を示す図である。高解像度画像データにおける当該圃場の位置や向きは、衛星の撮影時の位置や角度により異なる。そこで、当該圃場の位置や向きが一致するよう、抽出した当該圃場の画像を平行移動、回転、拡大縮小により座標変換し、例えば、高解像度画像１と高解像度画像２とで、当該圃場の最北部と最西部が同じ座標になるようにする。座標変換の結果ピクセル値の再計算が必要になる場合は、最近傍法などの手法で再サンプリングを行う。
【０１５５】
次いで、ステップＳ７−８に示すように、ステップＳ７−７の座標変換処理後の当該圃場の画像の各ピクセルＰに対して、ピクセルＰにおけるピクセル別植生指数ＶＩ（Ｓ，Ｐ）を計算する（図２６参照）。そして、予測対象作付の作付ＩＤとピクセルＰのＩＤを付加し、ピクセル別分析結果テーブル１１−８における生育ステージＳのピクセル別植生指数のカラムにＶＩ（Ｓ，Ｐ）を登録する。
【０１５６】
ここで、ピクセル別植生指数には、低解像度画像データ分析部２１−２で用いる植生指数とは別の指標を用いてもよい。特に（７）式で定義されるGRVIは、葉の緑色の濃さに反応しやすい指標であり、欠株部分の識別に好適である。また、一般的に可視光の衛星画像は近赤外線よりも解像度の高いパンシャープン画像が取得できる場合が多い。GRVIは可視光の反射率のみから算出される指標であるので、近赤外線を使用したNDVIよりも高い解像度で分析できるという利点がある。
【０１５７】
＜正常株率推定部２１−８における処理、図１９、図２７Ａ、図２７Ｂ＞
ステップＳ８は、詳細には、以下の処理を含む。
【０１５８】
ステップＳ８−１において、正常株率推定部２１−８は、生育ステージテーブル１１−３から、高解像度画像データ解析対象のカラムが“Ｙｅｓ”となっているレコードを抽出することで、高解像度画像データ解析対象の生育ステージＳ１，…，ＳＮを取得する。ここで、Ｎは抽出したレコード数であり、高解像度画像データ解析対象の生育ステージの数である。
【０１５９】
ステップＳ８−２において、ピクセル別分析結果テーブル１１−８から、予測対象作付の作付ＩＤをキーとしてレコードを抽出し、ピクセルＩＤのカラムの情報から、抽出したレコードに登録されているピクセル {Ｐｍ}ｍ＝１，…，Ｍを取得する。ここで、Ｍは抽出したレコード数であり、高解像度画像から抽出した予測対象作付の圃場の画像全体のピクセル数である。
【０１６０】
ステップＳ８−３において、ピクセル別分析結果テーブル１１−８から、予測対象作付の作付ＩＤをキーとしてレコードを抽出し、各ピクセル別植生指数のカラムの情報から、以下の数式で定義されるＮ次元のピクセル別植生指数のベクトルの集合Ｘを取得する。
【０１６１】
【数１２】
0006639300.tif 000013
【０１６２】
ステップＳ８−４において、前記Ｎ次元のピクセル別植生指数のベクトルの集合Ｘをクラスタリング手法によりＫ個のクラスタＣ１，…，Ｃｋに分割する。ここで、Ｋはクラスタの分割数であり、少なくとも2以上の整数値として事前に登録されているものとする。
【０１６３】
クラスタリング手法には、Fuzzy c-means 法などのファジークラスタリング手法を用いてもよい。Fuzzy c-means 法の場合の出力は、各クラスタの中心点ｃｋ(k = 1, … , K) と、各ベクトル x ∈ X の各クラスタＣｋへの帰属度 ｗｋ(x) である。ここで、帰属度 ｗｋ(x) は0以上1以下の実数である。
【０１６４】
図２７Ａに示すように、ステップＳ８−５において、各クラスタＣｋ（ｋ＝１，…，Ｋ）に対して、クラスタＣｋ に属するベクトルの平均を計算して、クラスタの中心点ｃｋを算出する。クラスタＣｋ の要素数をＭｋ とすると、クラスタの中心点ｃｋは、以下の数式で算出される。
【０１６５】
【数１３】
0006639300.tif 000014
【０１６６】
なお、Fuzzy c-means 法では出力として既に得られているので、算出する必要はない。
【０１６７】
ステップＳ８−６において、クラスタの中心点ｃ１，…，ｃｋに対する条件から、正常株のピクセルのクラスタＣｋを特定する。
【０１６８】
一例として、「植生停止」の生育ステージにおける植生指数の値が最も大きいｃｋを選択することで、正常株のピクセルのクラスタＣｋを特定する方法を用いることができる。
【０１６９】
正常な株では、「植生停止」の生育ステージまで植生指数が増加していくのに対して、欠株部分では「植生増加」の生育ステージの途中から植生指数が増加しなくなる。また、株が作付されていない土壌部分は「定植直後」の生育ステージから継続して植生指数が低いままとなる。
【０１７０】
このことから、「植生停止」の生育ステージにおける植生指数の値が最も大きいｃｋは、正常株のピクセルのクラスタの中心点であると考えられる。「植生停止」の生育ステージがＳＮであり、cj = (c1j, … , cNj) (j = 1, … , K) であったとすると、上記の方法を用いる場合に特定されるクラスタの添え字ｋは、以下の数式で算出できる。
【０１７１】
【数１４】
0006639300.tif 000015
【０１７２】
このように、本実施の形態においては、図２７Ａに示すように、正常株の部分のピクセルが分類されるクラスタを所定の条件により特定することが可能であり、従って、教師なし学習を使用しても正常株を精度良く推定することができるため、分類モデルを学習するために教師データを用意する必要がないという利点がある。
【０１７３】
図２７Ｂは、単一の生育ステージの情報だけでは欠株部分と正常株のピクセルを識別できない場合の例を示す図である。
【０１７４】
入力装置を介して設定された閾値ηに対して、下記の条件を満たすckを選択することで、正常株のピクセルのクラスタＣｋを特定する方法がある。ここでηは正の定数であり、SL (1 ≦L < N) は「植生増加」の生育ステージであるとする。
【０１７５】
全ての j = 1, … , K に対して、次の（A）または（B）が成り立つ：
（A）cNk− η ≧cNj
（B）|cNk− cNj| ≦ η かつ cLk ≧cLj 
この例は、「植生停止」の時期に欠株部分が周囲の株の葉で覆われ、欠株部分のピクセルと正常株のピクセルとで植生指数に差が生じない場合などに適した例である。このような場合、先の例では「植生停止」の生育ステージにおける植生指数のみ使用するため、誤ったクラスタが選択されてしまう可能性がある。そこで、「植生停止」の生育ステージにおける植生指数の値にη以上の差がないときは、「植生増加」の生育ステージにおける植生指数の値で比較し、その値が最も大きくなるckを選択する。
【０１７６】
図２７Ｂに示す例のように、多次元ベクトルを用いて複数の生育ステージにおける植生指数を同時に考慮することで、単一の生育ステージの情報だけでは識別できない欠株部分と正常株のピクセルを、適切に識別できるようになる。
【０１７７】
なお、本実施例では正常株のピクセルのクラスタは一つが特定されるものとしたが、二つ以上の複数のクラスタが正常株のピクセルのクラスタとして特定されてもよい。
【０１７８】
ステップＳ８−７において、特定したクラスタＣｋの要素数Ｍｋをカウントし、予測対象作付における推定正常株率Ｒを、Ｒ=Ｍｋ/Ｍにより算出する。
【０１７９】
Fuzzy c-means 法では、MkはクラスタＣｋに対する各ベクトルの帰属度の総和であり、以下の数式で算出される。
【０１８０】
【数１５】
0006639300.tif 000016
【０１８１】
ステップＳ８−８において、推定正常株率テーブル１１−９に、予測対象作付の作付ＩＤを付加して、推定正常株率Ｒを登録する。
この処理までを行うのが、正常株率推定装置である。
【０１８２】
＜収穫量予測部２１−９：Ｓ９＞
収穫量予測部２１−９では、生育度推定部２１−５で推定した推定生育度と、正常株率推定部２１−８で推定した正常株率を説明変数に含む収穫量予測モデルで、予測対象作付の収穫量を予測する。一例として、以下のような数式の収穫量予測モデルで、収穫量Ｙを予測することができる。
【０１８３】
【数１６】
0006639300.tif 000017
【０１８４】
収穫量予測モデルのパラメータ a, b, c, d は、過去に実施した作付の実績収穫量をもとに推定することができる。過去に実施した作付の情報が作付情報テーブル１１−１に登録されている場合、上記と同様の方法で、過去に実施した作付に対しても式（１６）による収穫量予測値を算出することができる。このようにして算出した収穫量予測値と、実績収穫量テーブル１１−１１に登録されている実績収穫量との誤差の二乗和が最小になるようにパラメータを学習することで、収穫量予測モデルのパラメータを推定することができる。
【０１８５】
ここで、a, b, c, d は収穫量予測モデルのパラメータであり、Nは作付情報テーブル１１−１から取得した予測対象作付の作付株数、Rは正常株率推定部で推定した正常株率、Vは生育度推定部で推定した推定生育度、Tは気象情報取得手段から取得した最寄観測所における平均気温の生育度集計対象の生育ステージでの合計値、Sは気象情報取得手段から取得した最寄観測所における日射量の生育度集計対象の生育ステージでの合計値である。
【０１８６】
尚、将来の平均気温や日射量は、気象情報取得手段１１により取得した、１か月予報や３か月予報などの長期予報と過去の気象観測結果から予測値を算出して使用することができる。予測値の算出方法の、一例を以下に示す。
【０１８７】
長期予報が「低い（少ない）」という予報であれば、例えば過去３０年の観測値を小さい順に並べて、小さいほうから順に１〜１０番目の平均値を取る。
「平年並み」という予報であれば、小さいほうから順に１１〜２０番目の平均値を取る。
「平年より高い（多い）」という予報であれば、小さいほうから順に２１〜３０番目の平均値を取る。
求めた予測収穫量の値Ｙを、予測収穫量テーブル１１−１０に格納する。
【０１８８】
＜予測収穫量出力部２１−１０：Ｓ１０＞
予測収穫量出力部２１−１０は、予測収穫量テーブル１１−１０に基づいて、出力装置５に予測された収穫量を出力することができる。
【０１８９】
図１２に示す実績収穫量テーブル１１−１１は、過去の作付における実績収穫量の情報が登録されているテーブルである。
【０１９０】
例えば、作付ＩＤ、圃場ＩＤ、作付開始日、当該作付の収穫が実施された日付である収穫実施日、当該作付の収穫量の実績値である実績収穫量を、図１１の予測収穫量と比較等することができる。
【０１９１】
以上に説明したように、本発明の実施の形態によれば正常株率推定部で推定した正常株率を説明変数として収穫量の予測を行う収穫量予測部を有しているため、欠株による収穫量の減少を考慮した予測ができ、高い精度で収穫量を予測することができる。
【０１９２】
また、正常株率推定部においては、クラスタリングによる教師なし学習を使用して正常株率を推定することができるため、分類モデルを学習するために教師データを用意する必要がない。従って、そのための事前の圃場調査や衛星画像の分析を行うことなく実施することができるという利点がある。
【０１９３】
上記の処理および制御は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）やＧＰＵ（Ｇｒａｐｈｉｃｓ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）によるソフトウェア処理、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）やＦＰＧＡ（Ｆｉｅｌｄ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）によるハードウェア処理によって実現することができる。
【０１９４】
また、上記の実施の形態において、添付図面に図示されている構成等については、これらに限定されるものではなく、本発明の効果を発揮する範囲内で適宜変更することが可能である。その他、本発明の目的の範囲を逸脱しない限りにおいて適宜変更して実施することが可能である。
【０１９５】
また、本発明の各構成要素は、任意に取捨選択することができ、取捨選択した構成を具備する発明も本発明に含まれるものである。
【０１９６】
例えば、本発明は、収穫量の予測の基になる正常株の推定技術に利用することも可能である。
【０１９７】
また、本実施の形態で説明した機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各部の処理を行ってもよい。尚、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
【０１９８】
また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１９９】
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また前記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。機能の少なくとも一部は、集積回路などのハードウェアで実現しても良い。
【産業上の利用可能性】
【０２００】
本発明は、農作物の収穫量予測装置として利用可能である。
【符号の説明】
【０２０１】
Ａ…収穫量予測装置、１…データベース群（補助記憶装置）、Ｃ…各種インタフェース部、３…入力装置、５…出力装置、７…中央演算装置（ＣＰＵ）（制御部）、１１…気象情報取得手段、２１…主記憶装置、２１−１…予測対象作付設定部、２１−２…低解像度画像データ分析部、２１−３…生育モデル適合部、２１−４…生育ステージ推定部、２１−５…生育度推定部、２１−６…高解像度画像データ選択部、２１−７…高解像度画像データ分析部、２１−８…正常株率推定部、２１−９…収穫量予測部、２１−１０…予測収穫量出力部。
【図１Ａ】
0006639300.tif 000018
【図１Ｂ】
0006639300.tif 000019
【図２】
0006639300.tif 000020
【図３】
0006639300.tif 000021
【図４】
0006639300.tif 000022
【図５】
0006639300.tif 000023
【図６】
0006639300.tif 000024
【図７】
0006639300.tif 000025
【図８】
0006639300.tif 000026
【図９】
0006639300.tif 000027
【図１０】
0006639300.tif 000028
【図１１】
0006639300.tif 000029
【図１２】
0006639300.tif 000030
【図１３】
0006639300.tif 000031
【図１４】
0006639300.tif 000032
【図１５】
0006639300.tif 000033
【図１６】
0006639300.tif 000034
【図１７Ａ】
0006639300.tif 000035
【図１７Ｂ】
0006639300.tif 000036
【図１８】
0006639300.tif 000037
【図１９】
0006639300.tif 000038
【図２０】
0006639300.tif 000039
【図２１】
0006639300.tif 000040
【図２２】
0006639300.tif 000041
【図２３】
0006639300.tif 000042
【図２４】
0006639300.tif 000043
【図２５】
0006639300.tif 000044
【図２６】
0006639300.tif 000045
【図２７Ａ】
0006639300.tif 000046
【図２７Ｂ】
0006639300.tif 000047