メインコンテンツまでスキップ

リモートセンシング

概説​

リモートセンシングとは、非接触で対象物の電磁波の反射、放射特性を計測する技術を指します。リモートセンシングという言葉は1960年代にアメリカでつくられた技術用語で、1972年最初の地球観測衛星Landsatが打ち上げられてから急速に普及しました。リモートセンシングによって取得された人工衛星画像や航空写真からは、特定の地物の二次元空間での分布状況が、またそれらのステレオ画像から生成される三次元データからは三次元空間での分布状況が把握できます。

リモートセンシングで使用されるセンサーは、対象物から反射または放射される電磁波を検知し、それが画像化されます。本節では、**地理情報システム（GIS）**を用いた解析でリモートセンシング画像を使用する状況を念頭に置き、必要と思われる知識を紹介します。

基礎知識​

リモートセンシングでは広範囲を対象とするため、センサーを移動体（プラットフォーム）に載せた観測が行われます。本節ではプラットフォームを衛星に限定して話を進めます。

軌道・高度​

衛星では飛行する軌道があらかじめ計画されています。地球観測衛星の軌道は太陽同期準回帰軌道と静止軌道に大別できます。前者は北極、南極の両極の近辺を通る約600-800kmの高度の極軌道であり、衛星は1周約90-100分で周回します。後者は赤道上空の高度約36000kmの軌道で、衛星が地球の自転に同期して飛行することで地球からは静止しているように見えます。

波長​

土木計画学の分野で利用されるリモートセンシングでは、可視光線から赤外線の範囲の波長と、マイクロ波の波長の電磁波が用いられます。前者の電磁波を利用するセンサーを光学センサーと呼びます。

光学センサーでは単位時間、単位面積、単位立体角当りの放射エネルギーである放射輝度が計測されます。放射輝度は入射光のエネルギーによって大きく変動します。その変動を抑えるために、入射光の放射輝度で基準化して得られる反射率が解析に用いられることが多いです。

横軸に観測波長を、縦軸に反射率を記した曲線を分光反射率曲線と呼びます。対象物によって、また対象物の状態の違いによって分光反射率曲線は変化します。例えば、植物は赤の波長帯で反射率が低く、近赤外の波長帯で反射率が高いことが知られています。

図-2：土地被覆と分光反射率曲線

すべての波長帯で反射率を計測できれば対象物を識別できる可能性が高まりますが、大掛かりで費用がかさむセンサーになってしまいます。そのため、光学センサーは、物体の識別に有効といわれている特定の波長帯において計測します。この計測用に設定された波長帯をバンドと呼びます。

光学センサーでは一般的に、4から7程度のバンドが設計されていることが多いのに対し、100以上の多数の狭い波長帯で計測されるハイパースペクトルセンサーというセンサーも運用されています。分光反射率をより細かく検知できるため、物体の識別能力の向上が期待され、ハイパースペクトル画像のさまざまな実利用例が今後増えていくと予想されます。

また、可視光線から赤外線の波長帯の電磁波は大気中の雲によって散乱されるため、光学センサーでは雲が存在すると地表面を観測できません。光学センサーは太陽光を頼りとして計測するため、夜間の撮影は難しくなります。

一方、マイクロ波を出力して、地上の散乱体で散乱された後方散乱の強度と位相を計測して画像を生成するセンサーを**合成開口レーダー（SAR）**と呼びます。レーダーは電波を使った検出と測距を行う技術を指しますが、電波そのものを指して使われることもあります。直下からの角度を意味するオフナディア角が約20-40度の斜め下方向にマイクロ波が発射されます。光学センサーと異なり、大半のSARでは単一の波長での計測が行われます。

一般的にレーダーは雲を透過するため、夜間でも悪天候時でも計測が可能である点が光学センサーとの大きな違いの一つです。ちなみに、レーダー気象学という分野もあるように、マイクロ波が大気中の雨滴で散乱される性質を利用することで、レーダーは雨滴の存在を検出できます。

空間分解能​

