JP2022103139A - 作物の生育診断方法、及び作物の生育診断プログラム - Google Patents

Abstract

【課題】作物の生育状態を迅速に判定し、出力する。【解決手段】第１指標値算出部３２は、同一地点（区画）の作物に対して実施した２回の生体調査の結果（画像から得らえる投影葉面積）に基づいて、作物の生育速度を示す第１指標値（Ｉｄ）を算出する（Ｓ１６）。また、第２指標値算出部３６は、２回の生体調査の結果と、生体調査の間の期間における日射量（Ｓｒn）と、に基づいて、日射による生育への影響を示す第２指標値（光利用効率ＬＵＥ1→2）を算出する（Ｓ１８～Ｓ２４）。そして、生育状態判定部３８は、第１指標値（Ｉｄ）と第２指標値（ＬＵＥ1→2）とに基づいて、作物の生育が順調か否か及び作物の生育と日射との関係を判定し（Ｓ２６）、出力部４２は、判定結果を出力する（Ｓ３６）。【選択図】図４

Description

本発明は、作物の生育診断方法、及び作物の生育診断プログラムに関する。

作物の収量を最大限に得るには、作物を順調に生育させることが重要である。しかし、現状においては、作物の生育が順調であるかは、生産者が経験や勘により判断していることがほとんどである。一方、生育診断として汁液等を採取して養分濃度を計測する技術はあるものの、専用の分析機器が必要である。

なお、従来、生育状態に影響する環境因子の変動データの相関から生育状態を診断する技術が知られている（例えば、特許文献１等参照）。

特開２０１７－０５１１２５号公報

しかしながら、従来においては、作物の生育状態を迅速に判断することは難しい。

本発明は、作物の生育状態を迅速に判定し、出力することが可能な作物の生育診断方法、及び作物の生育診断プログラムを提供することを目的とする。

本発明の作物の生育診断方法は、同一地点の作物に対して実施した２回の生体調査の結果に基づいて、前記作物の生育速度を示す第１指標値を算出し、前記２回の生体調査の結果と、前記生体調査の間の期間における日射に関する情報と、に基づいて、日射による生育への影響を示す第２指標値を算出し、前記第１指標値と前記第２指標値とに基づいて、前記作物の生育が順調か否か及び前記作物の生育と日射との関係を判定し、前記判定する処理の判定結果を出力する、処理をコンピュータが実行する作物の生育診断方法である。

本発明の作物の生育診断方法、及び作物の生育診断プログラムは、作物の生育状態を迅速に判定し、出力することができるという効果を奏する。

一実施形態に係る農業システムの構成を示す図である。 図１のユーザ端末のハードウェア構成を示す図である。 図１のユーザ端末の機能ブロック図である。 ユーザ端末の処理を示すフローチャートである。 図４のステップＳ３０の詳細処理を示すフローチャートである。 カメラが撮影した画像の一例を示す図である。 投影葉面積と乾物の関係を示すグラフである。 図８（ａ）は、ｎ日目の投影葉面積を求める方法について説明するための図であり、図８（ｂ）は、投影葉面積と葉面積指数の関係を示すグラフである。 生育状態判定部が判定に用いるテーブルの一例を示す図である。 図１０（ａ）～図１０（ｃ）は、スコア化テーブルの例を示す図である。 図１１（ａ）～図１１（ｃ）は、スコア化テーブルを用いて気温スコア、湿度スコア、土壌水分スコアを求めた例を説明するための図である。 図１２（ａ）は、生育状態表示画面の一例を示す図であり、図１２（ｂ）は、抑制要因表示画面の一例を示す図である。 変形例に係るユーザ端末の処理を示すフローチャートである。 変形例における第１指標値を用いた判定結果の一例を示す図である。 図１５（ａ）、図１５（ｂ）は、１回目及び２回目の撮影により得られた画像であり、図１５（ｃ）、図１５（ｄ）は、１回目及び２回目の撮影により得られた画像を二値化したものである。 １回目の撮影画像と２回目の撮影画像とから得られた情報を示す図である。 閾値を求める際に用いた情報と、閾値を示す図である。 図１８（ａ）～図１８（ｄ）は、気温スコアと土壌水分スコアについて説明するための図である。 図１９（ａ）は、気温と土壌水分の期間平均値を示す図であり、図１９（ｂ）は、気温スコアと土壌水分スコアを示す図である。 図２０（ａ）～図２０（ｃ）は、生育抑制の要因の特定方法について説明するための図である。

以下、農業システムの一実施形態について、図１～図１２に基づいて詳細に説明する。図１には、一実施形態に係る農業システム１００の構成が概略的に示されている。本実施形態の農業システム１００は、ブロッコリーやキャベツなどの露地野菜を栽培する生産者（以下、ユーザと呼ぶ）が利用可能なシステムである。

農業システム１００は、図１に示すように、ユーザ端末１０と、カメラ７０と、センサ７２と、を備える。ユーザ端末１０と、カメラ７０と、センサ７２は、インターネットやＬＡＮ（Local Area Network）などのネットワーク８０に接続されており、ユーザ端末１０は、カメラ７０によって撮像された画像やセンサ７２によって検出された情報を取得することが可能となっている。

カメラ７０は、圃場に設けられ、作物（株）を上方（所定高さ）から撮影する。カメラ７０が撮影した画像は、ユーザ端末１０に送信される。なお、カメラ７０は、１つでもよいし、複数でもよい。

センサ７２は、圃場の複数区画における日射量や土壌水分を検出するセンサ、圃場の気温や湿度を検出するセンサ等を含む。センサ７２が検出した情報は、ユーザ端末１０に送信される。なお、各区画の土壌水分については、例えば、降水量や気温から検出される圃場全体の平均的な土壌水分を、画像から得られる各区画の土壌の色や、表面温度、電導度等に基づいて補正することで、得ることができる。

ユーザ端末１０は、圃場に設置されたカメラ７０により所定時間間隔で撮影された画像から得られる作物（ブロッコリーとする）の葉に関する情報や、センサ７２により検出された情報を取得する。また、ユーザ端末１０は、取得した情報に基づいて、生育診断（作物の生育が順調であるか否かや、順調でなければその原因が何であるのかの判定）を行い、診断結果（判定結果）を出力する。

