南極地球物理学ノート No. 6 (2011.11.20)

干渉合成開口レーダーで決めた氷床接地線１


澁谷和雄・小澤　拓

Keyword: 氷床接地線 (Grounding Line), 合成開口レーダー (SAR), 干渉SAR, 棚氷の海洋潮汐変形,
Zubchatty棚氷, エンダビーランド



１．氷床接地線 (Grounding Line)

空中写真、衛星画像を見ただけでは、南極大陸の海陸境界は判然としない場合が多い。棚氷が広がり、一面真っ白な世界で海の色を識別できるのがはるか沖合いだったりするからである。しかし、太陽・月など天体が及ぼすニュートン力の現れ方は海と陸では厳然と異なる。昭和基地付近の場合、陸上での重力潮汐による上下変位は～40 cmであるが、海での潮位変化は～2 mに達する。そのため、陸に着底している氷床(grounded ice sheet)と海に浮いている棚氷(floating ice shelf)の境界域では、上下変位の遷移域(transition zone)が現れるはずである。

この遷移帯をリモートセンシングで同定できないだろうか？そして、海陸境界を決められないだろうか？

衛星合成開口レーダー（SAR: Synthetic Aperture Radar）データの干渉処理で遷移帯の同定が可能なことをGoldstein et al. (1993)がロンネ棚氷のラトフォード氷河において、初めて示した。一方、Rignot (1996)はグリーンランドのぺテルマン氷河を例にとり、海陸境界の最も内陸よりを氷床接地線(GL)と定義した。

このノートは東南極エンダビーランドの棚氷での我々の干渉SAR適用結果である。



２．合成開口レーダーを用いた遷移帯の検出

2.1. SAR

SARは、1951年6月、アメリカのGoodyear Aircraft Corporationの技師（数学者）であったCarl A. Wileyが最初に考案し、1952年にはその概念検証システムを製作したとされる能動型(active方式)マイクロ波センサーシステムである。飛翔体から地表面に向けてマイクロ波を照射し、その跳ね返り（後方散乱波の強度と位相）を同一の飛翔体で受信・記録し、解析する仕組みである。SARの画像生成原理は複雑なので、また、本題ではないので、地球の曲率を無視しうる低高度を飛ぶ航空機SARを例に観測の原理を図示する（図１）。