空間分解能とは識別可能な地物の大きさを反映する能力を指します。高い空間分解能を有する光学センサーでは数十cmから1m程度の解像度の画像が得られます。一般的に空間分解能が高いと、センサーが観測できる幅が狭くなります。そのため時間分解能が低くなり、特定の地点を撮影した画像を取得する頻度が少なくなります。

SARの空間分解能は光学センサーのものと性質が大きく異なります。SARではレーダーの進行方向をレンジ方向、レンジ方向に直交するセンサーの飛行方向をアジマス方向と呼びます。レンジ方向に走査することで一次元のデータが得られ、それをアジマス方向に飛行しながら繰り返すことで二次元の画像が生成されます。高い空間分解能を達成するためにレンジ方向にはパルス圧縮と呼ばれる技術が、アジマス方向には仮想的に大きなアンテナで観測したように同一地点を観測したデータを統合するアジマス圧縮が適用されます。ここで述べられている圧縮は、情報量を減らすような圧縮という意味ではなく、特定の地点が含まれるデータを畳み込んで合成する処理を指します。このような処理で生成された画像を合成開口レーダー画像と呼びます。圧縮前は数kmから十数kmと低い解像度の画像ですが、圧縮後は数mから十数mの高い解像度を持つ画像が生成されます。

時間分解能​

時間分解能は観測される時間間隔を表し、観測間隔が短いと時間分解能が高いといわれます。

一般的に衛星リモートセンシングでは、数日から数十日で完全に同じ軌道に戻り、この日数を回帰日数と呼びます。観測幅が狭い高空間分解能センサーでは、この回帰日数と同じ間隔で同一地点を撮影できます。一方、空間分解能が低いセンサーでは観測幅が大きく、センサーから見た角度は変動するものの同一地点を毎日観測できます。

米国航空宇宙局 (NASA) の MODISはTerraとAquaというこつの衛星に搭載されている光学センサーです。前者は10:30頃，後者は13:30頃に観測するよう設計されています。MODIS画像の解像度は 250 m、500 m、1 km と決して高くはないものの、毎日昼夜合わせて複数回観測する高い時間分解能に魅力があります。NASAのWeb サイト では MODIS の様々なプロダクトがダウンロードでき、土木計画学の観点では、大気補正済み反射率、地表面温度、植生指数等のプロダクトが利用価値が高いと思われます。

一般的に衛星の飛行軌道は固定されているのですが，災害時等に視線方向を変更できるセンサーもあります。この機能により本来の時間分解能よりも高い頻度で特定地域を集中的に観測できます。

解析に必要な補正や処理​

衛星リモートセンシングでは広範囲を対象とするため、また地球の自転の影響を受けるため、画像に幾何学的なひずみが生じます。また異なる時期に計測された画像を比較する際には、大気状態の違いの影響を除去しなければなりません。本項では、そのような補正や処理を説明します。

物理量変換​

衛星画像には、例えば0から255までの離散的な数値であるDN(digital number)が記録されています。光学画像では、センサー固有の変換係数を適用して放射輝度や、さらに入射光に対する反射光の割合を意味する反射率へ変換できます。

SAR画像では、DNに変換係数を適用して対数化し後方散乱係数を生成できます。次項で述べる二次元主題図作成の際にDNを用いても処理できますが、DN自体は意味がないため、放射輝度や反射率といった物理量に変換してから処理が行われることが多いです。

大気補正​

光学画像では雲が存在すると写ってしまい、その下の地表面が観測できません。可視光線から赤外線の波長では、電磁波は大気分子やエアロゾルによって吸収や散乱されます。水蒸気や窒素、オゾン等の大気分子の粒径は波長に比べて短く、大気分子による散乱は波長の4乗に反比例するレイリー散乱と呼ばれます。そのため、可視光線の中では赤や緑のバンドより青のバンドでの散乱の影響が強くなります。

一方、霧等の水滴や大気中のゴミ等のエアロゾルの粒径は電磁波の波長よりも大きいため、ミー散乱という波長依存性が小さい散乱が生じます。これらの大気による吸収、散乱成分を除去するために、散乱過程を考慮した大気補正処理が行われます。