図２には、ユーザ端末１０のハードウェア構成が概略的に示されている。図２に示すように、ユーザ端末１０は、ＣＰＵ（Central Processing Unit）９０、ＲＯＭ（Read Only Memory）９２、ＲＡＭ（Random access memory）９４、記憶部（例えば、ＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）など）９６、ネットワークインタフェース９７、表示部９３、入力部９５、及び可搬型記憶媒体用ドライブ９９等を備えている。表示部９３は、液晶ディスプレイ等を含み、入力部９５は、キーボードやマウス、タッチパネル等を含む。これらユーザ端末１０の構成各部は、バス９８に接続されている。ユーザ端末１０では、ＲＯＭ９２あるいはＨＤＤ９６に格納されているプログラム（作物の生育診断プログラムを含む）、或いは可搬型記憶媒体用ドライブ９９が可搬型記憶媒体９１から読み取ったプログラム（作物の生育診断プログラムを含む）をＣＰＵ９０が実行することにより、図３に示す各部の機能が実現される。なお、図３の各部の機能は、例えば、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等の集積回路により実現されてもよい。

図３には、ユーザ端末１０の機能ブロック図が示されている。ユーザ端末１０においては、ＣＰＵ９０がプログラムを実行することにより、図３に示すように、画像取得部３０、第１指標値算出部３２、センサ情報取得部３４、第２指標値算出部３６、生育状態判定部３８、抑制要因特定部４０、出力部４２、としての機能が実現されている。なお、図３には、生育状態判定部３８による判定結果や抑制要因特定部４０の特定結果を格納する判定結果テーブル４８も図示されている。

画像取得部３０は、カメラ７０が作物（株）の直上の任意の高さから撮影した株の投影画像を取得し、第１指標値算出部３２に受け渡す。

第１指標値算出部３２は、所定時間間隔をあけて同一の株を撮影した２枚の画像（すなわち、圃場の同一の区画を撮影した２枚の画像）を取得し、各画像を二値化したり、ＡＩ解析等を行うことで、各画像の土壌の部分と葉の部分を区別する。また、第１指標値算出部３２は、各画像における葉の部分の面積、すなわち株の投影葉面積を算出する。なお、第１指標値算出部３２が算出した各画像の投影葉面積の値は、２回の生体調査の結果であるといえる。更に、第１指標値算出部３２は、２回の生体調査の結果（投影葉面積）を用いて、作物の生育速度を示す指標（第１指標値）を算出する。

センサ情報取得部３４は、センサ７２が検出した、圃場の複数区画における日射量や土壌水分、圃場の気温や湿度の値を取得し、第２指標値算出部３６や抑制要因特定部４０に受け渡す。

第２指標値算出部３６は、第１指標値算出部３２が算出した投影葉面積を用いて、２回の生体調査を行ったときの乾物重や２回の生体調査の間における各日の葉面積指数（ＬＡＩ：Leaf Area Index）を推定する。また、第２指標値算出部３６は、推定した乾物重や葉面積指数、センサ情報取得部３４が取得した日射量等を用いて、第２指標値として、２回の生体調査の間における光利用効率（ＬＵＥ：Light Use Efficiency）を算出する。

生育状態判定部３８は、第１指標値（作物の生育速度を示す指標）と第２指標値（光利用効率）とを用いて、作物（株）の生育が順調であるか否か、順調でないならばその要因が日射であるかを判定する。

抑制要因特定部４０は、生育状態判定部３８の判定の結果、作物（株）の生育が順調でなく、その要因が日射量でない場合に、センサ情報取得部３４が取得した、各区画における土壌水分や、圃場の気温や湿度の値を用いて、抑制要因を特定する。

なお、生育状態判定部３８の判定結果や、抑制要因特定部４０の特定結果は、判定結果テーブル４８に格納される。

出力部４２は、判定結果テーブル４８に格納されている情報（生育状態判定部３８の判定結果や、抑制要因特定部４０の特定結果）を表示する画面を生成し、表示部９３上に出力する。

（ユーザ端末１０の処理について）
以下、ユーザ端末１０の処理について、図４、図５のフローチャートに沿って、その他図面を適宜参照しつつ、詳細に説明する。なお、本実施形態では、圃場で生育しているブロッコリーの生育が順調であるか否か、順調でなければその要因が何であるのかを判定し、判定結果を出力する例について説明する。

図４の処理の前提として、予めユーザによって、画像の撮影タイミング（１回目、２回目）が設定されているものとする。なお、１回目の撮影と２回目の撮影の間は、所定日数空いているものとする。ユーザは、作物の定植日や播種日等に基づいて、撮影タイミングを決定する。また、センサ情報取得部３４は、センサ７２が所定タイミングで検出する各種情報を適宜取得しているものとする。

図４の処理では、まず、ステップＳ１０において、画像取得部３０が、各カメラ７０によって撮影された画像（１回目）を取得する。このステップＳ１０では、図６に示すような画像が取得できるようになっている。なお、図６は、圃場内の複数のカメラ７０によって撮影された画像を並べたものであり、圃場全体を撮影した画像であるといえる。この場合、１つのカメラ７０で複数の区画（破線で区切られた範囲）を撮影してもよいし、１つのカメラ７０で１つの区画を撮影してもよい。また、１つのカメラ７０で圃場全体を撮影してもよい。いずれの場合であっても、画像取得部３０は、図６のような圃場全体を撮影した画像を、所定の方法で複数の区画に分割する。なお、各区画は、概ね１つの株が含まれる範囲である。

次いで、ステップＳ１２では、画像取得部３０が、各カメラ７０によって撮影された画像（２回目）を取得する。この場合も、画像取得部３０は、図６と同様の画像を取得する。

次いで、ステップＳ１４では、第１指標値算出部３２が、１区画（１つの区画）を特定する。例えば、第１指標値算出部３２は、図６の左上の区画を特定するものとする。

次いで、ステップＳ１６では、第１指標値算出部３２が、特定した区画の株の生育速度を示す指標である第１指標値（Ｉｄ）を算出する。具体的には、第１指標値算出部３２は、まず、特定した区画を１回目に撮影した画像を二値化したり、ＡＩ解析等を行うことで、土壌と葉の部分を区別し、葉の部分の面積（投影葉面積）Ａ₁を算出する。次いで、第１指標値算出部３２は、特定した区画を２回目に撮影した画像を二値化したり、ＡＩ解析等を行うことで、土壌と葉の部分を区別し、葉の部分の面積（投影葉面積）Ａ₂を算出する。そして、第１指標値算出部３２は、投影葉面積Ａ₁、Ａ₂と、第１回目の撮影と第２回目の撮影との間の日数ｄを用いて、次式（１）より、第１指標値Ｉｄを算出する。
Ｉｄ＝（ｌｎＡ₂－ｌｎＡ_１）／ｄ …（１）