図１において茶色いハッチのかかっている部分が地表面である。この範囲において、航空機（黄色）の進行方向をアジマス(Azimuth)方向、それに直交して電波の放射 される方向をスラントレンジ(Slant-range)方向と呼ぶ。今の場合、進行方向左側に電波が放射されるので Left Looking Radar である。 アジマス方向に ℓ m長のアンテナを持ち、斜め下向きθ°(オフナディア角と呼ぶ）の方向に電波を照射して地表からの後方散乱波を受信する場合を考察する。

news06_01図
図１. (SARの観測原理)

SARの画像化における主要なパラメーターは(1)空間分解能、(2)観測波長、(3)入射角である。アジマス方向の最高分解能はアンテナ帳の半分 ℓ /2mであるが、ふつう、ルック処理という重ね合わせノイズ低減化を行うので分解能は1/3 - 1/6に低下する。スラントレンジ方向分解能ΔRは、信号の帯域幅fBで決まり

ΔR　= c / (2fB)　　　　　　　　　　　　　　　(1)

である。ΔRに対応したグランドレンジ方向の分解能Δｘは

Δｘ　= ΔR　/ sin i　(i は入射角)　　　 　(1)'

である。アジマス分解能とレンジ分解能で決まる画像の最小単位を画素あるいはピクセルと呼ぶ。散乱強度と位相はピクセル単位で記録される。


2.2. 干渉SAR (InSAR)

干渉SAR (SAR Interferometry) もしくはInSARは1970年代にアメリカで開発された技術である。図２にその概念図を示す。図中、航空機の進行方向は紙面に垂直で、A1が１回目の、A2が２回目のSAR観測位置である。また、地上の標高Ｚの位置にターゲットTが存在すると仮定する。Tが１回目観測と２回目観測の間にT1 → T2へ移動した（変動があった）として、A1 －T1間の見通し距離（slant range）ρ1がA2 －T2間のslant range ρ2に変わったとする。InSARによって得られるのはslant rangeの差、すなわち、δ =ρ1 ― ρ2である。この差は観測波長に応じた位相差(φ)という形で最終的に表され、

φ = 4π (B sin (θ - α) + dD)/λ         (2)

Z = H - ρ1cosθ       　　　　　　　　 　(3)

である。

図２のべクトル成分図においてdDは

dD = δ2 + δ3 = Dx sinθ - Dz cosθ      (2)'

(Dx, Dzはxz平面内の変動量成分)

で、１回目と２回目の観測の間に地表面で発生した変動量の視線方向成分である。また、Hは航空機の飛行高度、λは波長、Bは軌道間距離(Baseline distance)、αはBが地表面（水平面）のなす角度である。地表が平面近似できる航空機SARの場合、地表への入射角i は近似的にオフナディア角θに等しい。

news06_図02
図2.（InSARの原理）

φに含まれるλ, B, H, αは既知情報とみなせるので、φの全微分は

  news06_数式01 (4)

と書ける。ここで、第一項を軌道縞（orbit fringe）、第二項を地形縞（topographic fringe）、第三項を視線方向変動縞（ここでは単純にdisplacement fringe）と呼ぶ。

軌道縞は

news06_数式02          (5)

と書けて、計算で求めることができる。Bpは図２に示すように軌道間距離の垂直成分を意味する。もし２つのSAR画像を取得した時の軌道が互いに平行であれば、航空機の進行方向に平行な軌道縞が現れる。

地形縞は次の式で示される。

 news06_数式03     (6)

(6)式の両辺をｚで積分した場合を考えると、地表面の地形（すなわち標高ｚ）の変化に応じて位相変化（すなわちφの変化）が現れることがわかる。これを空間的に求めるとDEMが生成されるが、その際、Bpが地形縞１フリンジあたりの変化量（高度変化量）を決めるパラメーターになる。逆にDEMがあれば、(6)式による位相変化をシミュレーションで計算できる。地形縞は地表面の標高に比例した位相変化によって生じるから等高線に類似した縞模様になる。

以後は飛翔体が高高度を飛ぶERS-1/-2  (European Remote Sensing satellite -1, -2) 衛星の場合を考えよう。この場合、地球の曲率を考慮しなければならない。図３においてslant rangeの差ρ2,1 ―ρ1とρ2,2 ―ρ1は等しくならない。地形縞は(6)式の代わりに

news06_数式04    (6)’

となり標高Zの代わりに楕円体高hを、オフナディア角θの代りに入射角としてρ1に対する値i を適用しなければならない。

(2)(3)式から(5), (6)あるいは(6)’式により軌道縞と地形縞を除去すると次の変動縞成分が残る。

  news06_数式05      (7)

この成分は地殻変動や地形変化に相当した成分である。

ERS-1/-2衛星の諸元を表１にまとめた。ERS-1/-2 SARは Right Looking Radar で、アンテナ長 ℓ =10 m, オフナディア角θ = 20˚、入射角i = 23˚、fB = 15.55 MHz、λ= 5.7 cm、飛行高度H = 785 kmである。これらの数値を基にノイズ低減化やフィルター操作に伴う劣化を考慮して、Azimuth 及び Ground-range分解能は、ともに30 mと与えられている。

news06_図03
図3. (地球の曲率が無視できない場合の地形縞表現。BPERPは基線長の直交"perpendicular"成分で、さらに省略してBpとも書く。)

(7)式において分母が波長 (5.7 cm) であるから、地表面がマイクロ波の視線方向に半波長2.8 cm変化すると干渉SAR画像の位相差は２π変化することがわかる。変動縞は等変動量線を示す。

実際の干渉SAR処理では、個々のSAR画像についてピクセル単位で位置合わせを行い、ピクセルごとに受信散乱波の位相差を計算する。得られた位相差は、基本的には上記３成分の和になる。しかし、実際に得られる干渉SAR画像では、各干渉成分にノイズが重なっている。あるピクセルからの散乱波強度が弱かったり、散乱強度を特徴づける誘電率の変化が大きすぎたりして位相の相関が失われると地球物理学的に意味のある特徴が得られなくなる。

同一画面であっても変動は小さいが標高差が大きい場所、標高差は小さいが変動が大きい場所で見かけのフリンジパターンが似かよる場合がある。干渉SAR画像はパッチワークのように部分、部分に着目する必要があり、先見情報がゼロだと解釈に無理が生じやすい。また、干渉処理によって得られる位相差は－π～＋πで計算されるので、変動縞成分を求める時、２πの整数倍の不確定性を持つ。この不確定値はアンラッピング(unwrapping)という処理で計算し、位相復元しなければならないが、ここでは触れない。



３．入り組んだ氷床・氷海域での実際の干渉SAR適用例

SAR干渉処理は専門ソフトウェア会社が作成したVexcel 3-D SAR プログラムを用いて行った。この処理ではlevel 0 CEOSフォーマト、SLC処理、マルチルック処理といった専門用語がからむ複雑な操作が必要になるが、それらの説明は専門書に譲り、得られた結果の解釈に焦点をしぼって説明する。


3.1. 適用域とデータ

図4は、解析を行った南極・エンダビーランド、アムンゼン湾域の地形図である。ここには露岩域に挟まれてWyers Ice Shelf, Zubchatyy Ice Shelf という棚氷が発達している。この地域の黒枠部分をERS-1の３日回帰周期で観測した。図５(a)が1991年12月７日、図５(b)が1991年12月10日撮像のデータから散乱強度を取り出し、濃淡表示したものである。最強の散乱強度に白、最弱の散乱強度に黒を割り当て、グレースケール表示してある（このような色の割り当て方が一般的である）。

news06_図04
図4. （エンダビーランドでのSAR解析適用例）

衛星軌道情報を用いて軌道縞を除去したが、含まれる誤差のため、取りきれない平行線状の軌道縞が残ったので、露岩域の海岸線が海抜０ｍであるとして、そこでの位相が同じになるように試行錯誤的に軌道縞を決定し除去した。

図５(a), (b)の取得データに対応した衛星位置間のBp = 67 mであった。これに対応して地形縞は楕円体高が139 m変化するごとに位相が２πラジアン変わることになる。高い分解能を持つDEMが存在すれば地形縞を予測計算し、除去することが可能である。最近（～2010年）、ASTER GDEMと呼ばれる1秒角分解能（～30 m）のDEMが南極沿岸域でも使用できる目途がつき、今後の活用が期待されている。

news06_図05
図５. （３日回帰パスでの散乱強度変化、Wyers棚氷での変化が著しい）

図６は軌道縞を除去して得られた(a)干渉SAR画像と、(b)２枚のSAR画像の複素相関値を色相（青が０－ピンクは0.5―黄は1.0）に分けて表現した相関値図である。相関が良い(0.5以上)黄色及びピンク系統の地域では明瞭な干渉縞が現れているが、0.2以下の青色系統地域では不明瞭である。 Zubchatyy棚氷ではある程度良い相関値が得られているが、Wyers棚氷は散乱強度の変化が大きく、1991年12月10日の撮像データでは散乱が殆ど得られなくなってしまった。

news06_図06
図６ (a) 干渉SAR画像、と(b) ２シーンの複素相関値の色相表現

散乱強度は雪の含水率変化によって大きく変わる。Zubchatyy棚氷が日射を受けにくい山地の南側にあるのに対して、Wyers棚氷は日中の日射を長時間受ける山地北側に位置している。夏場の３日間で散乱特性が急激に劣化してしまい、著しく相関が低下してしまったことが示唆される。

図６(a)の干渉SAR画像では地形縞と変動縞が重なって現れている。露岩域の地形が３日間で変動してしまうことはない。露岩域での変動縞は殆どが地形縞起源で、等高線に類似した干渉縞になっている。露岩域はB―M―Aで区画された左上の地域が主に対応していて、６フリンジ分の地形縞があるとすると標高差は～830 mになる。

Zubchatyy棚氷上で得られた干渉縞を拡大してみよう（図７）。この棚氷はオーストラリア測量局が作成した地形図（1:1,000,000）を参考にすると、高低差は40 m程度と思われる。40 mの高度変化では地形縞の位相はπ/2ラジアン (1/4サイクル)の変化におさまるが、実際にはZubchatyy棚氷中央部で８サイクル以上の干渉縞が得られていて、これらは概ね、棚氷の変動に起因する変動縞とみなせる。

news06_図07
図７. （Zubchatyy棚氷上の干渉縞）

そこで位相のアンラッピング処理を行いピクセル間の相対的な位相差を復元した。復元位相に含まれる誤差は各要因を、大きく見積もっても合計で1/2 – 1波長（2.8 – 5.7 cm）である。得られた位相差（変動縞）を地表が衛星の視線方向に2.8 cm（半波長）変化するごとにπラジアン変化することを用いて変動量に変換した。水色から水色までの色調の変化が２πラジアン、すなわち１波長5.7 cmのslant range変化に対応している。白線が海陸境界と思われ、そこに近づくほど、色調バンドの繰り返しが密になっていることがわかる。つまり、短い水平距離の間に大きな上下変位が生じている。中央部では間隔が広く、南北方向にゆるやかな傾斜の変動になることが見て取れる。画面を４分割したときの右下付近に北東―南西に連なる密な環状のフリンジパターンが見られるが、これは島を取り巻く海岸線での変動がゼロになるように、海洋潮汐に伴う変動縞が局所的に集中した結果と思われる。


3.2. 海洋潮汐による棚氷の変形モデル

Holdsworth (1977)は図８のように半無限長、厚さ2Hで単純化した氷板上で、接地線から x の距離にある点X の時刻 t における変位 X(t, x) が

X(t, x) = waζ(t) { 1 - e-βx(cosβx +　sinβx)}        (9)

であらわされることを示した。ここで、waζ(t)は最大振幅をwaとする海洋潮汐変位、βは氷の物性定数、氷厚、重力加速度などによって決まるdamping factorで、氷板を弾性(Elastic)モデルで近似すると、

βE4 = 3ρw g (1- μ2)/ 8EH3 (10)

である。ここで棚氷の尤もらしい物理定数として、μ = 0.3, ρw = 1.025 Mg/m3, g = 9.8 m/s2, E = 2.7 x 103 MN/m2 を代入し、HをパラメーターとしてβE を見積もると、

H = 200 m →βE = 3.55 x 10-3 1/m →1/βE = 282 m
H = 1000 m →βE = 1.06 x 10-3 1/m → 1/βE = 943 m

が得られる。1/βE が遷移帯（transition zone）の特徴的長さである。

Holdsworth (1977)は氷板を遷移クリープモデルで近似した場合、(10)式はもっと多くのパラメーターを含む複雑な式に置き換わることを示し、βC (Creep model)を推定したが、ここでは立ち入らない。

news06_図08
図８. （Holdsworthによる棚氷変動モデル）

3.3. 海洋潮汐による棚氷変形から求めた氷床接地線

図７の最大潮位差を示す点Mでは潮汐による変形がゼロ（damping free）と考え、Mを通る２本のプロファイルAA’, BB’において変動量を見積もることにする。図９(a)の実線はプロファイルAA’に対応した変動で、図９(b)の実線はBB’に対応した変動である。そして交点Mでの上下変位は41.5 cmと求められた。一方、図９(c)は点MにおけるORI96海洋潮汐モデル(Matsumoto et al., 1995) による20日間の理論潮位変動である。SARデータを観測した1991年12月７日05:17UT（図の点P）と1991年12月10日05:17UT （図の点Q）を破線でマークしたが、その理論潮位差は35.2 cmである。ORI96モデルの元になったTOPEX/POSEIDON衛星のレーダー高度計が65˚S以南をカバーしていないにも関わらず、InSAR推定潮位差（41.5 cm）と海洋潮汐モデル予測値（35.2 cm）の一致度が良いことは注目に値する。

news06_図09
図9 (a) プロファイルAA'に対応した棚氷変形、(b) プロファイルBB'に対応した棚氷変形、(c)ORI96海洋潮汐モデルによるM点での理論潮位変動

実際の干渉SAR画像についてprofile AA’, BB’でのtransition zoneの形態をあらわす曲線を、(9)式の弾性体でモデル化した尤もらしい氷板変形のdamping factor βで現わしてみると、βa = 6.5 x 10-4, βb = 1.0 x 10-3となる。この場合、βa に対応する弾性氷板モデルの特徴的長さは1540 mでそれに対応する氷厚2Hは3846 m, βb に対応する特徴的長さは1000 mでそれに対応する氷厚2Hは2166 mになる。



４．本研究の意義と限界について

本ノートは1999年に測地学会誌に発表された小澤拓らの「干渉合成開口レーダから得られた東南極Zubchatyy棚氷の海洋潮汐による変形」をもとにまとめたものである。 我々のInSAR研究の端緒の時期でもあり、現時点で考えると初歩的解析に見えることは否めない。

この研究の限界は大きく３つあって、(1)地形縞と変動縞の分離が出来なかったこと、(2)棚氷のgeometryはおそらく１次元の半無限平板均一弾性（あるいは遷移クリープ）モデルでは単純化できないこと、(3)視線方向の変位を上下変位とみなし水平変位を無視したこと、である。

(1) について、南極域では現在でも同一地域で繰り返しascending, descendingの複数パスを受信できることはまれで、データ取得面からは2011年になっても顕著な改善はみられない（狭い地域を繰り返し観測するより、回数は少なくともできるだけ広い地域を観測することが優先されてきた）。しかし、地形縞と変動縞の不分離は氷床接地線の同定という意味では大きな障害ではない。画像全体で１－２サイクル位相誤差が生じてもフリンジバンドの幅、位置に大きな修正は不要だからである。
なお、Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)という光学センサーで撮像した衛星ステレオ画像を用いたGlobal DEM (GDEM) の整備が進んでいる。GDEMは水平分解能約30 m、高さ精度約20 mで、南極域でも使用できることが最近、Doi et al. (2011)によりわかったので、今後の氷床変動研究への活用が期待できる。

