干渉SARにおける軌道推定と大気遅延分布の特定 Refinement of the Satellite's State Vector and Identification of the Atmospheric Delay Pattern in the SAR Interferometry 小林茂樹(九州東海大学・工学部・宇宙地球情報工学科)・大塚彰（農業研究センター）・瀬古弘（気象研究所） Shigeki KOBAYASHI* ( Dept. of Space and Earth Information Technology, Kyushu Tokai University ), Akira OTSUKA ( National Agriculture Research Center ) and Hiromu SEKO ( Meteorological Research Institute ) shigeki@ktmail.ktokai-u.ac.jp

Abstract: L-band SAR interferograms are easily obtained over the wide range of height despite vegetation or slope. In general, however, remarkable phase pattern due to the atmospheric excess path delay remains near the summit of high mountains on differential interferograms, even if errors in the satellite’s state vector are precisely corrected. Atmospheric effect contains two components. Pattern 1) Residual topographic phase pattern correlated with height in the first order, which is due to the difference in the static structure of water vapor at two SAR observation times: This pattern can be verified by calculation of relationships between height and the atmospheric delay using radiosonde data. Pattern 2) Complicate phase pattern (e.g. mountain wave, foehn phenomenon) which not correlated with topography: For example, strong atmospheric noises, which appeared as phase undulation with ~15 km in wavelength, were observed by using SAR data acquired just passing or approaching of developing cold front. Such wavelike patterns of phase noises may be caused by mountain wave blowing over high mountains. We developed a meteorological method to quantify the water vapor effect by 3-dimensional numerical simulation of atmospheric dynamics. The initial atmospheric condition was given using the radiozonde data of the nearest weather station. The results show that a mountain wave, a standing wind, appeared on the lee of the mountain. The downward (upward) winds in the mountain wave corresponded well to the relatively dry (wet) area on the differential phase pattern. Downward wind decreases the path delay, upward increases. １．はじめに 干渉SAR画像の中から大気遅延による位相パターンを特定するためには、以下の条件が不可欠である。 (1)　衛星軌道が精密に決定されていること（できれば、干渉縞を用いずに軌道情報が決められるとよい） (2)　軌道縞や地形縞を補正する段階において、衛星軌道情報と地球楕円体およびDEMとの幾何学的関係が精密に考慮されていること しかし、実際には、特にJERS-1干渉SARの場合には以下の問題点がある。 (1)　JERS-1の位置情報は数10ｍ以上の誤差を含んでいること。特に衛星時刻の誤差が著しい場合があるし、（リーダーファイルに記載された）速度情報の定義も不完全であること。 (2)　通常の差分干渉処理ではベッセル楕円体に基づくDEMが用いられるが、ジオイド高の考慮はせいぜい解析領域内の１点でゲタを考慮する程度に留まっていること。 　これまでのJERS-1干渉SARでは、もともと軌道情報の精度が悪いので時刻誤差や地球の形状に関する厳密な考慮があまりなされてこなかった。しかし、衛星軌道精度が格段に向上すると期待されるALOS時代の干渉SAR解析に向けて、可能な限り???の問題に忠実な解析を追及すべきである。ALOS時代には、?の問題でさえ解析誤差の”いいわけ”にできなくなる可能性もある。逆に、軌道情報と整合性のある干渉処理を追求することで新しく見えてくるものがあるに違いない。その１つがまさに大気遅延による位相誤差である。 　本論文では、まずJERS-1干渉SAR の精密な軌道推定方法について議論し、軌道情報が悪いとされるデータペアの中から良質なものを判定する。次にそのような軌道精度のよいデータペアを使った干渉画像に見られる大気遅延パターンの特徴について議論する。最後に、高層気象データ（ゾンデデータ）を用いた大気運動の理論計算により水蒸気の分布量を求め、大気遅延による位相パターンの特定を行う。 ２．JERS-1干渉SARにおける軌道推定方法 以下の方法でJERS-1の軌道の信頼度を判定する。 (1) リーダーファイル内の軌道情報のうち、位置ベクトルのみを使用する。 (2) 位置ベクトルを何組か用いて、Orbit propagator 法（ケプラー運動に関する軌道要素の計算）を応用することで、Master、Slave双方の軌道情報を精密に補間し直す（例えば[1]）。これを初期条件とする。 (3) Master振幅画像とSlave振幅画像の相互相関を計算することでTiepointを決定する。と同時に、それぞれの軌道情報に基づいて、Tiepointの座標値のずれ（座標誤差）を求める。一般にJERS-1の場合には、時刻精度が悪いので（Azimuth方向の）座標誤差が10数ピクセル以上にもおよぶことがある。あるいは、そもそも衛星の姿勢が悪いとRange方向の座標誤差も大きくなるので注意が必要である。 (4) 座標誤差を最小にするように、Slaveの軌道を再決定する。ただし、軌道要素をすべて（MasterとSlave双方の位置ベクトル成分６個）決めようとしてもほとんど解は収束しない。実際的には、(1)未知数２個（時刻誤差とその分の地球自転量の補正）、(2)未知数３個（時刻誤差、軌道面の傾き、軌道半径）の場合を当てはめて、誤差の最も小さいモデルを採用する。?述べた座標誤差が大きい場合には、(1)(2) の方法を適用しても解の精度が改善されない。このような場合には、残念ながら、軌道に忠実な干渉処理はできないので、例えば軌道縞を多項式曲面で近似して人工的に取り除くことなどしかできず、特徴的な地殻変動のパターンを抽出できたとしても、例えば地形に相関するような残差地形縞成分が大気遅延によるものなのかどうかは判定できない。 (5) (4)で求められた軌道要素に従って、軌道縞と地形縞を理論計算するフォワードモデリングの手法により干渉画像を精密にフラットニングする（Slave 軌道をさらに微調整する）。このとき、基線長の決定誤差が、数％以下のもののみを信頼できるデータペアだと判定する（もしくは、?の過程を繰り返すことで誤差を数％以内に収束させる）。この段階では、原則的にフルシーンの軌道縞の勾配が極小になることをもって正しい軌道とする。具体的には、感度の大きい垂直基線長BperpとYaw角の２つを未知数にしてフラットニングを行う（平行基線長Bparaが影響を与え得る特殊な場合もある）。 ３．大気遅延による位相パターンの特徴 精密な軌道推定を行うことで初めて大気遅延による位相パターンが区別できるようになる手順について具体的に説明する。 今もし大気遅延を考えないとすると、地殻変動が全くない地域における軌道縞間隔と地形縞間隔は衛星と地球楕円体との幾何学的な関係のみで決まるはずである（実際には、投影基準とした地球楕円体の座標系のずれやジオイド高の誤差の影響があるがここでは無視する）。軌道縞間隔dx/dphiと地形縞間隔dh/dphiは、簡単な幾何学を仮定すると、 dx/dphi=λR/(2Bperp cosθ)=100000/Bperp, dh/dphi=λR sinθ/(2Bperp)=50000/Bperp [m/2π]で表される。ただし、λは波長で約23.5cm、Ｒはスラントレンジで約700km、θは入射角で約35-42度(JERS-1の場合)、Bperpは垂直基線長[m]である。一般には地表は凸凹しているので軌道縞間隔と地形縞間隔の両方を同時に使って（Ground Control Pointの位置・標高の情報を利用して）軌道を推定することになる。ここで、上の関係式からも容易に想像できるように、（標高差はどんなにかせいでも国内ではたかだか4000mにしかならないが）水平方向に広大な範囲を選べば軌道縞（もしくはTiepointの座標誤差）からの軌道決定感度が増す。これが一般に言われる「広大な領域で干渉解析すべし」の根拠になっている。逆に広大な領域を使ってフラットニングを行えば、（衛星軌道の無視できない２次成分が、画像再生時に悪影響を起こさない限り）ほとんど軌道縞間隔から真の値に近い軌道間距離が決まるはずである。その理由は、例えばBperpが500mならdx/dphi=200m/2πとなり、フルレンジ75kmに渡って300本以上の軌道縞が現れるが、ここに、たとえシーン全体でAzimuth方向の平行フリンジ１本分、すなわち２πの位相変化（片道の視線距離変動に換算して11.8cmの変化）を及ぼす水蒸気ノイズの広域トレンドが被っていたとしても推定される軌道間距離の誤差は致命的にならないと考えられるからである（厳密には、解析領域全体が傾動するような地殻変動があったり、軌道縞そっくりな曲率の水蒸気遅延が被っていたら軌道縞と区別できない可能性は残るし、軌道間距離が特に短いとフリンジ数が極端に少なくなるために本質的な位相誤差の影響が相対的に大きくなるので注意が必要である）。従って、もし軌道縞をフラットにする軌道要素で差分干渉処理をしても顕著な地形縞が残るのであれば、それは水蒸気の鉛直構造の差による見かけの残差である可能性が強い（つまり軌道推定誤差ではなく、水蒸気遅延の本質が映っているものと考えられる）。我々は、アルゴリズムの全く異なる３つの軌道推定法を同じデータセットに適応し、残差地形縞の現われ方がすべて同じであることを確認した。すなわち、いずれの手法でも軌道推定自体は正しく行われており、共通して残った顕著な残差地形縞は、本質的に水蒸気の影響と考えてよいと確信した。この残差地形縞のうち標高に比例した残差位相は、標高に応じて長さが変わる大気伝播経路中の水蒸気量が２つの観測時期で異なるための差の効果でうまく説明できることが多い。それは、後述するように高層気象データを用いて、大気遅延の鉛直プロファイルを計算することで検証できる。それゆえ地殻変動がないと仮定できる領域を選んで、残差地形縞から標高に比例した成分を除去する簡易補正法が一般に行われている。 Fig.1は、差分干渉処理のそれぞれの段階における位相画像を表している。この事例は、93年8月20日と93年11月16日の富士山周辺のJERS-1　SARデータを用いたものである。 (a)はいわゆる軌道縞も地形縞も取り除いていない全く生の干渉位相図である。ここではフルレンジ75kmあたり軌道縞が150サイクル以上現れている。ここで上に書いた様々な精密軌道推定を行い、軌道縞を取り除いたもの（すなわち、干渉処理で作成した、いわゆる“地形縞”）が(c)図である。軌道の推定結果に基づき、DEMを用いて地形縞をシミュレーションしたものが(b)図である。さらに（b）図と(c)図の差を求め、差分干渉図を求めた結果が(d)図である。上で説明したように、一般に標高とよく似た位相画像が得られる。この事例では富士山の標高3776mに対して、残差位相は約20cm（視線距離の見かけの変動量）に及ぶことが分かる。(e)図は、(d)図の中から標高（DEM）と相関した成分だけを抽出したものである。そして(d)図から(e)図を引いた結果が(f)図である。 Fig.1 A series of JERS-1 SAR interferograms around Mt. Fuji (height = 3776m) acquired on August 20, 1993 and November15, 1993 (Bperp=230m) showing (a) raw interferogram, (b) simulated topographic phase, (c) topographic fringes, (d) differential interfrogram without atmospheric correction, (e) 1st order atmospheric phase and (f) atmospheric phase corrected differential interferogram. Fig.2は、富士山の事例と同様に岩手山周辺の結果（98年4月30日と98年7月27日の組み合わせ）である。左図は軌道推定の結果に得られる差分干渉図であり、標高と相関した成分が残存している（Fig.1(d)に対応する）。岩手山頂の標高2038ｍに対して約14cmの残差位相が見られる。ここから標高に相関した成分(Fig.2中図)を差し引いたものがFig.2右図である。 Fig.1やFig.2の一連の位相図から、水蒸気遅延による位相パターンを以下の２つに大別することができる。 パターン１：標高に相関する位相成分（Fig.1(e)、Fig.2 中図） パターン２：標高との相関しない不規則な位相成分（Fig.1(f)） Fig.2 A series of JERS-1 SAR interfrograms around Mt. Iwate (height=2038m) acquired on April 30, 1998 and July 27, 1998 (Bperp=879m) showing raw interferogram ( left figure ), 1st order atmospheric phase ( middle ) and atmospheric phase corrected differential interferogram ( right ). ４．静的な大気構造による電波遅延 標高と比例する位相成分（パターン１）について考えてみる。このような位相成分は一般に、（残差）地形縞と呼ばれているものである。この残差地形縞が「静的な大気（構造）による電波遅延の効果」と見なせるかについて、気象学的に検証してみた。まず観測地域に近いゾンデ観測点の高層気象データを用いて大気遅延量の鉛直プロファイルを求め、パターン１の位相プロファイルと比較した。これは気象学的には、風なし水平一様場を仮定したことになる。 大気遅延Delay(m)は以下の式で表される。ここでは干渉SARの結果に合わせるため、片道分の遅延量に換算してある。 ここでPは気圧(hPa)、Tは気温(K)、eは水蒸気分圧(hPa)で、dlはスラントレンジ方向への積分経路を意味する。右辺第１項が湿潤大気遅延を、第２項が乾燥大気遅延を意味する。 Fig.3は富士山の干渉SARの事例に対応させて、館野のゾンデデータから大気遅延量（の差）を求めた結果である。横軸は標高で、縦軸は大気遅延量Delayの差（93年8月20日の値から11月16日の値を差し引いたもの）である。図では乾燥大気遅延を■で、湿潤大気遅延を▲で、両方を合わせた全大気遅延量を◆で示した。大気遅延の大半は湿潤大気遅延、すなわち水蒸気による電波遅延の効果であることがわかる。富士山山頂の標高3776mにおける大気遅延量の差は約14cmとなった。一方で、差分干渉図（Fig.1(d)）に見られる約20cmの残差位相の大半を説明できる量である。館野観測点が富士山から80kmも離れていることを考えればよく一致しているといえるだろう。 Fig.3 Relationships between height and the difference in atmospheric path delay at two SAR observations in case of Fig.1 (Mt. Fuji). Atmospheric delay is calculated by using radiosonde data at Tateno station. Open triangles show the wet atmospheric delay term, closed squares the dry term and closed diamonds the total mount of each term. Fig.4は岩手山の事例（98年4月30日と98年7月27日）に対応させた、仙台のゾンデデータから大気遅延量の差を求めた結果である。岩手山の標高2038ｍにおける大気遅延量の差は約13cmで、差分干渉SARによる結果とかなりよく一致していることが分かる。 Fig.4 Relationships between height and the difference in atmospheric path delay at two SAR observations in case of Fig.2 (Mt. Iwate). Atmospheric delay is calculated by using radiosonde data at Sendai station. 同様に、 Fig.5に、岩手山の別の事例（98年8月24日と98年10月7日：Ascending軌道からの夜の観測）に対応した仙台のゾンデデータによる計算結果を示す。この場合、干渉SARではほとんど残差地形縞が見られなかった（岩手山で2cm以下）が、ゾンデデータを用いた計算でもたかだか3cm程度であることが見て取れる。 以上の結果から、差分干渉 SARの処理で一般的に見られる（標高に相関するような）残存地形縞は、静的な大気による電波遅延の差で大半が説明できることがわかった。このことは、干渉SARの処理における残差地形縞が、（例えば、軌道推定における）処理上の人工的な誤差なのではなく、本質的に電波遅延の効果を表していることを意味している。一方で、第一近似的には、残差地形縞から標高との比例成分を簡易補正することの妥当性を示している。 Fig.5 Relationships between height and the difference in atmospheric path delay at two SAR observations in case of Fig.2 (Mt. Iwate acquired on August 24, 1998 and October 7, 1998). Atmospheric delay is calculated by using radiosonde data at Sendai station. ５．大気のダイナミックな運動による大気遅延分布 次に、標高と相関のある位相成分を取り除いた後にも見られる不規則な位相パターン２について考える。我々は、複数のデータ組み合わせの差分干渉図を見比べることで、ある特定のデータを含んだ場合に共通して大気遅延量の特徴的なパターンが現れることに経験的に気づいている。そのような複雑な位相パターンは、特定の日における大気のダイナミックな（風の）効果によるものと直観できる。例えば、フェーン現象による、岩手山系の峰を境にした位相トレンド（ 98年8月24日）や伊豆半島上に現れた山岳波による位相の波状構造（98年6月15日）、Fig.11(f)の富士山・箱根周辺の差分干渉図のうち、小田原平野に見られる北西-南東方向に伸びた波状の位相パターン（93年8月20日）などがこの効果である可能性が強い。しかもそのような特定の観測日は、ほとんど例外なく前線が接近、もしくは通過直後であり、特定の方向からの湿った風が山地に向かって吹き続けていた場合が多いことを天気図や地上気象データとの比較により確かめている。そこでそのような特に不均質な水蒸気遅延分布を生じさせたと推定される日の大気構造を数値計算により再現させ、大気遅延量の空間分布を理論計算して差分干渉図と比較検討してみた。 用いた大気モデルは気象研究所モデル (MRI-NHM モデル)で、鉛直方向には27kmを38層に分割し、水平方向には123km四方の領域を1.5kmメッシュに分割し、地形も同じメッシュで与えた。初期条件として、SAR観測領域に近い高層気象観測点のゾンデデータを与え（富士山・箱根・伊豆半島の場合には館野観測点、岩手山の場合には秋田観測点のゾンデデータを用いた）、時間間隔５秒で約3000ステップ（総時間４時間）分の大気運動を計算する。大気運動が定常状態になったところで、SARの視線方向の水蒸気遅延量を次のように計算して空間分布を求め、差分干渉図と比較した。標高と相関した残差地形縞からの偏差に注目するために、大気遅延の水平方向の平均値からの偏差dDelayを求めた。 ここでqvは比湿、ｄは密度、k1,k2は定数、dlはスラントレンジ方向への積分経路を意味し、添字AVGは大気構造モデルの各層で水平平均を取ることを意味している。なお、可降水量1.5mmが（片道の）視線距離変動1cmに相当する。 ６．岩手山周辺のフェーン現象 Fig.6に岩手山周辺の大気運動計算で再現された大気遅延分布を示す。SARデータの組み合わせは98年7月11日と98年8月24日（観測時刻は日本時で22時30分）である。左が差分干渉図、右が大気運動計算による大気遅延分布である。Fig.7に観測前後の天気図を示す。8/24夜に日本海から前線が接近しており、画像領域の北西側に位置する低気圧に向かって湿った南〜南南東の風が岩手山・秋田駒ケ岳に吹きつけていたことが天気図・地上気象データから読み取れる。Fig.6左図で、秋田駒ケ岳南東山麓などで赤色に見えるのは8/24側で水蒸気量の増加があったことを表している。岩手山頂から北側の緑色は大気の乾燥を意味する。右図は、秋田のゾンデデータを用いて大気運動を再現した場合の（衛星方向の）水蒸気遅延量分布である。秋田駒ケ岳の南東山麓に水蒸気が集まっている効果がよく再現されている。秋田駒ケ岳山頂を超えた北西の風下側での大気の乾燥効果も差分干渉SARとよく一致している。少なくともこれらの領域に見られる差分干渉位相のパターンは地殻変動を表しているのではないと判断できるであろう。 Fig.6 JERS-1 differential interferogram around Mt. Iwate (on July 11, 1998 and August 24, 1998 ) corrected for the atmospheric path delay in the first order ( left ) and result of meteorological simulation of atmospheric dynamics using radiosonde data on August 24, 1998 at Akita Station showing the distribution of the atmospheric path delay ( right ). Fig.7 Weather charts at just before and after each SAR observation in the case of Fig.6. ７．伊豆半島を超える山岳波 98年6月15日に伊豆半島で観測された波状の差分干渉位相のパターン（Fig.8左図）についても大気運動の効果を再現してみた。当日の天気図によると、発達中の前線が通過した直後に、さらに次の低気圧が日本海から接近しているという、かなり悪い天候下でSAR観測が行われており、伊豆半島東部の地上でも風速10m/秒以上の湿った南西風が吹き荒れていた。この場合、伊豆半島の達磨山−天城山系に当たった南西風によって、風下側の伊東〜熱海地域に北西−南東方向の軸を持ち上下に振動する山岳波が現れることが大気運動計算でも再現された。Fig.8右図は、計算された高度約3km付近の鉛直流速のコンター図である。SAR位相パターンと高度約3km付近の大気層内の鉛直流速とがよく対応しているのが分かる。すなわち、境界層内の振動により、下降（上昇）流の領域では上層（下層）の乾いた（湿った）空気の進入により大気遅延量が減少（増加）したと説明できる。 Fig.8 JERS-1 differential interferograms around the Izu peninsula (on June 15 ? March 19 and June 15 ? May 2, 1998) corrected for the atmospheric path delay in the first order ( left ) and result of meteorological simulation of atmospheric dynamics using radiosonde data on June 15, 1998 at Tateno showing the vertical velocity at a height of 3.14km ( right ). ８．富士山・箱根周辺のフェーン現象 Fig.9は富士山・箱根周辺における差分干渉SAR図（93年8月20日と93年11月16日）と大気運動計算の結果である。またFig.10に２つのSAR観測における天気図を示す。この事例では、11/16の大気が乾燥していたのに対して、8/20は強い湿った南西風が伊豆半島に吹きつけていた。左図で、富士山や愛鷹山、箱根外輪山の東斜面の風下側で差分位相が青から赤に急変しているのは、8/20に水蒸気遅延がマイナス（11/16に対して減少）に働いていたセンスである。同様に、西側山麓斜面の青から緑の位相変化は、水蒸気遅延がプラスに働いたセンスである。右図は、館野のゾンデデータを用いて大気運動を再現した場合の（衛星方向の）大気遅延量である。差分干渉図の細かいスケールのパターンがかなりよく再現できていることが分かる。例えば、左図では乾燥大気の固まりによる赤色の位相が富士山頂の東から愛鷹山山頂の北東側へ伸びている様がよく再現されている。また、小田原平野に北西-南東方向に軸を持った波状の位相パターンが見られていたが、大気運動計算の結果でも同様なパターンが再現された。これは湿った強い南西風が箱根を吹き抜けた後に形成された山岳波である。 Fig.9 JERS-1 differential interferogram around Mt. Fuji (on August 20, 1993 and November 16, 1993 ) corrected for the atmospheric path delay in the first order ( left ) and result of meteorological simulation of atmospheric dynamics using radiosonde data on August 20, 1993 at Tateno Station showing the distribution of the atmospheric path delay ( right ). Fig.10 Weather charts at each SAR observation in the case of Fig.9. 一方、Fig.11は、湿潤大気遅延の分布と乾燥大気遅延の分布とを分けたものである。大気遅延の大半は、湿潤大気、すなわち水蒸気による電波遅延の効果であることがわかる。 １点のゾンデデータを用いていることを考えれば、水蒸気分布の推定は第一近似としては成功しており、大気遅延によるパターンを区別するにはかなり有効な手法であると言えよう。しかし、定量的には差分干渉SARの大気遅延量と気象計算結果とで３〜５倍程度異なる場合もある。例えば、前述した伊豆半島の山岳波の振幅（98年6月15日の事例）に関して、差分干渉SARと大気運動計算結果の比は６cm：2.8cmであった。 Fig.11 Wet atmospheric delay term (upper) and dry term (lower) in the case of Fig.9 calculated by the meteorological simulation of atmospheric dynamics using radiosonde data. ９．まとめ 衛星軌道を精密に推定することにより、干渉SARに現れる大気遅延による位相パターンを抽出することができた。大気遅延には静的な成分と動的な成分とがあり、それぞれをラジオゾンデデータに基づく「風なし一様場モデル」と「風あり大気のダイナミックモデル」とによって検証した。干渉SARでしばしば現れる顕著な「標高に相関した残差地形縞」は大気遅延の静的な成分(特に水蒸気による湿潤遅延)として定量的に説明できる。このことは逆に、干渉SAR解析における軌道推定法（特に衛星-地球間の幾何学計算、軌道誤差の補正計算）が正しいことをも意味する。山岳波などによる大気遅延分布を予測・再現させることで、干渉SARによる（特に山岳地域の）地殻変動検出の信頼性・精度を向上させることが可能である。 参考文献 [1] Klinklad, H., 1989, A semi-analytical theory for precise single orbit predictions of ERS-1, Mecanique Spatial Space Dynamics, Proceedings of the CNES Symposium, 1989, France, Nov. 1989, pp575-600.