なお、本実施形態では、撮影タイミングが、ブロッコリーの投影葉面積がシグモイド曲線に沿って大きくなる時期（ステージ）であるため、上式（１）のような式を用いたが、これに限られるものではない。作物の種類や生育ステージに応じて、式を適宜変更してもよい。

次いで、ステップＳ１８では、第２指標値算出部３６が、特定した区画の、１回目の撮影と２回目の撮影の間の各日（ｎ日）の日射量Ｓｒ_nをセンサ情報取得部３４から取得する。

次いで、ステップＳ２０では、第２指標値算出部３６が、１回目の撮影の際と２回目の撮影の際の乾物重Ｄｍ₁、Ｄｍ₂を推定する。具体的には、第２指標値算出部３６は、次式（２）に基づいて、乾物重Ｄｍ₁、Ｄｍ₂を算出する。なお、ｘは撮影回数を示し、１又は２である。また、ａ、ｂは係数である。
Ｄｍ_x＝ａ・Ａ_x＋ｂ …（２）

例えば、ブロッコリーの場合、予め行った試験の結果、投影葉面積Ａと乾物重Ｄｍとの間には、図７に示すような関係が見出されている。したがって、第２指標値算出部３６は、上式（２）の係数をａ＝２０１．６、ｂ＝－５．０として、投影葉面積Ａ₁、Ａ₂を上式（２）に代入することで、乾物重Ｄｍ₁、Ｄｍ₂を精度よく算出することができる。なお、投影葉面積は葉の重なりが増えてくると、ＬＡＩとの間に直線関係が失われてくるため、０を最小値とする二次曲線で近似したほうがよい場合もある。このような場合には、上式（２）に代えて、二次曲線の式を用いることとしてもよい。

更に、ステップＳ２０では、第２指標値算出部３６は、１回目の撮影と２回目の撮影の間の各日の葉面積指数ＬＡＩ（Ｌａ_n）を推定する。具体的には、第２指標値算出部３６は、次式（３）に基づいて、葉面積指数Ｌａ_nを算出する。なお、ｃ、ｄは係数である。
Ｌａ_n＝ｃ・Ａ_n＋ｄ …（３）

ここで、Ａ_nは、１回目の撮影と２回目の撮影の間のｎ日目の投影葉面積を意味する。ｎ日目の投影葉面積Ａ_nを求める場合、まず、図８（ａ）に示すように、１回目の撮影画像から得られた投影葉面積Ａ₁と２回目の撮影画像から得られた投影葉面積Ａ₂の間を直線で補完して一次直線の推定式を作成する。そして、１日あたりの投影葉面積増加量（（Ａ₂－Ａ₁）／ｄ）とｎの積からｎ日間の投影葉面積増加量を求め、１回目の投影葉面積Ａ₁に加算することで、ｎ日目の投影葉面積（Ａ_n）を求めることができる。これを式に表すと、次式（４）のようになる。
Ａ_n＝Ａ₁＋｛（Ａ₂－Ａ₁）／ｄ｝×ｎ …（４）

なお、図８（ａ）では、一次直線の推定式を作成する場合について説明したが、生育状況によっては、１回目と２回目の投影葉面積の間を二次曲線などで補完する場合もある。

例えば、ブロッコリーの場合、予め行った試験の結果、投影葉面積Ａと葉面積指数ＬＡＩとの間には、図８（ｂ）に示すような関係が見出されている。したがって、第２指標値算出部３６は、上式（３）の係数をｃ＝３．７２、ｃ＝０とし、Ａ_nに上式（４）で得られる値を代入する。これにより、ｎ日目の葉面積指数Ｌａ_nを精度よく算出することができる。

図４に戻り、次のステップＳ２２では、第２指標値算出部３６が、推定したｎ日目の葉面積指数Ｌａ_nと、ｎ日目の日射量Ｓｒ_nから、ｎ日目における株の受光量Ｌｉ_nを推定する。ここでは、第２指標値算出部３６は、次式（５）を用いることで、ｎ日目における株の受光量Ｌｉ_nを推定する。
Ｌｉ_n＝（１－ｅ^{-k Lan}）・Ｓｒ_n …（５）

次いで、ステップＳ２４では、第２指標値算出部３６が、１回目の撮影と２回目の撮影の間（撮影間）の総受光量と、撮影間の乾物重の変化とから、撮影間の光利用効率ＬＵＥ（第２指標値）を算出する。

具体的には、第２指標値算出部３６は、撮影間の総受光量Ｌｉ_1→2を次式（６）から求める。
Ｌｉ_1→2＝ΣＬi_n …（６）

なお、上式（６）のΣＬｉ_nは、撮影間の各日の受光量Ｌiの合計を意味する。

そして、第２指標値算出部３６は、撮影間の総受光量Ｌｉ_1→2と、撮影間の乾物重の変化（Ｄｍ₂－Ｄｍ₁）を用いて、次式（７）から撮影間の第２指標値（光利用効率）ＬＵＥ_1→2を算出する。
ＬＵＥ_1→2＝（Ｄｍ₂－Ｄｍ₁）／Ｌｉ_1→2 …（７）

次いで、ステップＳ２６では、生育状態判定部３８が、特定した区画の作物（株）の生育が順調か否かを判定する。具体的には、生育状態判定部３８は、図９に示すようなテーブルに基づいて、判定を実行する。ここで、図９のＴｈ_Idは、第１指標値Ｉｄの大小を判断する閾値であり、Ｔｈ_LUEは、第２指標値ＬＵＥ_1→2の大小を判断する閾値である。生育状態判定部３８は、第１指標値Ｉｄと閾値Ｔｈ_Idの大小関係と、第２指標値ＬＵＥ_1→2と閾値Ｔｈ_LUEの大小関係とから、以下のように生育状態の判定を実行する。