(2)については、海洋潮汐モデルの検証が目的ではないので、理論予測（35.2 cm）がSAR変動実測（41.5 cm）とほぼ同等の値を示したことで良しとすべきであろう。つまり、密なフリンジバンドは間違いなくtransition zoneを示すと考えて良い。しかし、 氷厚2Hの2000 ～4000 mはどのように考えても過大見積もりであろう。後のノートで示す白瀬氷河の氷舌部のGL付近でもせいぜい300～800 mである。

(3)については、棚氷は大陸氷床が海に流入する（あるいは接地線が内陸側に後退する）ことで発達するが、流入する氷床に押されて棚氷内部で水平移流が生じる可能性が高い。図７(a)ではdD ～Dzとしたが、本来は(2)'式において、DxとDzを分離しなければならない。もし、変動量dDがDzではなくすべて水平成分Dxに起因するとしたら、M点における流動量のrange方向成分は画像の右方向（東方向）に97 cm/ 3days = 32 cm/dayとなり、有り得ない値ではない。Slant range方向の変動量を、上下、水平成分に分離するためには、Interferogramを２ペア解析する必要があるため、同一エリアについて適切な４シーンのデータ取得が必要である。

以上、弱点はあるがしかし、大局的には大陸氷床縁辺部のInSAR画像にどこでも狭いフリンジバンドが出現するので、これが接地線付近のtransition zoneを示すことに疑いはない。