例えば、画像撮影日時や対象地点の緯度経度、気候タイプ等のパラメーターを利用者が入力することで、あらかじめ用意された大気モデルを通して大気による吸収、散乱成分が計算されます。その影響を除去することで、衛星画像から直接得られる大気上端の反射率から、地上での反射率を推定できます。

また厳密な大気補正ではないものの、光学画像では大気の影響を緩和するためにコンポジット処理が適用されることがあります。

コンポジット処理とは、特定の目的のために、一定期間内に取得された複数枚の画像を合成して1枚の画像を生成する処理を指します。光学画像のコンポジット処理では、一般的に大気の影響を受けると反射率が低下する現象を利用して、雲の影響が最小のときに最高の反射率が観測されると仮定します。例えば1箇月単位でコンポジット処理する場合に、ある画素における1箇月内での最高の反射率を採用して画像を生成します。このようにして各画素の反射率を決定することで、例えば熱帯モンスーン地域でも雲がほとんどない画像が生成できます。その反面、画素ごとに違う日の観測データが記録されていることになります。また、当該期間内で土地被覆に大きな変動が生じないという仮定が前提となっています。そのため、8-10日間と短い期間でコンポジットすることもあります。

一方、SAR画像では雲の影響を受けにくいと述べましたが、後ほど述べる干渉処理を応用して、数mmmm程度の高精度で地盤変動を推定する場合には大気補正が必要となります。高精度推定では2枚の画像から算出される位相差が重要になり、電磁波の伝播の遅延は無視できなくなります。そのため、大気の水蒸気量を見積もった位相の補正が行われることもあります。

幾何補正​

受信信号から生成される衛星画像は幾何学的にひずんでおり、GIS上で他の画像やデータと重ね合わせて解析するには、幾何補正が必要です。単に画像間で位置を合わせる場合と、特定の地図座標系と投影法に合わせる場合に大別できます。前者では基準となる画像を決めて、幾何補正する画像との間で共通して写る地物であるGCP（ground control point）を自動あるいは手動で複数取得します。一連のGCPを用いて投影変換係数を決定し、基準となる画像の形状に合わせた画像が得られます。射影変換式やそれから派生した変換式が用いられることが多いです。

衛星画像を入手すると画像本体のファイル以外に、衛星画像の画素ごとの緯度、経度を算出できる補助データが記述されたファイルが含まれています。このデータを読み取ることで衛星画像を特定の座標系に自動的に変換できます。ただしこの場合でも既存の他の画像やデータと重ね合わせると、幾何学的なずれが残っていることが多いです。そのため、厳密に重ね合わせをする場合にはGCPを取得した精密な幾何補正を実施することが多いです。

衛星画像の1シーンは数十kmkmから百数十kmkmの幅の領域に限られることが多いです。より広範囲を網羅する画像を生成するにはモザイク処理と呼ばれる複数シーンの画像を接合して処理を行います。各シーンは一部分が重なるように撮影されており、GCP検出と同様に対応する地点が検出されます。

図-3：モザイク処理の例

参考リンク

パンシャープン処理​

前述したように、光学センサーには可視光線から赤外線の波長帯において複数のバンドが設計されています。そのような細かいバンドとは異なり、例えば可視光線全体を計測するパンクロマチックというバンドが設けられているセンサーもあります。一般的に、パンクロマチックバンドは波長帯域が広いためSN比が高くなり、高い空間分解能を期待できます。

例えば、2014年に打ち上げられた商業衛星WorldView-3では、青、緑、赤、近赤外等のマルチスペクトルバンド画像では1.24mの解像度である一方、0.45~0.8μmの波長帯を観測するパンクロマチックバンド画像は0.31mの高い解像度を有します。

パンシャープン処理と呼ばれる画像処理の技術を通して、パンクロマチックバンドの解像度を有するマルチスペクトルバンドのカラー画像が生成されます。