ｉ）Ｉｄ≧Ｔｈ_IdかつＬＵＥ_1→2≧Ｔｈ_LUEのとき、
生育は順調で日射量に見合った生育を示している。
ii)Ｉｄ＜Ｔｈ_IdかつＬＵＥ_1→2≧Ｔｈ_LUEのとき、
生育は抑制されているが、その抑制要因は主に日射量不足である。
iii）Ｉｄ＜Ｔｈ_IdかつＬＵＥ_1→2＜Ｔｈ_LUEのとき、
生育は抑制されているが、日射量は十分あるため、他の抑制要因がある。
iv）Ｉｄ≧Ｔｈ_IdかつＬＵＥ_1→2＜Ｔｈ_LUEのとき、
何らかの抑制要因により、日射量に見合った程の生育は得られておらず、生育は抑制されている。たとえば、生育は抑制されているものの、日射量の多さで生育が補償されている状態である。

次いで、ステップＳ２８では、生育状態判定部３８が、生育が順調であるか否かを判断する。すなわち、上記ｉ）と判定された場合には、ステップＳ２８の判断は肯定され、ステップＳ３２に移行する。一方、上記ｉ）以外と判定された場合、ステップＳ２８の判断は否定され、ステップＳ３０に移行する。

生育が抑制されている場合（ステップＳ２８の判断が否定された場合）、ステップＳ３０に移行し、抑制要因特定部４０は、抑制要因特定処理を実行する。具体的には、抑制要因特定部４０は、図５のフローチャートに沿った処理を実行する。

図５の処理においては、まず、ステップＳ５０において、抑制要因特定部４０が、抑制要因が日射量であるか否かを判断する。すなわち、抑制要因特定部４０は、ステップＳ２６において上記ii）と判定されたか否かを判断する。このステップＳ５０の判断が肯定された場合には、図５の全処理を終了して、図４のステップＳ３２に移行する。一方、ステップＳ５０の判断が否定された場合には、ステップＳ５２に移行する。

ステップＳ５２に移行すると、抑制要因特定部４０は、センサ情報取得部３４から複数種類の環境要素の値を取得する。本実施形態では、例えば撮影間の期間における圃場の期間平均気温、圃場の期間平均相対湿度、特定した区画の期間平均土壌水分を取得する。

次いで、ステップＳ５４では、抑制要因特定部４０が、取得した複数種類の環境要素をスコア化する。スコア化においては、図１０（ａ）～図１０（ｃ）に示すようなスコア化テーブルを用いる。

図１０（ａ）のスコア化テーブルは、期間平均気温とスコア（気温スコア）との関係を示すテーブルである。図１０（ａ）のテーブルは、以下のようにして、作物の品種ごとに事前に取得しているものとする。例えば、環境制御が可能なチャンバ内において、基準気温（例えば２２℃）に設定したうえで、作物を所定時間生育させ、そのときの作物の重量変化を基準（気温スコア＝１００％）として取得する。そして、チャンバ内を他の気温に設定して作物を所定時間生育させたときの作物の重量変化を測定し、当該重量変化の基準に対する割合（百分率）を当該他の気温における気温スコアとして図１０（ａ）の座標系上にプロットする。この処理をチャンバ内の気温を異ならせながら繰り返すことで、気温スコアを複数プロットし、最小二乗法等を用いて直線近似等することで、図１０（ａ）の関係を取得する。なお、図１０（ａ）の関係を取得する手順は、上記手順に限られるものではない。

なお、図１０（ｂ）のスコア化テーブルは、期間平均相対湿度とスコア（湿度スコア）との関係を示すテーブルである。また、図１０（ｃ）のスコア化テーブルは、期間平均土壌水分とスコア（土壌水分スコア）との関係を示すテーブルである。これら、図１０（ｂ）、図１０（ｃ）のテーブルについても、図１０（ａ）と同様にして作成することができる。

抑制要因特定部４０は、以下のａ）～ｄ）のように期間平均気温Ｔ（℃）を場合分けして、気温スコアＳｃ_Tを算出する。

ａ) Ｔ＜６℃、またはＴ＞４０℃のとき
抑制要因特定部４０は、気温スコアＳｃ_T＝０（％）とする。

b) ６℃≦Ｔ＜１２℃のとき
抑制要因特定部４０は、次式（８）から、気温スコアＳｔを求める。
Ｓｃ_T＝１６．７×Ｔ－１００（％） …（８）

c) １２℃≦Ｔ＜２８℃のとき
抑制要因特定部４０は、気温スコアＳｃ_T＝１００（％）とする。

d) ２８℃≦Ｔ<４０℃のとき
抑制要因特定部４０は、次式（９）から、気温スコアＳｔを求める。
Ｓｃ_T＝－８．３×Ｔ＋３３３（％） …（９）

なお、抑制要因特定部４０は、上記と同様に、図１０（ｂ）のテーブルに期間平均相対湿度を当てはめて、湿度スコアＳｃ_Hを求めるとともに、図１０（ｃ）のテーブルに期間平均土壌水分を当てはめて、土壌水分スコアＳｃ_Wを求める。

以下、２つの例について説明する。例えば、例１（期間平均気温＝２２．７℃、期間平均相対湿度＝７７．５％、期間平均土壌水分＝４４．８ｍ³／ｍ³）であれば、図１１（ａ）～図１１（ｃ）において破線で示すように、気温スコアＳｃ_T＝１００％、湿度スコアＳｃ_H＝１００％、土壌水分スコアＳｃ_W＝１００％と求められる。また、例２（期間平均気温＝２０．８℃、期間平均相対湿度＝６６．１％、期間平均土壌水分＝３９．９ｍ³／ｍ³）であれば、図１１（ａ）～図１１（ｃ）において一点鎖線で示すように、気温スコアＳｃ_T＝１００％、湿度スコアＳｃ_H＝１００％、土壌水分スコアＳｃ_W＝６０％と求められる。

次いで、ステップＳ５６では、抑制要因特定部４０が、スコアに基づいて抑制要因を特定する。

抑制要因特定部４０は、例えば、次式（１０）にステップＳ５４で得られたスコアを代入することで、ＬＵＥ_1→2を算出する。
ＬＵＥ_1→2＝ＬＵＥ_p×ｆ（Ｓｃ_T×Ｓｃ_H×Ｓｃ_W） …（１０）