引用文献

小澤拓・土井浩一郎・渋谷和雄(1999): 干渉合成開口レーダーから得られた東南極
Zubchatyy棚氷の海洋潮汐による変形。測地学会誌、第45巻第３号、165-179頁。
Doi, K., Bassler, M., Dietrich, R., Shibuya, K., Yamanokuchi, T., Nakamura, K.,
Omura, M., Koike, K., 2011. Application of ASTER GDEM to estimating flow
rate over the Antarctic ice sheet. J. Geod. Soc. Jpn, 57(2), 61-70.
Goldstein, R.M., Engelhardt, H., Kamb, B., Flolich, R.M., 1993. Satellite radar 
interferometry for monitoring ice sheet motion: application to an Antarctic ice
stream. Science, 262, 1525-1530.
Holdsworth, G. (1977): Tidal interaction with ice shelves. Ann. Geophys., 33,
133-146.
Matsumoto, K., M. Ooe, T. Sato, J. Segawa (1995): Ocean tide model obtained from
TOPEX/POSEIDON altimetry data. J. Geophys. Res., 100, 25319-25330.
Rignot, E., 1996. Tidal motion, ice velocity and melt rate of Petermann Gletscher,
Greenland, measured from radar interferometry. J. Glaciol., 42, 476-485.