もともとのマルチスペクトルバンドから3バンドを選び、RGB画像を生成します。RGB画像からhue(色相)、saturation(彩度)、intensity(明度)で構成されるHSI画像へ変換します。このintensityをパンクロマチックバンド画像に置き換えて、HSI画像からRGB画像へ逆変換します。この結果、パンクロマチックバンド画像の解像度を有する、視覚的に効果のあるカラー画像が得られます。

図-3：パンシャープン処理

参考リンク

二次元主題図作成・地物抽出の原理​

リモートセンシング画像から二次元の主題図を生成するために、特徴空間をなんらかの分類基準に基づき分割し、比較的等質な画素のグループにラベルを付ける分類という処理が行われます。主題図を作成するにはこの分類という処理が欠かせません。

リモートセンシング画像からは土地被覆別に分類できます。一方、土木計画学の分野では、人間の利用形態という観点から整理される土地利用別での主題図が求められるかもしれません。

例えば、空港を例に挙げると、空港自体はアスファルト舗装の道路やコンクリート製の建物等の人工構造物と、草地、裸地、場合によっては樹林帯のような非人工構造物から構成されます。コンピューターを通してリモートセンシング画像から直接分類されるのは、コンクリート構造物、草地、樹林帯のような土地被覆です。複数の土地被覆の集合体である土地利用を生成するには、土地利用を定義した上で高度な処理が求められます。よって、土地利用図が必要な場合には、後述する手順にならって利用者が先に土地被覆図を生成し、その後、土地利用図へ編集する必要があります。

分類手法は、大別して事前に指定するデータの有無と分類する単位の二つの観点から整理できます。まず前者では、分類あるいは抽出の際に対象とする典型的な土地被覆や地物を利用者が指定することがあります。指定されたデータは教師データと呼ばれます。分類手法は教師付き分類と教師なし分類に分けられます。教師付き分類では、与えられた分類クラスを代表する部分を画像からサンプリングして教師データとし、それらから母集団の特徴を推定します。教師なし分類では、特徴の類似した画素等を自動的にグループ化し、それらから特徴が推定されます。教師付き分類手法には最尤分類法、教師なし分類手法にはクラスタリングが代表例として挙げられます。

一方、分類単位という観点では、分類手法は画素単独のデータのみを利用する画素単位の分類手法と、対象画素の一定範囲の周囲を含めたデータを利用するオブジェクト単位の分類手法に分けられます。オブジェクト単位の分類手法の一つである領域分割は、画素周辺の一定範囲における輝度値や反射率、場合によっては形状も含めた同質性を判定します。画素単位の分類手法よりも誤差や局所的な陰影、例えば屋根の部分的な輝度値の違い等にも頑健に同一クラスを生成できます。領域分割によって多少の誤抽出も発生するものの、例えば建物だけを抽出することもでき、効率的に衛星画像を処理する手法として種々の分野で利用されています。

三次元データ生成の原理​

三次元データを生成する原理は光学画像とSAR画像で大きく異なります。

光学画像​

写真測量と同じ原理で、重複して撮影された2枚以上の写真に共通して撮影されている点を特定することで、三次元データを復元できます。ここでは、宇宙航空研究開発機構(JAXA)が打ち上げたALOS衛星に搭載されたPRISM(Panchromatic Remote-sensing Instrument for Stereo Mapping)を例に取り上げて、処理の大まかな手順を示します。2006年から2011年まで運用されたPRISMは観測幅35km、地上分解能2.5m(直下視)で前方、直下、後方の3方向を観測しながら飛行し、数秒間で特定の地域を3方向から観測した画像群(トリプレット)が提供されています。以下、このPRISMのトリプレット画像から三次元データを復元する流れを説明します。

  1. 地上座標が既知の基準点1点を使って、RPCモデル適用時に必要となるバイアス項を決定
  2. 特定の地物を3枚の画像中で検出し、画像座標を決定(ステレオマッチング)
  3. 三つの画像座標と各画像のRPCモデルを用いて地上座標を推定
  4. 3枚の画像が重複する範囲に2、3を繰り返して三次元データを持つ画像を生成