なお、ＬＵＥ_pは、ポテンシャルＬＵＥであり、生育抑制のない場合の最大ＬＵＥを意味する。

そして、抑制要因特定部４０は、上式（１０）においてスコアＳｃ_T、Ｓｃ_H、Ｓｃ_Wのそれぞれの影響を順次除し、ＬＵＥ_1→2≒ＬＵＥ_pとなったときに除していたスコア（環境要素）を抑制要因と特定する。例えば、図１１の例２の場合、気温スコアＳｃ_T＝１００％＝１、湿度スコアＳｃ_H＝１００％＝１、土壌水分スコアＳｃ_W＝６０％＝０．６であるので、上式（１０）は次式（１０）’のようになるが、
ＬＵＥ_1→2＝ＬＵＥ_p×ｆ（１×１×０．６） …（１０）’

このうち、土壌水分スコアＳｃ_W＝０．６の影響を除したときに、ＬＵＥ_1→2≒ＬＵＥ_pとなるのであれば、土壌水分が抑制要因であると特定することができる。

一方、抑制要因特定部４０は、スコアＳｃ_T、Ｓｃ_H、Ｓｃ_Wのそれぞれの影響を順次除しても、ＬＵＥ_1→2＜ＬＵＥ_pとなった場合には、その他の抑制要因があると特定する。その他の抑制要因には、例えば、肥料養分の不足や病虫害の被害が発生していること、などが挙げられる。

以上のようにして図５の処理（ステップＳ３０）が行われた後は、図４のステップＳ３２に移行する。なお、ステップＳ３２に移行する段階では、特定されている区画の生育が順調であるか否かが判定され、成長が抑制されている場合には、その抑制要因（日射量、気温、湿度、土壌水分、その他）が特定されているので、判定結果テーブル４８には、これらの判定結果及び特定結果が、各区画の位置情報と対応付けて格納される。

図４のステップＳ３２に移行すると、生育状態判定部３８は、全区画の処理が終了したか否かを判断する。このステップＳ３２の判断が否定されると、ステップＳ３４に移行し、第１指標値算出部３２が、圃場内の次の区画を特定する。その後は、ステップＳ１６に戻り、上述した処理を繰り返す。一方、ステップＳ３２の判断が肯定されると、ステップＳ３６に移行する。

ステップＳ３６に移行すると、出力部４２は、判定結果テーブル４８に格納されている情報に基づいて画面を生成し、表示部９３に出力する。例えば、出力部４２は、図１２（ａ）に示すように、圃場の区画ごとに、生育状態判定部３８の判定結果を色分けして示す生育情報表示画面を生成し、表示部９３に出力してもよい。この場合、出力部４２は、各区画にマウスカーソルが合わせられたときに、各区画の情報（例えば、どの環境要素に注意すべきかや、各スコアの値など）を表示するようにしてもよい。なお、日射量不足の箇所は、建物や樹木によって日陰になっている箇所であると考えられる。また、例えば、ユーザが土壌水分スコアの表示要求を入力した場合には、図１２（ｂ）に示すような各区画の土壌水分スコアを表示する画面（抑制要因表示画面）を表示部９３に出力してもよい。なお、土壌水分スコア以外の表示要求があった場合にも、図１２（ｂ）と同様の画面を出力すればよい。

以上、詳細に説明したように、本実施形態によると、第１指標値算出部３２は、同一地点（区画）の作物に対して実施した２回の生体調査の結果（２枚の撮影画像それぞれから得られる投影葉面積）に基づいて、作物の生育速度を示す第１指標値（Ｉｄ）を算出する（Ｓ１６）。また、第２指標値算出部３６は、２回の生体調査の結果と、生体調査の間の期間における日射量（Ｓｒ_n）と、に基づいて、日射による生育への影響を示す第２指標値（光利用効率ＬＵＥ_1→2）を算出する（Ｓ１８～Ｓ２４）。そして、生育状態判定部３８は、第１指標値（Ｉｄ）と第２指標値（ＬＵＥ_1→2）とに基づいて、作物の生育が順調か否か及び作物の生育と日射との関係を判定し（Ｓ２６）、出力部４２は、生育状態判定部３８による判定結果を表示部９３に出力する（Ｓ３６）。これにより、本実施形態では、専用の分析機器等を用いずに、作物の生育状態を迅速に判定し、出力することができる。また、従来のように経験や勘に基づいて生育状態を判断する場合と比べ、精度よく作物の生育状態を判定することができる。

また、本実施形態では、生育状態判定部３８は、第１指標値（Ｉｄ）と予め定めた第１閾値（Ｔｈ_Id）との大小関係と、第２指標値（ＬＵＥ_1→2）と予め定めた第２閾値（Ｔｈ_LUE）との大小関係と、に基づいて、図９のテーブルから、作物の生育が順調か否か及び作物の生育と日射との関係を判定する。これにより、作物の生育状態の判定処理を簡素化することができる。

また、本実施形態では、生育状態判定部３８は、第２指標値が第２閾値よりも小さい場合（ステップＳ２８：否定、かつステップＳ５０：否定）に、２回の生体調査の間の期間における複数種類の環境要素（気温、湿度、土壌水分等）に関する情報を取得し（Ｓ５２）、取得した情報に基づいて、作物の生育に対して影響を与えている環境要素の種類を特定する（Ｓ５４、Ｓ５６）。これにより、作物の生育に対して影響を与えている環境要素の種類を実際の環境要素の値に基づいて精度よく特定することができる。また、第２指標値が第２閾値よりも大きい場合には、ステップＳ５２～Ｓ５６の処理を実行しないため、処理の簡素化を図ることができる。

また、本実施形態では、作物を所定高さから撮影した画像から作物の葉に関する情報（投影要面積）を取得することにより、生体調査を行う。これにより、作物に直接触れず、かつ作物を破壊することなく、生体調査を行うことができる。また、画像を用いることで、ユーザが葉の大きさ等を計測しなくてもよいため、ユーザの手間を軽減することができる。

また、本実施形態では、出力部４２は、圃場内の複数の作物（株）に関する判定結果を、作物（株）の位置と関連付けて表示する画面（図１２（ａ））を生成し、出力する。これにより、ユーザは、どの位置において、どのように環境を調整すべきかを容易に把握することができる。また、ユーザは、生育状態を把握することで、収量を適切に予測することができる。なお、ユーザ端末１０は、本実施形態により得られる生育状態判定部３８の判定結果や抑制要因特定部４０の特定結果を用いて収量予測を行い、予測結果を表示部９３に出力するようにしてもよい。