本ノートで省略したSAR, InSAR解析の詳細は以下の教科書に詳しい。

大内和夫、2004. 合成開口レーダの基礎―リモートセンシングのための、東京電機大学
出版局、東京。ISBN4-501-32330-2.
Massonnet, D., Feigl, K. L., 1998. Radar interferometry and its application to
changes in the Earth’s surface. Rev. Geophys., 36(4), 441-500.
doi:10.1029/97RG03139.
Hanssen, R. F., 2001. Radar interferometry: Data interpretation and error analysis.
Kluwer Academic, Amsterdam. ISBN 0-7923-6945-9.

SARは軍事技術として開発された。初期の開発状況が、Wikipedia
http://en.wikipedia.org/wiki/Synthetic_aperture_radar /の第４章Origin and early
development (ca 1950-1975)に記載されていて、発明者としてCarl A. Wileyの名がある。
(2012.01.09参照)。 ASTER GDEMについては下記の検証レポートを参照のこと。

ASTER GDEM Validation Team, 2009. ASTER Global DEM Validation Summary Report,
http://www.ersdac.or.jp/GDEM/E/3.html (2012.01.09参照)



Q and A

Q1:  2.1.にアジマス方向の最高分解能はアンテナ長さの半分という記述がありますが、アンテナ長を短くして分解能を上げた方が良いのではないですか？
A1:  SARのパラドクスとして知られる事柄ですが、アンテナ長が短くなるとビーム幅を鋭く出来ず、先端出力も弱くなります。そのため、後方散乱強度も弱くなり、信号が雑音に埋もれてしまうので、サイズを小さくするのには現実的な限界があります。



