Ⅰ 緒 言 中山間地の傾斜地水田の圃場整備では，大規模な 地形改変を伴う基盤造整備が行われることが多い． 一筆圃場内の不均一な生育障害は圃場整備前の地形 と密接に関係していることがある．例えば下層土に 起因する地力ムラや，切土・盛土部での土壌の変化， かつて水路であったところの地盤沈下，浅礫層の分 布，盛土部での湿害の発生など2, 8, 9, 11, 13），過去の地 形に起因する圃場内の土壌環境の不均一性に関する 報告がされており，基盤造成を行った圃場で生じた 不均一な生育不良を考える上で，基盤造整備された 圃場の過去の地形を考慮することは重要である． しかし，過去の地形を復元するためには圃場整備 前の地形図を入手する必要がある．国土地理院発行 の 2．5 万分の 1 の旧版地形図は全国で容易に入手 できるが，等高線間隔が 20 m でかつ圃場スケール を対象とするには解像度が足りない．圃場整備時に 作成された圃場整備計画図は，1,000 分の 1 程度の 縮尺で圃場整備前と整備後の計画地形の両方が記載 されており，地形改変箇所を復元するために必要な 情報が記載されている（第 1 図）．しかし，これら の地図は数十年以上も前に発行されたものが多く， 土地改良区等で保管されるが散逸していることもあ り，入手できても不鮮明な紙地図であることが多く デジタル標高モデル（Digital Elevation Model; DEM） の作成は困難である．このため，圃場整備計画図が 入手できないことを想定し，また DEM の構築を容 易にするため，別の汎用性のあるデータを用いた復 2019 年 7 月 31 日受領 2019 年 10 月 10 日受理 Correspondence: shimizuy@affrc.go.jp 〔原著論文〕
SfM 多視点ステレオ写真測量による過去の空中写真からの
三次元地形モデルの構築
清水裕太
＊，＊＊・松森堅治
＊ ＊国立研究開発法人 農業・食品産業技術総合研究機構 西日本農業研究センター ＊＊国立研究開発法人 農業・食品産業技術総合研究機構 農業情報研究センターReconstruction of digital elevation model based on historical aerial photos
using SfM-MVS technique
Yuta Shimizu
＊, ＊＊and Kenji Matsumori
＊＊Western Region Agricultural Research Center, NARO ＊＊Research Center for Agricultural Information Technology, NARO
第 1 図 圃場整備計画図の例．
圃場整備前の地形および測量された標高点と圃場整備後の圃場区 画および計画標高が記載されている．
い．そこで，本研究では 1970 年代の空中写真と SfM-MVS ソフトウェアを使用して，圃場整備前の 地形を集落スケールで構築し，その再現精度を検証 することを目的とした． Ⅱ 材料および方法 本研究では，広島県東広島市の中山間地に位置す る農業集落を対象とした（第 2 図）．この地域では 国土地理院によって 1970 年代に空中写真が撮影さ れており，1974 年に撮影された 4 枚を使用した（第 1 表）．写真はツァイス社製 RMKA カメラによって 高度 2,000 m 程度から撮影されたもので，その後 1,200 dpi の解像度でデジタル化され画像データとし て提供されている．地理情報と重ね合わせるために は，幾何補正された合成写真（オルソモザイク）に 地理座標を付与する必要がある．通常は，地上の目 印となる場所（地上標定点）の既知の緯度経度およ び標高をあらかじめ測量等によって調べておくが， 過去の地形を復元する場合は，既に地形が変化して いる場所では測量による位置情報を取得することが 困難である． 本研究では過去と現在の空中写真に共通して写っ ている比較的頑強な構造物の標高は変化しないと仮 定して（第 3 図），空中写真から 34 地点を選定し， 17 地点を地上標定点に，残りの 17 地点を検証点と 元手法を確立することが必要である．国土地理院に よって航空機から撮影された空中写真は 1960 年代 から入手可能であり，その中でも 1974 年から国土 情報整備事業のために撮影された空中写真は他の年 代よりも高い解像度で日本全国を撮影しており10）， これらはインターネットを通して撮影したデータが 入手可能であることから有用である．空中写真をス テレオペア画像として実体視する三次元化は従来か ら行われてきたが，近年ではドローンに代表される UAV の著しい進歩と共に，写真測量用のソフトウェ アも格段に進歩しており，従来の手法よりも簡便に 三次元モデルを作成することができるようになっ た．このこともあり，UAV を用いた精細な地形モデ ルの作成や精度の検討などは近年行われ4, 12），空中 写真を大量に撮影すれば，ほぼ半自動で幾何補正か らモザイク合成まで行える多視点ステレオ写真測量 （Structure from Motion/Multi View Stereo; SfM-MVS）
ソフトウェアを用いると，比較的容易に三次元モデ ルを作成することができる．この手法を用いた地形 に関わる先行研究としては，地すべりモニタリング など詳細な地形の再現だけでなく，時系列的な変化 を検出や5），棚田を対象とした研究事例などが報告 されている3, 7）．過去の写真から三次元モデルを復 元した事例として，植生や地形変化のモニタリング を行った研究1）などが報告されているが，耕区から 農区程度の比較的小規模なスケールでの検証はな 第 2 図 地域概要図． 黒丸点は地上標定点を，三角点は検証点を，黒枠は再現性検証地区を，灰色のエリアは写真の撮影 範囲をそれぞれ表す．等高線の間隔は 10m で，基盤地図標高 10m メッシュより作成した．
Align Photos 機能による写真の撮影位置の推定後に， Build Dense cloud 機能によって三次元のポイントク ラウド作成の手順で行った．Align Photos は，GCP 情報がある場合は，その情報を基に写真を初期配置 し，その後に各写真間に共通して写っている特徴点 の位置的関係から最終的な撮影位置が決定される． ポイントクラウド作成に用いた各種パラメータは第 2 表に示す．作成されたポイントクラウドには，ノ イズが含まれており，Agisoft Photoscan 自体でも除 去は可能であるが，集落スケールでノイズをほぼ手 動で除去することは現実的では無いため，地理空間 解析ソフトウェア ArcGIS 10.5 および点群データ解 析ソフトウェア Cloud Compare version2.10.2（https:// www.cloudcompare.org）を使用してノイズの除去を 行った．圃場整備計画図に記載された圃場整備前の 等高線を読み取り（第 1 図），ArcGIS を用いてデジ タル標高モデルを作成した． した（第 2 図）．これらの地点において，ネットワー ク型 RTK 測量機材を用いて計測した座標および標 高点を真値と仮定して解析に使用した．本研究では 写真解析を行う Agisoft PhotoScan version 1.4.3 build 6529 を使用した．手順の概要は，写真の読み込み， 複数の写真に共通する特徴点（タイポイントあるい はパスポイント）を割り出して撮影時のカメラの位 置や向きを推定してから，ステレオペア画像として 空中三角測量による三次元的な密な点群（ポイント クラウド）を生成し，それらをもとに DEM を作成 するという流れである．入力前の前処理として，写 真のフレームの除去を行い，写真中心点の座標およ び撮影高度を入力した．中心点の座標は地理院地図 の過去のオルソモザイク画像地図から世界測地系の 座標を割り出した．撮影高度は写真に写されている 高度計の値を読み取り入力した．撮影時のカメラの 三次元的な姿勢を表す Yaw，Pitch，Roll については 写真に写されている水準器から，ほぼ水平であると 仮定して入力した．これらの数値はフィッティング パラメータの初期値として与えているが，最終的に は写真と地上基準点との位置的な関係が最もフィッ トするように調整される．Agisoft Photoscan では 第 3 図 地上標定点の例． （a）1974 年撮影，（b）2010 年撮影．α - α’およびβ - β’は同 じ位置である． 機能 パラメータ 採用した値
Camera Calibration Pixcel size (mm) 0.020 151.7 Align Photos Accuracy Ultra High s e Y ri a P
Key point limit 40,000 Tie point limit 4,000 Aligned Photos 4 Mean photo quality 1.2 Build Dense Cloud Quality Ultra High Depth filtering Aggressive
第 2 表 Photoscan 解析時のパラメータ設定． コース 写真 撮影高度* 撮影 カメラ 焦点距離 番号 番号 (m) 縮尺 名称 (mm) 東距(m) 北距(m) CCG748 C14A 70 1975/3/4 1,950 10,000 RMKA 153.29 63077.9 -165332.9 CCG748 C14A 71 1975/3/4 1,950 10,000 RMKA 153.29 63950.9 -165278.6 CCG748 C14A 72 1975/3/4 1,950 10,000 RMKA 153.29 64804.5 -165261.7 CCG748 C14A 73 1975/3/4 1,950 10,000 RMKA 153.29 65717.1 -165243.0 中心座標** 整理番号 撮影年月日 第 1 表 解析に使用した空中写真． ＊ 写真内の高度計の値を使用． ＊＊ 地理院地図より中心座標を求めた．平面直角座標系第Ⅲ系による投影座標系．
画像しかないため，限られた情報から除去すること が必要である．圃場のみを対象とする目的の場合に は，これらの比較的大きな地物の影響は小さいと考 えられるが，棚田と山林の境界や，畦畔に樹木が植 えられている場合は樹枝が圃場内に迫り出している こともあり，この場合はこれらの影響を除外する必 要がある．撮影時期によっては圃場内に背丈の高い 作物あるいは一時的な構造物などが映り込んでいる 場合があるが，それらの影響が比較的少ない写真を 使用することが推奨される． 2 地表面ノイズ処理の必要性 今回の主要な解析の対象としている棚田において は，撮影時期が冬季であったため大部分の圃場には 植生はほとんど見られなかったが，田面に稲わらが 積み上げられているようなところが確認された（第 5 図 a）．解析結果を拡大すると地表面に細かなノイ ズが多数出現しており，稲わら等の地物の場所とは 異なる場所においてもノイズが出現していたことか ら（第 5 図 b），地物の影響というよりは写真自体の ノイズあるいは三次元の点群を作成する際の欠落点 に起因するものと考えられた．このため，圃場の平 均的な地表面の標高を再現する目的の場合はそれら 地物の影響を無視することができるが，前処理とし てノイズの除去が必要であると考えられた．小さな 距離測定誤差によって引き起こされるデータの不規 則性は，統計分析で取り除くことが難しいが，追加 の写真を取得することが不可能である状況において も解決する方法を考える必要がある．そこで，ノイ ズ除去のために，等高線の長さによるフィルタリン グおよび移動最小二乗平滑化の 2 通りの処理を試み た．まず，等高線の長さによるフィルタリングにつ いて，構築された DSM をポリラインの等高線とし （第 5 図 b），ArcGIS を用いて等高線の各線分の長さ をジオメトリ演算によって求め，閾値以下の長さを 持つ等高線を除去した．閾値を 100 m にしてフィル タリングを行った結果，棚田境界部の等高線が鋸状 をしているところも見られるが，全体的に細かなノ イズが除去されており，棚田の地形を再現していた （第 5 図 c）．次に，Cloud Compare の移動最小二乗近 似平滑化・再構築（Moving Least Squares Smoothing and Reconstruction）ツールを使用してノイズの除去 Ⅲ 結果と考察 1 地表表層モデルの構築 解析の結果，空中写真から圃場レベルの精度を持 つ集落スケールでの三次元地形を構築した（第 4 図 a）．1,200 dpi でデジタル化された縮尺 1：10,000 の 空中写真 4 枚から構築できた表層モデル（Digital Surface Model: DSM）の解像度は 1 ピクセル当たり 44 cm で，圃場スケールの地形を表現するための解 像度を持つことが確認された．写真周辺部ではオー バーラップ率が低く，空中写真から判別される地上 の特徴点が不足したため，ぼやけてうまく構築でき なかったエリアが一部で見られた．これらのエリア を除いて，構築できた DSM の面積は約 4 km2であっ た．解析後に出力される三次元データは厳密には地 上の標高モデル（DEM）ではなく，オルソモザイク 画像（第 4 図 b）からも判読できるように樹木や家 屋等の地物を含んだ表層モデル（DSM）となってい る．このため，地表面の標高に変換するためには地 物の影響を考慮して除去することが必要となる．こ れまで，地上の樹木を除去する手法が考案され，例 えば可視画像と近赤外画像を組み合わせた検出方法 なども提案されているが6），過去の空中写真は可視 第 4 図 空中写真より復元された 1974 年当時の三次元地 形（a）およびオルソモザイク画像（b）．
以上より，空中写真をそのまま解析に使用した場 合，地表面に多くのノイズが出現するため，ノイズ を低減させる前処理が必要であることが示された． 本研究では 2 通りのノイズ除去方法を試したが， これらノイズは撮影時期やカメラ，飛行高度，撮影 対象の地形によって変化することが予想されるた め，本研究で行ったノイズ処理方法を一律的に適用 するだけでは対処できない可能性も考えられる．今 後は事例を増やして検証する必要がある． 3 地上基準点および既知標高点での標高精度の検 証 本研究では，探索半径 10.0 m の移動最小二乗平滑 化処理によってノイズが除去された DSM の再現精 度を検証するため，過去から現在にかけて地表面が 変化していないと考えられる 17 地点を対象とした （第 2 図）．検証点における残差の平均二乗誤差（Root
Mean Square Error; RMSE）は 2.01 m，ピクセル当た りでは 0.032 pix であった．高さ方向の誤差は最大で 1.45 m となり，集落中央部では 0.23 m 程度であっ たが，平地ではない起伏の大きい傾斜地の特に山際 において誤差が大きく過小評価することがあった． これらの誤差は，集落を囲む山林内に基準点が無 かったことや，写真周辺部であったため，相対的に を試みた．ノイズ除去フィルターの感度分析のため， パラメータの一つである探索半径（Search Radius） の設定を，5.0 m，7.5 m，10.0 m の 3 通りに変化さ せ最適な値を探索した（第 3 表）． 探索半径を 5.0 にした場合では，無処理の結果と 比較して細かなノイズが除去されたが，ノイズの多 数は残ったままであった（第 5 図 d）．探索半径を拡 大していくにつれて細かなノイズは除去されていき （第 5 図 e，f），探索半径を 10.0 m にした時に概ね ノイズは除去され DSM の等高線がクリアに表現で きた．ただし，平滑化処理を受けているため，前述 の等高線の長さによるフィルタリングと比較すると 等高線が滑らかになりすぎているところも見られ た． 第 5 図 地表面のノイズ除去結果． （a）オルソモザイク画像，（b）無処理，（c）線分延長 100m 以下の等高線のフィルタリング，（d）移動最 小近似平滑化処理の探索半径 5.0m，（e）探索半径 7.5m，（f）探索半径 10.0m．（b）～（f）は 0.5m 間隔の 等高線． 値 た し 用 採 タ ー メ ラ パ 0 . 0 1 , 5 . 7 , 0 . 5 s u i d a R h c r a e S s e Y s l a m o N e t u p m o C s e Y l a i m o n y l o P e s U 2 r e d r O l a i m o n y l o P
Squred Gaussian Parameter 25.00, 56.25, 100.00 e n a l P l a c o L e l p m a S d o h t e M g n il p m a s p U 0 s u i d a R g n il p m a s p U 1 0 . 0 e z i S p e t S g n il p m a s p U 第 3 表 ノイズ除去パラメータ．
高線として表現されているが，一部の樹高の高い樹 木群や家屋等による地物の影響を受けた場所を除く と，概ね圃場の地形を再現できた（第 6 図 a）．検証 地区 B では，南から北に向かって棚田が傾斜してお り，検証地区 A と比較して山際に近い場所に立地し ており，各棚田の面積も狭小であった（第 6 図 b）． 基盤整備前の既知の標高点 34 地点の検証結果は， RMSE で 0.22 m であった．圃場の面積が狭く写真か ら地形の特徴点を捉えられなかったためか，各棚田 の境界と一致しない箇所も見られたが，圃場の地形 は大まかには捉えられていた． Ⅳ 摘 要 本研究では 1970 年代の空中写真と写真測量用ソ フトウェアを使用して，圃場整備前の地形を集落ス ケールで復元し，その再現精度を検証することを目 的とした．過去の空中写真を用いて，SfM-MVS 手 法による過去の三次元地形モデルを構築した結果， 地表面には写真あるいは点群作成時に生じたと思わ れる不規則なノイズが多数出現し，除去が前処理と して必要であることが明らかとなった．本研究では， GIS ソフト上で測定した等高線の線分長が 100 m 以 下のものを除去するフィルタリングによるノイズ除 去，および点群編集ソフトウェア上で移動最小二乗 近似平滑化・再構築によるノイズ除去を試みた．前 者では，細かなノイズが除去され棚田の地形を再現 しているように見えたが，等高線が鋸状になってい る場所も見られた．一方，後者でも細かなノイズが 基準点が多い集落中心部の圃場と比べて高さ方向の 精度が悪くなった可能性が考えられる．水平方向の 誤差は高さ方向よりも小さかった．これらの誤差は 地上標定点の精度や，数の不足も考えられるが，多 すぎることもオーバーフィッティング問題があるた め，適切な数で調整する必要がある．また，今回用 いた空中写真は撮影区域の最南端の撮影コースで あったため，東西方向に撮影された空中写真のみの 解析および検証となった．南方向に隣接する空中写 真は撮影時期が異なることから，今回は解析に使用 することができなかったが，東西方向だけでなく南 北方向を含めた全方向からオーバーラップできる空 中写真が利用可能な場合は今回とは異なる精度にな る可能性もある．しかし，東西方向の空中写真のみ でも対象集落の地形は概ね再現できており，高さ方 向の誤差が最大で 1.45 m でも許容できる目的には使 用できると考えられた． 次に，再現された過去の地形の圃場の区画につい て検証する．本研究で主に対象としている圃場の標 高の再現性を確認するため，過去の地形図から得ら れた標高と今回の空中写真から得られた標高を比較 した．対象集落内の圃場整備計画図を基に作成され た DEM マップが存在する 2 か所の地区を対象に検 証を行った．まず，検証地区 A は北西から南東に向 かって棚田が傾斜しており，地区内の中央部には水 路を挟んで東側は小高い丘が位置していた．圃場整 備計画図に記載された基盤整備前の地形を測量した 48 地点における検証の結果，RMSE は 0.30 m であっ た．空中写真から復元された地形は 0.5 m 間隔の等 第 6 図 空中写真および圃場整備計画図から復元された三次元地形の比較． 検証地区 A（a）および検証地区 B（b）．等高線が空中写真から，背景の濃淡が圃場整備計画 図から復元された地形を，三角点は標高検証地点をそれぞれ表す．検証地区 B は圃場整備計 画図より復元された標高と重複するエリアのみを示している．
2） 渕山律子・大田 健・草場 敬・福原道一：赤 外カラー空中写真と基盤整備前の地形図を用い た浅礫層分布の推定手法，日本土壌肥料学雑誌， 82，134-138，2011． 3） 栗田英治・福本昌人：小型 UAV 空撮・三次元 形状復元技術による傾斜地農地環境の把握，農 業農村工学会誌，84，753-756，2016． 4） 小花和宏之・坂上清一・八木隆徳：RTK-UAV を用いた地形計測の測位性能および省力効果， 地形，40，125-134，2019．
5） Sirbu Flavius, Lucian Dragu, Takashi Oguchi, Yuichi Hayakawa, and Mihai Micu: Scaling land-surface variables for landslide detection, Progress in Earth and Planetary Science, 6, 44, Jun, 2019.
6） Skarlatos, D. and M. Vlachos: Vegetation removal from UAV derived DSMs, using combination of RGB and NIR imagery, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, IV-2, 255-262, 2018.
7） 栗田英治：小型 UAV 空撮・三次元化技術を用 いた農村ランドスケープの復元，沙漠研究， 28，25-31，2018． 8） 小原 洋・浜崎忠雄：基盤整備田の土壌分布と 整備後の微地形との関係，日本土壌肥料学会講 演要旨集，32，125，1986． 9） 清水裕太・松森堅治・望月秀俊：中山間地水田 輪作圃場におけるダイズ・ムギの湿害に及ぼす 過去の地形の影響，日本土壌肥料学会講演要旨 集，64，6，2018． 10） 西村蹊二：空中写真の手引き，財団法人 日本 地図センター，1-54，1977． 11） 石塚直樹・岩崎亘典・坂本利弘：マルチコプタ 型 UAV による熊本地震被災水田の不陸計測と 不 陸 発 生 の 素 因 解 明， シ ス テ ム 農 学，34， 41-47，2018． 12） 早川裕弌・小花和宏之・齋藤 仁・内山庄一郎： SfM 多視点ステレオ写真測量の地形学的応用， 地形，37，321-343，2016． 13） 鳥山和伸：フィールドから展開される土壌肥料 学：新たな視点でデータを採る・見る：大区画 水田における地力窒素ムラと水稲生育，日本土 壌肥料学雑誌，72，453-458，2001． 除去されたが，平滑化処理により等高線が滑らかに なりすぎる結果となった．目的に応じてノイズ除去 方法の選択とパラメータの設定を行う必要が示され た．検証点 17 地点での誤差の検証の結果，高さ方 向は最大 1.45 m となった．平地では 0.23 m と比較 的誤差は小さかったが，起伏の大きい傾斜地の特に 山際において過小評価することがあった．圃場整備 計画図との比較では，棚田群の地形は概ね再現でき ていた．これらの結果，精度の向上は今後必要であ るが，過去の地形図が手に入らない場合に，過去の 地形を復元する方法の一つとして本研究で示した方 法は有用であると考えられた．また，過去の地形図 が存在していた場合でも，紙の劣化が激しいことが 予想され，手作業で判読しデジタル化に要する時間 は非常に大きい．この点についても容易に入手可能 な空中写真から半自動的に地形を復元することは作 業時間も従来法より短縮でき，かつ汎用性も高い． ただし，高さ方向の誤差は最大で 1.45 m となり，こ れらの誤差範囲で問題ない用途には用いることがで きる． 謝 辞 本研究を行うにあたり，地図資料を提供して頂い た，共和の郷・おだの皆様，（農）ファーム・おだ の森原敏昭氏，上野貴彦氏から多大なご協力を頂い た． また，現地調査においては，農研機構西日本農業 研究センター営農生産体系研究領域の川崎洋平氏， 同センター水田作研究領域の浅見秀則氏の協力を得 た．ここに記して謝意を表する．本研究の一部は JSPS 科研費 18K14549 の助成を受けたものである． 引 用 文 献
1） Christopher, Gomez, Hayakawa Yuichi, and Obanawa Hiroyuki: A study of Japanese landscapes using structure from motion derived DSMs and DEMs based on historical aerial photographs: New opportunities for vegetation monitoring and diachronic geomorphology, Geomorphology, 242, 11-20, 2015.