なお、上記実施形態では、出力部４２は、図１２（ａ）などの画面を生成し、表示部９３に出力する場合について説明したが、これに限られるものではない。出力部４２は、生育状態判定部３８による判定結果や、抑制要因特定部４０による特定結果に基づいて、環境制御装置に対して制御指示を出力することとしてもよい。環境制御装置は、例えば屋外の圃場であれば、給水装置などである。また、温室内であれば、環境制御装置は、ヒートポンプ、換気窓、暖房機、遮光・保温カーテン、細霧装置等である。

なお、上記実施形態では、生育状態判定部３８は、第２指標値（ＬＵＥ_1→2）と閾値（Ｔｈ_LUE）との大小関係に基づいて、生育状態を判定する場合について説明したが、これに限られるものではない。例えば、０＜ＬＵＥ_1→2＜０．３×Ｔｈ_LUEの場合、日射量以外の生育抑制が大きい、０．３×Ｔｈ_LUE≦ＬＵＥ_1→2＜０．７×Ｔｈ_LUEの場合、日射量以外の生育抑制が中くらい、０．７×Ｔｈ_LUE≦ＬＵＥ_1→2＜Ｔｈ_LUEの場合、日射量以外の生育抑制が小さい、というように、複数段階で判定を行うこととしてもよい。

なお、上記実施形態では、ユーザ端末１０は、図４、図５のフローチャートに沿って処理を実行する場合について説明したが、これに限らず、図１３のフローチャートに沿った処理を実行することとしてもよい。図１３の処理は、ステップＳ１０、Ｓ１２、Ｓ１４、Ｓ１６は図４と同様であるが、ステップＳ１６において第１指標値（Ｉｄ）が算出された後に、ステップＳ１７が実行される。このステップＳ１７においては、第２指標値算出部３６が、第１指標値（Ｉｄ）が閾値（Ｔｈ_Id）よりも大きいか否かを判断する。

ステップＳ１７の判断が否定された場合には、第２指標値算出部３６は、ステップＳ１８、Ｓ２０、Ｓ２２、Ｓ２４の処理を実行して、第２指標値（ＬＵＥ_1→2）を算出する。そして、ステップＳ２４の後は、ステップＳ３０に移行し、抑制要因特定部４０が抑制要因特定処理（図５の処理）を実行する。なお、図５のステップＳ５０においては、抑制要因特定部４０は、算出した第２指標値（ＬＵＥ_1→2）が閾値Ｔｈ_LUE以上であるかを判断するものとする。その後、図５の処理（ステップＳ３０）が終了すると、ステップＳ３２に移行する。

一方、ステップＳ１７の判断が肯定された場合には、ステップＳ１８～Ｓ２４、Ｓ３０を経ずに、ステップＳ３２に移行する。このように、図１３の処理では、第１指標値（Ｉｄ）が閾値（Ｔｈ_Id）である場合にのみ、第２指標値（ＬＵＥ_1→2）の算出（Ｓ１８～Ｓ２４）や抑制要因特定処理（Ｓ３０）を行うこととしている。例えば、第１指標値（Ｉｄ）による判定の結果、図１４において黒塗りで示すように、生育が順調と判定された区画については、第２指標値（ＬＵＥ_1→2）を算出せず、図１４においてハッチングを付して示すように、生育が抑制されていると判定された区画についてのみ、第２指標値（ＬＵＥ_1→2）を算出する。このようにすることで、第１指標値Ｉｄと第２指標値ＬＵＥ_1→2を毎回算出する場合と比べ、処理量を低減することができる。なお、図１３の処理の結果としては、図１４に示すような画面を表示することができる。この場合、図１４においてハッチングを付して示す区画にマウスカーソルが合わせられたときに、図１２（ａ）と同様、各区画の情報（例えば、どの環境要素に注意すべきかや、各スコアの値など）を表示するようにしてもよい。

なお、上記実施形態では、撮影した画像を用いて生体調査を行う場合について説明したが、これに限らず、生体調査としてユーザが作物の各種寸法（葉長や葉幅、葉数など）を測定してもよい。また、生体調査により得られる情報は、投影葉面積に限らず、立体的な植物体の大きさや草丈、株はりなどであってもよい。

なお、上記実施形態では、図４、図５の処理をユーザ端末１０が行う場合について説明したが、これに限らず、ユーザ端末１０とネットワーク接続されている情報処理装置（サーバ等）が図４、図５の処理を行うこととしてもよい。また、複数の情報処理装置が、図４、図５の処理を分担して処理することとしてもよい。

なお、上記実施形態では、ブロッコリーやキャベツなどの露地野菜の生育状態を判定する場合について説明したが、これに限らず、トマトなどの施設野菜や、その他の作物の生育状態を判定する場合にも上記実施形態の農業システム１００を利用することができる。

（実施例）
以下、実施例について説明する。本実施例は、上述した実施形態の処理を行うことで、適切に生育診断を行うことができることを検証したものである。本実施例における生育診断の対象作物は、キャベツである。本発明者は、試験圃場（露地）において、慣行法に従い施肥、耕耘を行って、２０２１年３月３１日にキャベツの苗を定植した。キャベツの供試品種は「おきな」である。

（画像の撮影）
本実施例においては、同一のキャベツ個体について１回目の撮影を２０２１年４月１３日に行い、２回目の撮影を２０２１年４月２７日に行った（図４のＳ１０，Ｓ１２に対応）。この撮影により得られた画像が、図１５（ａ）、図１５（ｂ）に示されている。撮影は、対象株の直上１４０ｃｍの高さから行った。撮影の際には、投影葉面積を計算するための既知の大きさ（投影面積）の指標物（図１５（ａ）、図１５（ｂ）ではボール）をキャベツ個体の近傍に置いた。

（第１指標値の算出）
本発明者は、第１指標値Ｉｄを算出するため（図４のステップＳ１６に対応）、図１５（ａ）、図１５（ｂ）の撮影画像を画像解析ソフト（Image J：米国NIH (National Institute of Health)開発のサイエンス向け画像処理ソフト）を用いて二値化した（図１５（ｃ）、図１５（ｄ）参照）。そして、土壌と区別された葉の部分のピクセル数と指標物のピクセル数とを用いることで、葉の投影葉面積（Ａ₁、Ａ₂）を算出した。