Q2:  氷床接地線はデータベース化しないと意味がないのでは？
A2:  おっしゃる通りです。干渉SAR画像では緯度・経度情報が消えてしまうので、その復元が必要です。データベース作成は南極地球物理学ノートNo. 7で取り上げます。



Q3: ３頁に「図5(a), (b)の取得データに対応した衛星位置間のBpは67 mであったが、これに対応して地形縞は楕円体高が139 m変化するごとに位相は２πラジアン変わることになる。」とありますが、どうして139 mなのかがわかりません。
A3: ２πラジアン位相変化に対応する高度変化は(6)'式によりdh = λρ1sin i　/2Bpで計算されます。Off-nadir角θが20°なので、ρ1 = 785 km/cos 20˚ = 835382 mです。入射角i　が23°なので、dh = 0.057 m x sin 23˚ x 835382 m / (2 x 67 m) = 138.8～139 mです。



表１.　ERS-1. ERS-2 搭載SARの諸元

Launch 

April 1991

Orbital altitude

785 km

Orbital inclination

98.5˚

Frequency

 5.3 GHz (C-band)

Wavelength 

 5.7 cm

Polarization  

VV

Off-nadir angle

 20˚

Incidence angle

 23˚

Swath width 

100 km

Azimuth resolution 

30 m

Range resolution 

30 m

Peak power

4.8 kW

Chirp bandwidth

15.55 MHz

Antenna size      

 1 x 10 m