写真測量では、実空間上の点と投影中心、CCD等の像面上の点は同一直線上に存在するという条件、共線条件を仮定します。実際にはカメラ内部のレンズひずみ等のために共線条件が満たされず、内部標定という処理を通して共線条件が満足されるようにセンサーに関係するパラメーターが調整されます。

ステレオ撮影された衛星光学画像を用いた三次元データ復元の現在の主流は、共線条件の概念を発展させた有理多項式モデル(RPCモデル) を用いる手法です。RPCモデルは特定の点の緯度、経度、楕円体高を与えることで二次元画像上に投影される座標を推定します。緯度ϕ\phi、経度λ\lambda、楕円体高hhに対し正規化された緯度ϕn\phi_n、経度λn\lambda_n、楕円体高hnh_nが式(15.21)で、また画像座標uu、vvに対し正規化画像座標unu_n、vnv_nが式(15.22)で定義されるとします。

ϕn=ϕ−ϕ0ϕs,λn=λ−λ0λs,hn=h−h0hs(15.21)\phi_n = \frac{\phi - \phi_0}{\phi_s}, \qquad \lambda_n = \frac{\lambda - \lambda_0}{\lambda_s}, \qquad h_n = \frac{h - h_0}{h_s} \tag{15.21}

un=u−u0us,vn=v−v0vs(15.22)u_n = \frac{u - u_0}{u_s}, \qquad v_n = \frac{v - v_0}{v_s} \tag{15.22}

ここで各変数に対し、ss、00の下付きはスケール項、オフセット項を意味します。式(15.21)、(15.22)を用いて、RPCモデルは式(15.23)で表現されます。

un=f1(ϕn,λn,hn)f2(ϕn,λn,hn),vn=f3(ϕn,λn,hn)f4(ϕn,λn,hn)(15.23)u_n = \frac{f_1(\phi_n, \lambda_n, h_n)}{f_2(\phi_n, \lambda_n, h_n)}, \qquad v_n = \frac{f_3(\phi_n, \lambda_n, h_n)}{f_4(\phi_n, \lambda_n, h_n)} \tag{15.23}

fi(ϕn,λn,hn)f_i(\phi_n, \lambda_n, h_n)は、高々三次の項から構成される多項式で20の係数を持ちます。式(15.23)全体で80個に上る係数はGPSとスタートラッカーを用いて推定される衛星の姿勢と軌道情報から計算され、これらは衛星画像販売元から提供されます。

RPCモデルを用いて三次元座標を推定すると、衛星の位置や姿勢という外部標定要素の誤差に基づくバイアス項が含まれることが知られています。そのため、1.で既知の基準点を用いてバイアス項を決定します。

2.のステレオマッチングでは種々の手法が提唱されていますが、最も基本的な方法の一つが面積相関法です。2枚の画像おのおのにおいて一定の大きさを持つ窓を配置し、窓枠の範囲に含まれる画素の輝度値から両者の類似度を計算し、最も類似度の高い画素を選びます。3.では、ここで得られた画像座標と、RPCモデルから推定される画像座標との残差の二乗和を最小化して地上座標を推定します。

2枚の画像を取得したときの撮影位置間の距離を基線長と呼びます。左右から物を見ることで対象物までの距離が特定しやすいように、光学画像では基線長が一定値以上であることが望まれます。基線長が短いと、推定高さの精度が低下します。

SAR画像​

SAR画像では、ほぼ同じ位置、同じオフナディア角で撮影した2枚の画像を干渉処理することで干渉縞画像を作成し、対象物の相対的な高さを推定します。

図-5に示すように、レーダーは異なる時期にM、Sの二つの位置で観測を行うものとします。衛星SAR画像を用いた干渉処理では、一つのアンテナで繰り返し観測された2枚の画像を利用することが多いです。基線長BB（数百mm～1km1km程度）に比べてセンサー-対象物間の距離RRは大きい（700700～800km800km）ので、二つの入射角は等しいと仮定できます。