投影葉面積Ａ₁、Ａ₂は、図１６に示すように、Ａ₁＝４２．４ｃｍ²／株であり、Ａ₁＝２８７．４ｃｍ²／株であった。この結果を上式（１）に代入することにより、第１指標値Ｉｄとして、次式（１）’のような値が得られた。なお、式（１）のｄは１４（日）である。
Ｉｄ＝（In287.4-ln42.4）/14＝0.136 …（１）’

（第２指標値の算出）
本発明者は、第２指標値ＬＵＥ_1→2を算出するにあたり、乾物重Ｄｍ₁、Ｄｍ₂を上式（２）から求めた（図４のステップＳ２０に対応）。なお、本実施例において、上式（２）の係数ａ、ｂは、ａ＝１／１４７．３、ｂ＝０である。したがって、乾物重Ｄｍ₁、Ｄｍ₂は次式（２）’、（２）”のように求めることができた。
Ｄｍ₁＝(1/147.3)×Ａ₁＝（1/147.3）×42.4＝0.29（g/株） …（２）’
Ｄｍ₂＝(1/147.3)×Ａ₂＝（1/147.3）×287.4＝1.95（g/株） …（２）”

更に、本発明者は、撮影間の総受光量Ｌｉ_1→2を上式（６）から求めた（図４のステップＳ２４に対応）。上式（６）のＬｉ_nは、１回目の撮影日（１日目：ｎ＝１）と２回目の撮影日（１４日目：ｎ＝１４）までの各日（ｎ日目）の日射量（実測値）をＳｒ_nとすると、次式（１１）のように表すことができる。
Ｌｉ_n＝Ａ_n×Ｓｒ_n …（１１）

また、外葉肥大期のキャベツの投影葉面積は指数関数的な増加を示すので、Ａ₁とＡ₂の間の投影葉面積補完値Ａ_nは次式（１２）の関数により求めた。
Ａ_n＝Ａ₁×EXP（（lnＡ₂－lnＡ₁）／ｄ×ｎ） …（１２）

この結果、１回目と２回目の撮影の間の総受光量Ｌｉ_1→2は、図１６に示すように、３．８２（ＭＪ／株）と求められた。

そして、本発明者は、１回目と２回目の撮影の間の総受光量Ｌｉ_1→2と、乾物重の変化（Ｄｍ₂－Ｄｍ₁）を用いて、上式（７）に基づいて、次式（７）’のように１回目と２回目の撮影の間の第２指標値（光利用効率）ＬＵＥ_1→2を算出した。
ＬＵＥ_1→2＝（1.95－0.29）／3.82＝0.435（g/MJ） …（７）’

（閾値Ｔｈ_Id、Ｔｈ_LUEの設定）
本発明者は、気温、土壌水分が抑制要因にはならない条件（光利用効率を減少させない範囲）で、日射量実測値を既存の生育モデルに入力して投影葉面積（Ａ_sim）を推定した。そして、本発明者は、推定した投影葉面積（Ａ_sim）を用いて、Ｔｈ_Id、Ｔｈ_LUEを設定した。すなわち、Ｔｈ_Id、Ｔｈ_LUEはＩｄ、ＬＵＥそれぞれのポテンシャル値であるといえる。具体的には、本発明者は、Ｔｈ_Idについては、上式（１）と同様の方法で、投影葉面積Ａ₁、Ａ_simを用いて算出した。また、Ｔｈ_LUEは、上式（７）と同様の方法で、投影葉面積Ａ₁、Ａ_simから求められる乾物重Ｄｍ₁、Ｄｍ_simと、総受光量Ｌｉ_1→2と、を用いて算出した。なお、図１７の投影葉面積Ａ₁、乾物重Ｄｍ₁は、図１６と同値を用いた。図１７の例では、Ｔｈ_Id＝0.176（cm²/cm²/日）、Ｔｈ_LUE＝0.560（g/MJ）となった。

（生育判定）
本発明者は、上述した１回目の撮影画像と２回目の撮影画像から得られた第１指標値Ｉｄ（＝0.136）及び第２指標値ＬＵＥ_1→2（＝0.435）を、Ｔｈ_Id（＝0.176）及びＴｈ_LUE（＝0.560）と比較した（図４のステップＳ２６に対応）。その結果、Ｉｄ＜Ｔｈ_Id、ＬＵＥ_1→2＜Ｔｈ_LUEとなった。この比較結果は、図９のiiiに対応するため、「生育は抑制されている。日射量は十分であるため、他の抑制要因がある。」と判定（診断）できた。本発明者は、この判定結果を受け、以下のように抑制要因の特定処理（ステップＳ５２～Ｓ５６）を行った。

（抑制要因の特定）
抑制要因の特定に際し、本発明者は、図１８（ａ）、図１８（ｂ）に示すようなスコア化テーブルを準備した。図１８（ａ）のスコア化テーブルは、期間平均気温とスコア（気温スコア）との関係を示すテーブルである。また、図１８（ｂ）のスコア化テーブルは、期間平均土壌水分とスコア（土壌水分スコア）との関係を示すテーブルである。なお、これらのスコア化テーブルは、２０１８年秋、２０１９年春及び秋に栽培試験を実施して得られた値を基に、光利用効率と気温、土壌水分との関係性を解析して、得たものである。

ここで、スコアについては、図１８（ｃ）に示すような環境計測値とスコアとの関係がある場合において、次式（１３）に示すような式で表すことができる。

上式（１３）のＢは図１８（ｃ）に示す最小値、Ａは図１８（ｃ）に示すスコア値の幅（最大値と最小値の差）を意味する。また、ａは図１８（ｃ）に示す天面幅（スコアが最大値を示す環境計測値の幅）、ｃは図１８（ｃ）に示す天面幅の中央値（スコアが最大値を示す環境計測値の中央値）、ｂは図１８（ｃ）に示す傾斜を意味する。これらの値は、複数の栽培試験で得られた日射利用効率と、日射利用効率を試算した２つのサンプリング日間の環境計測値の平均値を式（１３）に代入して得られる値との残渣平方和が最小となるように式（１３）の係数を推定し、さらにスコアの最小値を０、最大値を１にスケーリングすることで、求めることができる。なお、図１８（ａ）、図１８（ｂ）のスコア化テーブルの各値ａ～ｃ、Ａ、Ｂは、図１８（ｄ）に示す通りである。

