原著論文
大規模均一栽培圃場における生育調査のドローン空撮による代替手法の検討
2023 年 2023 巻 15 号 p. 53-60
詳細
2023 年 2023 巻 15 号 p. 53-60
大規模均一栽培圃場での生育調査は，大きな労力を要する．それらをドローン空撮および画像解析により代替する手法を検討した．複数の画像から対象物の形状を３次元的に構成する手法に SfM がある．ドローン空撮画像を SfM 処理する手法は広く使われているが，空撮の際に地上基準点（GCP）を適切に配置することでその精度を高める必要がある．大規模均一栽培圃場では全面を作物で被覆しているため，GCP の設置に制約がある．特に圃場内部への GCP の設置は作業強度が高く，かつ圃場均一性への影響もあるため極力避けたい．検証の結果，圃場の 4 つの角および 4 つの辺の中程の計 8 箇所に GCP を設置すれば，圃場内部にも GCP を設置した場合と比較しても遜色ない精度で草高を推定できることが分かった．また，ドローンの飛行時間は飛行高度の 2 乗に反比例するため，より高い高度で飛行することが時間短縮につながるが，高度を上げればそれだけ空間分解能は低下する．本試験の条件下では，4 ha 程度の圃場の場合，前述の GCP を 8 個設置し飛行高度を 50 m とすることで十分な精度が得られた．これらにより，ドローン空撮により従来の生育調査をより省力的にかつ高精度に行うことができることを示した．
試験研究圃場での精密な試験の前提条件として，地力を維持し，かつその均一度を向上させることは極めて重要である．作物の均一栽培および青刈鋤込みにより地力の向上および均一化を図りつつ，草丈の測定により不均一度を把握し圃場設計に活用する手法が広く用いられている（遠山 1988）．草丈の測定は，圃場内をくまなく巡回し，設定した調査地点にて手動で計測を行う方法が一般的である．この作業は，特に大規模圃場では大きな労力を要する．また，できる限り多くの調査地点で測定を行うことが望ましいが，調査地点の数に比例して労力も増加するため，それらには限りがある．
現在，ドローン空撮により対象物の３次元的構造を把握する技術が様々な分野で用いられている．この手法は，Structure from Motion (SfM)と呼ばれ(Ullman 1979)，農業分野においても，草高やバイオマスの推定などに応用する事例が数多く存在する(Holman et al. 2016, Gilliot et al. 2021, 藤原ら 2022)．均一栽培圃場においても同様の手法を活用することで，より省力的にかつ精密に圃場の状態を取得することが期待できる．ドローン空撮画像に対する SfM 処理において，カメラのレンズ特性などの影響から生成された画像にドーム状の歪みが発生することが知られている（Lee and Han 2020, Zhou et al. 2020）．地上に複数の対空標識を適切な配置で設置し，それらを地上基準点(Ground Control Points, 以下 GCP)として用いることで，この歪みは補正することができる．しかしながら，均一栽培圃場では全面が作物により被覆されていることが一般的であるため，それらの設置場所には大きな制約がある．作付範囲の外側への GCP の設置は，作物の干渉を受けることがないため，また，車両等を用いた移動や機材の運搬も可能であるため，作業負担はそこまで大きくはない．しかし，圃場内部に GCP を設置する際は，作物を分け入って GCP を運搬し，設置箇所の作物を除去するなどの作業が必要なことから作業強度が高い．また，作物を除去することによる圃場の均一性への影響もあることから，極力内部への GCP の設置は回避しつつ適切な数を確保する必要がある．空撮をする際に，カメラを通常の真下向きの他に角度を付けて撮影することで，少数の GCP でも SfM の歪みが軽減されることが報告されている（Nesbit and Hugenholtz 2019）．しかし，角度を付けた撮影は，複数の方向からの撮影が必要になるため，大規模圃場においては労力の問題から必ずしも望ましい方法ではない．
飛行高度が低いほど得られる画像の地上分解能は向上するが，飛行時間は飛行高度の 2 乗に反比例するため，地上分解能と飛行に係る労力とのトレードオフを考慮した上で飛行計画を策定するための一定の指針が必要である．均一栽培圃場における調査の主な目的は，生育のばらつきを把握することであるため，例えば育種研究圃場のようなより精緻なデータの取得が求められる場面とは別の視点で条件設定を行うことで，省力的に大規模圃場に対応するといった工夫も求められる．本報告では，GCP の配置方法および飛行高度が SfM の結果に与える影響について検証し，大規模均一栽培圃場に適した空撮方法の提示を行う．
１．調査対象圃場と使用ドローン
調査対象として，農研機構北海道農業研究センター芽室研究拠点（北海道河西郡）内の均一栽培試験圃場を用いた．作物はダイズ（品種：ユキホマレ），作付面積は約 158 × 258 m（約 4.1 ha），2022 年 5 月 26 日に畦幅 66 cm，株間 18 cm にて播種，8 月 15 日にチョッパーにより細断，8 月 17 日にサブソイラーで心土破砕，8 月 18 日にディスクハローにより鋤込みおよび粗整地を行った後，8 月 22 日にプラウにて完全に整地した．従来手法による生育調査は，8 月 10 日に実施した（図 1）．生育調査は，生育のばらつきを把握することを目的とし，圃場全体にまんべんなく設定された複数の地点の草高を手動で計測することにより実施する．
ドローンは，DJI 社製の Phantom 4 Pro を用いた．飛行高度は，30 m，50m，100 m の 3 通り，カメラアングルは常に 90 度（鉛直真下方向）とし，圃場長辺方向に並行に飛行，画像間のオーバーラップがフロント，サイド共に 80% となるよう飛行経路，速度を調整して自動飛行させた．飛行高度 30，50，100 m で撮影した際に得られる画像は，撮影枚数がそれぞれ約 930，370，130 枚，地上分解能がそれぞれ約 0.79，1.31，2.74 cm/pixel である．空撮は，作物生育時期である生育調査の前日の 8 月 9 日と当日の 8 月 10 日の 2 日間，生育終了後である粗整地実施直後の 8 月 18 日と整地終了後の 9 月 7 日の 2 日間の計 4 日間実施した．
２．GCP の設置と測位
GCP は，圃場の 4 つの角に 1 つずつ，4 つの辺の中程に 1 つずつ，圃場の中心 1 つ，4 つの角と圃場中心を結ぶ直線上の中程に 1 つずつの計 13 個を設置した（図 2 のパターン a）．角と辺の 8 つについては，作物から約 1 m 離した位置に設置した．圃場内部の残りの 5 つについては，作物を 2 畦にわたり約 2 m 抜き取りを行った後，作物が除去された範囲の中心付近に設置した．設置した各 GCP の位置は GNSS レシーバー DG-PRO1RWS（BizStation Corp.）を用いて，RTK 精度（絶対位置の誤差は水平方向，鉛直方向ともに± 0.02 m 程度）で測位した．
３．オルソモザイク画像と数値表層モデル(DSM)画像の取得
オルソモザイク画像および数値表層モデル(Digital Surface Model，以下 DSM)画像は，市販の SfM ソフトウェアである Agisoft Metashape Professional, Ver. 1.7.3 により生成した．手順は以下の通りである．まず，空撮画像のアラインメントを行う．次に，GCP の測位データを読み込み，画像内の各 GCP との紐づけを行う．なお，本ソフトウェアは GCP を自動認識する機能を有しているが本試験では使用していない．その後，高密度クラウドを品質 High（日本語版表記：高）に設定して構築した後，オルソモザイク画像と DSM 画像をエクスポートする．これらの作業を 3 つの飛行高度それぞれに対して実施した．なお，飛行高度 100 m については高密度クラウドの品質を Ultrahigh（日本語版表記：最高）でも生成した（以下，本文や図表等において_uh と表記）．なお，高密度クラウドの生成においては，SfM に加えて MVS (Multi-View Stereo)処理が行われる（以下，MVS 処理を含む場合も SfM と表記）．
また，GCP の配置パターンが結果に与える影響を調べるために，使用する GCP の数を変化させて上記の処理を実施した．用いた配置パターンは，すべての GCP を処理に用いる配置（パターン a），角 4 つ辺 4 つ中心 1 つの配置（パターン sc），角と辺への配置（パターン s），長辺のみの配置（パターン ls），短辺のみの配置（パターン ss），角と中心への配置（パターンcc），角のみへの配置（パターン c）の 7 パターン（図 2）である．以下，飛行高度と GCP 配置パターンの組み合わせを用いて，例えば飛行高度 30 m のパターン a によるデータセットを 30 m-a と表記する．SfM の処理を経て得られる DSM に現れる歪み（以下，SfM の歪みと表記）が上記の各パターンでどの程度発生するかを調べるために，パターン a を除く 6 つのパターンについて，処理に用いられなかった GCP（以下，不使用 GCP）の位置が SfM によりどのように推定されたか調べることで，各配置パターンの SfM の歪みへの影響を検証した．水平方向の推定位置は，オルソモザイク画像上の GCP の座標から取得，垂直方向の推定位置は同地点の DSM 画像の値から取得し，実測値である GCP 測位データとの差異を検証した．
上記の手順を，4 つの空撮日に対し実施した．以下，各日付に対する表記をそれぞれ 20220809，20220810，20220818，20220907 とする．
１．GCP配置パターンの SfM の歪みへの影響
GCP 配置の各パターンにおける不使用 GCP の SfM 推定位置と測位データとの差異を，MAE（平均絶対値誤差）を用いて水平方向と垂直方向それぞれ算出した（図 3）．なお，パターン c については垂直方向の差異が突出して大きかったため，他のパターンの違いを明確にするために図から除外し，別途表 1 に示した．
水平方向の差異については，飛行高度が高いほど差異がやや大きくなる傾向が見られるが，同一飛行高度では GCP 配置パターン間で大きな違いは見られない．飛行高度 100 m に関しては，品質を Ultrahigh とした方が High と比較して差異が小さい．4 ha 規模の圃場における生育のばらつきを把握するという目的に照らせば，また，RTK 測位の精度自体が ±0.02 m 程度であることも考慮すれば，いずれの飛行高度，パターンにおける差異も許容範囲内だと考えられる．
一方，垂直方向の差異はパターン間で顕著な違いが見られた．前述の通り，パターン c はいずれの高度においても差異が他のパターンに比較して特に大きく，このパターンは SfM の歪みが著しいことが確認された．なお，パターン c については，生育調査を代替するという本試験の目的に合致しないことが明確であるため，以後考察の対象とはしない．飛行高度 100 m における品質 High の差異が特に 20220809 および 20220810 において大きくなる傾向が見られた．しかし，同じく飛行高度 100 m の品質 Ultrahigh についてはいずれのパターンも品質 High に比べて差異は低く抑えられている．品質を Ultrahigh にすることにより同じ GCP 配置であっても結果が大きく改善することから，品質 High で MAE が大きくなる要因は SfM の歪みというより，DSM の解像度が不足していることによることだと推測された．20220809 および 20220810 における内部の GCP は，周囲の作物を除去した影響により草高が大きく変化する場所にあり，それらの立体構造を品質 High では反映しきれていない可能性が考えられる．飛行高度 100 m，品質 High 以外の 3 種類のデータについて，パターン ss は他のパターンより差異がかなり大きく，パターン ss は SfM を大きく歪ませることが確認できた．パターン ls は，飛行高度 30 m および 50 m で大きな歪みを生じさせていることが分かる．パターン sc，パターン s，パターン cc の 3 つのパターンについては，いずれも MAE が 0～0.04 m 程度に抑えられており，この程度の差異であれば圃場全体の生育のばらつきを把握するという目的に照らせば許容範囲内だと考えられる．
２．1 メートルメッシュによる評価
最も飛行高度が低く GCP の数が多い 30 m-a のデータがより真値に近いと仮定し，そこから GCP を削減すること，また，飛行高度を上げることで DSM にどの程度の違いが現れるかを検証した．飛行高度ごとに解像度が異なる DSM を比較するために，作付範囲を 158 × 258 の 40,764 個の 1 メートルメッシュに分割し，楕円体高について各メッシュでの平均値（以下，1 メートルメッシュ高）を算出し，これを各パターンでの比較の基準とした．図 4（a）は，30 m-aと各データセットの間での 1 メートルメッシュ高の二乗平均平方根誤差（RMSE）を示す．また，草高に相当する値を算出するために作物生育期の 20220809 と 20220810 と裸地の 20220818 と 20220907 の間で組み合わせを変え，各組合せにおいて 1 メートルメッシュ高の差分を算出した．図 4（b）は，その差分を基準とした 30 m-a と他のデータセットの間での RMSE を算出した結果である．
図 4 において，(a)，(b)共に飛行高度 100 m の品質 High と Ultrahigh との差異はほとんど見られない．つまり，不使用 GCP の検証において現れていた局所的な誤差は，1 メートルメッシュで平均化することにより表出しなくなると考えられる．
すべての飛行高度において，RMSE はパターン ss で最も高く，ついでパターン ls が高い．また，それ以外のパターンでは低く抑えられている．これらは，前述の不使用 GCP と DSM の垂直方向の差異とも傾向が同じであり，SfM の歪みが 1 メートルメッシュ高全体にも同様の影響を及ぼしていることが分かる．
３．適切な飛行高度と GCP の配置
本試験の飛行高度 30 m，50 m，100 m のフライトにおいて離陸から着陸までに要した時間は，それぞれ 37 分（途中 1 回のバッテリー交換の時間を含む），14 分，7 分程度であった．50 m のフライト時間は 30 m のフライト時間の 4 割未満でありバッテリー交換も不要なことから，飛行高度 50 m はより省力的である．図 4（a）に示す通り，50 m-a，50 m-sc，50 m-s の 30 m-a と比較した１メートルメッシュ高の RMSE は 0.02 程度であり，これらは生育のばらつきを把握するという観点からも十分許容できる．また，留意すべき点として，均一栽培圃場や裸地の空撮により得られる画像は，画像間の類似度が高く特徴量の乏しいものになりやすいため，画像のアラインメントが上手くいかないなどの理由で SfM が失敗する可能性があり，それは高度が低い場合ほど発生しやすい．このことからも，また飛行高度 30 m と 50 m では飛行時間が大きく異なるが 1 メートルメッシュ高の差は小さいことも考慮すれば，敢えて飛行高度 30 m を選択する必要性は低いと考えられる．また，飛行高度 100 m においても，一部の GCP 配置を除いて 1 メートルメッシュ値の RMSE は十分許容できる範囲内だと考えられるが，フライト時間は 50 m と比較して半分程度になるがその差は 7 分程度であるため，データの精緻さ以上に一層の作業の省力化を求める場合以外は敢えて選択する必要はない．
飛行高度を 50 m とした際の GCP の配置について考える．圃場内部の GCP 設置は，作業強度が高く作物の均一性への影響があることから，周辺への設置と比較して優先度は下がる．圃場内部への GCP 設置のないパターンは，50 m-s，50 m-ls，50 m-ss の 3 つであるが，50 m-s 以外の 2 つのパターンは誤差が大きくなることから選択肢になり得ない．50 m-sは，50 m-a，50 m-sc と比較してもほとんど精度に差は見られないため，圃場内部への GCP 設置が必要な 50 m-a，50 m-sc を敢えて選択する理由はない．また，50 m-cc は 50 m-s と比較して GCP の数を 3 つ少なくできるが，内部に 1 つの GCP を設置するよりも辺に 4 つの GCP を設置する方が作業強度が低いこと，更に 50 m-s の方が精度が高いことからも，50 m-cc も選択肢にはなり得ない．これらのことは，飛行高度 100 m についてもほぼ同様に考えられる．
各パターンについて，さらに 2 つの比較データを示す．表 2 は，図 4（b）で用いたデータについて，30 m-a とその他の各パターンの 1 メートルメッシュ値の相関係数を示している．50 m-s は，すべての組み合わせにおいて相関係数が 0.97 を超えており，データの傾向が 30 m-a と比較してほとんど相違ないことが分かる．100 m-s_uh も50 m-s と比較するとやや劣るものの，いずれも相関係数が 0.95 を超えておりデータの傾向の把握には十分に機能するであろう．図 5 は，図 4（b）と同様のデータについて 1 メートルメッシュ値に基いて着色，可視化したものである．各画像は作付範囲における草高の平均値を表しており，値が高い地点ほど薄く描画されている．ここでも，パターン s が 30 m-a とほとんど相違なくばらつきの把握が行えることが視覚的に確認できる．
４．従来調査との比較
手動による生育調査は，作付範囲を約 15 m 間隔の格子に区切り，計 126 箇所の格子点において草丈（鉛直方向に作物を最大限引き伸ばした高さ）を計測した．各調査地点において，中心の 1 点とそこから四方に概ね 1.4 m 離れた 4 点の合計 5 点の作物を計測し 5 点の平均を結果とした．図 6 に，生育調査による実測草丈と DSM 算出値の関係を示す．DSM 算出値は，生育調査日である 20220810 と完全に整地が完了した後の 20220907 の 2 つの１メートルメッシュ高の差分について，生育調査地点を含む 3 × 3 m のメッシュ値を抽出しそれらの平均値を用いた．生育調査が作物を引き延ばした状態の高さを計測しているのに対して，DSM 算出値は作物の地面からの高さの平均を取っているため，前者の値の方が大きく出ているが傾向はよく一致している．つまり，DSM 算出値は，従来調査の代替になり得ると言える．なお，DSM 算出値は実測草丈より過小評価されるが，調査の目的は圃場全体のばらつきの把握であり個々の草丈の計測ではないため，特に従来手法の測定値への換算をせずとも DSM 算出値の提示のみで目的は達成できる．
生育調査に要する所要時間は概ね以下の通りである．調査地点を示す支柱（目印）の設置，調査終了後のそれらの撤収作業には，2 人組で 90 分を要した．草丈の計測は，3 人組で 180 分を要した．また，調査終了後の野帳から電子媒体への転記，統計処理，マップへの落とし込みに 1 人で 120 分を要した．以上の合計は，840（人・分）である．これらは一方，ドローン空撮および画像解析に要する所要時間は 50 m-sのパターンの場合で概ね以下の通りである．GCP 8 箇所の設置および測位，空撮終了後のそれらの撤収作業には，2 人組で 40 分を要した．ドローン空撮には，飛行前後の準備も含め 2 人組で 20 分を要した．SfM 処理にかかるファイルの整理，GCP データの生成と GCP 位置の確認に 1 人で 50 分を要した．その他 SfM ソフトウェアの操作時間として 1 人で 10 分を要した（ただし，人手を介さずにソフトウェアが自動計算を行うなどの過程の時間は考慮していない）．上記の作業は作物生育時期と裸地期の 2 回行う必要があるため，以上を合計しさらに 2 倍した結果は 360（人・分）である．なお，本試験ではこれらの処理に加えて DSM 算出値を検証するために1メートルメッシュごとの数値データやその可視化画像（図 5）の生成に自作のソフトウェアを用いたが，単に 2 つの異なる DSM ファイルの差分を可視化するだけであれば，例えばオープンソフトウェアである QGIS の標準機能で実現できる．その作業時間は，ソフトウェアの操作を習得していれば数分で完了する．
本試験で用いた機材により 4 ha 程度の圃場を空撮する場合，作業時間や精度を総合的に考えると 50 m-s が空撮に適した方法だと結論付ける．100 m-s_uh は，やや精度が劣るものの結果は比較的良好であるため，特に作業時間を重視する際には選択肢となり得る．
また，ドローン空撮の利用により従来の生育調査の作業時間を大きく削減できることが分かった．それ以上に，人手による調査では調査地点の場所や数を限定せざるを得ないが，空撮画像の解析は全体を面として捉えることができることの意義は大きい．
すべての著者は開示すべき利益相反はない．
本研究の一部は，AIP 加速課題“ビッグデータ駆動型 AI 農業創出のための CPS 基盤の研究”（JPMJCR21U3）の支援により実施された．
均一栽培圃場における調査に当たって，栽培管理および生育調査を統括していただくと共に，圃場内 GCP 設置等に関して貴重な助言をいただいた北海道技術支援センター北海道第 3 業務科の小川英明総括作業長に感謝の意を表する．