三次元データの復元に焦点を当てるため、2時期の間には有意な地盤変動は生じないと仮定して話を進めます。

図-5：SAR画像を用いた地物の高さ推定

以下、干渉処理の大まかな手順を示します。

  1. 2枚の画像の相対位置合せ
  2. 特定画素における位相差の計算
  3. 位相差から基準面での位相差を除去
  4. 位相差のアンラッピング
  5. 4.の位相差から基準面からの相対高を計算
  6. 2枚の画像が重複する範囲に5.を繰り返して三次元データを持つ画像を生成

まず1.では2枚のうち片方の画像を固定し、もう片方の画像の位置を合わせます。画像の付属データを用いて位置合せすることも可能であり、自動処理する商用ソフトウェアも提供されています。

以下、2.以降の処理を述べていきます。特定の地点に対するレーダーからの二つの距離の差（光路長差）は、レーダーと対象物間の往復を考えると図-5に示す距離δ\deltaの2倍です。

光路長差に対応する位相差φ\varphiは、基線長BB、波長λ\lambdaを用いて、下式で示されます。

ϕ=2⋅2πλδ=4πλBsin⁡(θ−α)(15.24)\phi = 2 \cdot \frac{2\pi}{\lambda}\delta = \frac{4\pi}{\lambda}B\sin(\theta - \alpha) \tag{15.24}

一方、図-5から、基準面からの特定地点の標高zz、マスター画像撮影時のレーダーの飛行高度HH、RRとの間には、次の関係式が成り立ちます。

z=H−Rcos⁡θ(15.25)z = H - R\cos\theta \tag{15.25}

φ\varphiはHHとRRを変数に持つ関数と考えて、式（15.26）の関係式を導けます。

dϕ=∂ϕ∂θdR+∂ϕ∂zdz(15.26)d\phi = \frac{\partial\phi}{\partial\theta}dR + \frac{\partial\phi}{\partial z}dz \tag{15.26}

式(15.26)の右辺第1項は軌道縞による位相の変動分、第2項は地形縞による位相の変動分を表します。

衛星の軌道情報から理論的に計算できる軌道縞を差し引くことで、第1項を差し引くことができます。

第2項に着目すると

∂ϕ∂z=∂ϕ∂θ∂θ∂z=4πBcos⁡(θ−α)λRsin⁡θ=4πB⊥λRsin⁡θ(15.27)\frac{\partial\phi}{\partial z} = \frac{\partial\phi}{\partial\theta}\frac{\partial\theta}{\partial z} = \frac{4\pi B\cos(\theta - \alpha)}{\lambda R\sin\theta} = \frac{4\pi B_{\perp}}{\lambda R\sin\theta} \tag{15.27}

と変形できます。ここでB⊥=Bcos⁡(θ−α)B_{\perp} = B\cos(\theta - \alpha)とし、垂直基線長と呼びます。

式(15.27)は位相の変化量と高さの変化量の関係を示しています。ある基準面のある地点との相対的な位相差をΔϕ\Delta\phi、高度差をΔz\Delta zとすると、式(15.27)から

Δz=λRsin⁡θ4πB⊥Δϕ(15.28)\Delta z = \frac{\lambda R\sin\theta}{4\pi B_{\perp}}\Delta\phi \tag{15.28}

が得られます。

つまり、3.で特定の基準面上の地点との相対的な位相差を計算した後で、式(15.28)に示すように2地点間の相対高さΔz\Delta zを2点間の位相差Δϕ\Delta\phiから推定できます。

しかし、式(15.28)のΔϕ\Delta\phiは本来の位相差とは限りません。SARでは反射波の強度と位相が記録されますが、位相は−π-\piからπ\piの範囲で表現されています。この変換をラッピングと呼び、ラップされた位相から元の位相を復元する処理をアンラッピングと呼びます。したがって高さの推定の前に4.のアンラッピングが必要になります。ただし、正確な位相を復元できればセンサーと対象物間の距離を推定できますが、実際にはアンラッピングの際に誤差が含まれてしまいます。