上述した２回の画像取得期間内（４月１３日～４月２７日）における気温と土壌水分（土壌含水率）の期間平均値は図１９（ａ）に示す通りであった。この場合、図１８（ａ）、図１８（ｂ）のスコア化テーブル又は上式（１３）を用いて気温スコアＳｃ_T、土壌水分スコアＳｃ_Wを求めた結果、図１９（ｂ）に示す通りとなった。

本発明者は、上記のように求められた気温スコアＳｃ_T、土壌水分スコアＳｃ_Wを用いて、生育抑制の要因が、気温や土壌水分にあるのかを判断した（図５のステップＳ５６に対応）。なお、本実施例における生育抑制要因の特定方法は、上記実施形態とは異なるが、以下のような方法で行うこととした。まず、本発明者は、図２０（ａ）に示すように、気温スコアＳｃ_Tと土壌水分スコアＳｃ_Wの積を求め、環境スコアＳとした。また、本発明者は、ＬＵＥ_1→2／Ｓを求め、Ｍとしし、更に、Ｍを閾値Ｔｈ_LUEで除した値（Ｍ／Ｔｈ_LUE）を求めた。そして、本発明者は、Ｍ／Ｔｈ_LUEの値が、図２０（ｂ）に示すいずれの条件を満たすかに基づいて、生育抑制の要因を判定した。

図２０（ａ）の場合、図２０（ｂ）の条件１を満たすため、「生育抑制は主にスコア化した環境要因による。」と判定された。

更に、本発明者は、図２０（ｃ）に示すように、各環境要因の影響度（＝１－環境スコア）を試算した。その結果、気温の影響度が０．１３３、土壌水分の影響度が０であることが分かった。

以上のように、本実施例によれば、対象作物がキャベツの場合において、上記実施形態で説明した生育診断や生育抑制要因の特定処理を適切に実施しうることが分かった。

なお、上記の処理機能は、コンピュータによって実現することができる。その場合、処理装置が有すべき機能の処理内容を記述したプログラムが提供される。そのプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。処理内容を記述したプログラムは、コンピュータで読み取り可能な記憶媒体（ただし、搬送波は除く）に記録しておくことができる。

プログラムを流通させる場合には、例えば、そのプログラムが記録されたＤＶＤ（Digital Versatile Disc）、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）などの可搬型記憶媒体の形態で販売される。また、プログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することもできる。

プログラムを実行するコンピュータは、例えば、可搬型記憶媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、自己の記憶装置に格納する。そして、コンピュータは、自己の記憶装置からプログラムを読み取り、プログラムに従った処理を実行する。なお、コンピュータは、可搬型記憶媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することもできる。また、コンピュータは、サーバコンピュータからプログラムが転送されるごとに、逐次、受け取ったプログラムに従った処理を実行することもできる。

上述した実施形態は本発明の好適な実施の例である。但し、これに限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々変形実施可能である。

１０ ユーザ端末
３２ 第１指標値算出部
３６ 第２指標値算出部
３８ 生育状態判定部
４０ 抑制要因特定部
４２ 出力部
７０ カメラ
７２ センサ

Claims (8)

1. 同一地点の作物に対して実施した２回の生体調査の結果に基づいて、前記作物の生育速度を示す第１指標値を算出し、
   前記２回の生体調査の結果と、前記生体調査の間の期間における日射に関する情報と、に基づいて、日射による生育への影響を示す第２指標値を算出し、
   前記第１指標値と前記第２指標値とに基づいて、前記作物の生育が順調か否か及び前記作物の生育と日射との関係を判定し、
   前記判定する処理の判定結果を出力する、
   処理をコンピュータが実行することを特徴とする作物の生育診断方法。
2. 前記判定する処理では、前記第１指標値と予め定めた第１閾値との大小関係と、前記第２指標値と予め定めた第２閾値との大小関係と、に基づいて、前記作物の生育が順調か否か及び前記作物の生育と日射との関係を判定する、ことを特徴とする請求項１に記載の作物の生育診断方法。
3. 前記判定する処理では、
   前記第２指標値が前記第２閾値よりも小さい場合に、前記２回の生体調査の間の期間における複数種類の環境要素に関する情報を取得し、
   取得した前記複数種類の環境要素に関する情報に基づいて、前記作物の生育に対して影響を与えている環境要素の種類を特定する、
   ことを特徴とする請求項２に記載の作物の生育診断方法。
4. 前記生体調査では、前記作物を所定高さから撮影した画像から前記作物の葉に関する情報を取得する、ことを特徴とする請求項１～３のいずれか一項に記載の作物の生育診断方法。
5. 前記生体調査では、前記作物の投影葉面積を取得し、
   前記第１指標値は、前記２回の生体調査で得られた前記作物の投影葉面積を用いて算出した生育速度を示す値であり、
   前記第２指標値は、前記２回の生体調査で得られた前記作物の投影葉面積を用いて算出した、前記２回の生体調査の間における前記作物の光利用効率である、ことを特徴とする請求項１～４のいずれか一項に記載の作物の生育診断方法。
6. 前記出力する処理では、所定範囲内の複数の作物に関する前記判定する処理の結果を、前記複数の作物の位置と関連付けて表示する画面を生成し、出力する、ことを特徴とする請求項１～５のいずれか一項に記載の作物の生育診断方法。
7. 前記出力する処理では、前記判定する処理の結果に基づいて、前記作物が生育する環境を制御する環境制御装置に対して制御指示を出力する、ことを特徴とする請求項１～６のいずれか一項に記載の作物の生育診断方法。
8. 同一地点の作物に対して実施した２回の生体調査の結果に基づいて、前記作物の生育速度を示す第１指標値を算出し、
   前記２回の生体調査の結果と、前記生体調査の間の期間における日射に関する情報と、に基づいて、日射による生育への影響を示す第２指標値を算出し、
   前記第１指標値と前記第２指標値とに基づいて、前記作物の生育が順調か否か及び前記作物の生育と日射との関係を判定し、
   前記判定する処理の判定結果を出力する、
   処理をコンピュータに実行させることを特徴とする作物の生育診断プログラム。

Images