SELENE 計画：月レーダサウンダ（LRS）による月面 SAR / InSAR 観測の計算機シミュレーション
小林 敬生１、小野 高幸１、大家 寛２
１２
東北大
福井工大
１．はじめに現在、我が国では SELENE 計画と名づけられた大型の月探査計画が 2005 年の探査機打ち上げを目指して進められている。SELENE 計画では２トン・クラスの飛翔体に合計 14の科学ミッションが搭載され、月の表面・地下の探査（月の科学）のみならず月軌道周辺のプラズマ、磁場等の宇宙環境計測（月での科学）や月軌道から地球観測（月からの科学）を行なう意欲的な計画となっている。最新の技術を利用して得られる詳細なデータは、月科学や地球科学に飛躍的な進歩をもたらすものと期待されている。
月の科学を行なう SELENE 搭載ミッションのひとつに月レーダサウンダ（Lunar RadarSounder : LRS）がある。LRS は HF 帯電波を利用したレーダであり、SELENE 計画においては月の地下構造を探査し月の誕生・進化の研究に資することを主目的としている。LRSのレーダ方式は FMCW 方式であり、送信波周波数は 200μsec の間に４MHz から６MHzまで直線掃引される。レンジ分解能はパルス圧縮処理後、真空中で 75m である。送受信アンテナは tip-to-tip 30 m のダイポールアンテナで、送信出力は 800W である。SELENE
は月面上空高度 100km の極軌道を周回し、1 年間をかけて月の全球観測をする計画である。
月地下を伝播する電磁波は強い減衰を受けるため観測される地下反射エコーは非常に微弱なものとなる。一方、月表面の大部分は無数の衝突クレータが作り出す複雑な地形を呈しているので、様々な方向、レンジからの表面反射波が微弱な地下反射エコーに重畳して受信される。したがって、LRS による地下探査の成否は、両者の分離・識別を行なうために地下反射エコーの表面反射波に対する S/N をいかに改善するかにかかってくる。我々はこの問題を解決するため、LRS 観測を模擬するシミュレーションコードを開発し、東北大学大型計算機センター（現：情報シナジーセンター）のスーパーコンピュータを利用してLRS 観測データの解析手法の研究開発を進めて、これまでに LRS 観測による月地下探査のための基本的なデータ解析手法を確立した（小林他 1999；小林他 2000）。
地下探査を目的とするデータ解析においては、月表面からの反射波は不要な雑音成分と見なされ、それをいかに取り除きあるいは抑圧して微弱な地下信号を取り出すかが主眼となる。しかし、そうして取り除かれて捨てられる表面反射波のすべてが月表面の情報を含んでいる。したがって、立場を変えて積極的に表面反射波成分の解析を行なえばそこから月表面情報を引き出すことができ、LRS 観測データを無駄なく有効利用することができるはずである。我々はこのような考えから LRS の観測データを利用した月表面探査の手法について検討し、LRS 観測データを用いて月表面の合成開口レーダ（SAR）観測および干渉計合成開口レーダ（InSAR）観測が可能であるとの結論を得た。
本稿では、LRS 観測シミュレーションによって実証された LRS 観測データを用いた月表面の SAR 観測および InSAR 観測、さらに SAR 観測の拡張である２D-SAR 観測について紹介する。

２．LRS 観測シミュレーション
2.1 シミュレーションコード
本稿で用いたシミュレーションコードは、これまでの LRS データ解析手法の研究において開発されてきたもので、LRS パルスの送信から月面における LRS パルスの反射・屈折および地中伝播パルスの地下構造境界面における反射、そして受信される反射波の波形の計算、受信データの解析までを統一的に扱う（Kobayashi 2000）。媒質境界面における LRSパルスの反射・屈折の計算を行なうコードの核となるサブルーチンは Kirchhoff 理論(Beckman and Spizzichino 1963)に基づいて設計された。月表面及び地下反射波データはパルス圧縮処理を受けた後複素 A-scope データとして保存される。この複素 A-scope データをもとに SAR 解析と InSAR 解析が行なわれる。
実観測では一つしかない月を多数のパルスを送信して観測する。シミュレーションでは月面地形を生成した後メモリの許す限り多数のパルスで観測を行なう（１回のランで 400パルス）ようコーディングし、この観測部分を並列化することにより計算時間を短縮した。
2.2 シミュレーションモデル
以下に、本研究で採用したシミュレーションモデルの主な特徴について記述する。
(ⅰ)レーダパルス LRS 送信パルス（5MHz）の電磁場は微小ダイポールアンテナによる放射電磁場でモデル化した。これは、LRS 送受信点が月表面から 100km 離れており、この距離が波長長さ 60m に比べて十分大きいと見なせるからである。
(ⅱ)月面 シミュレーションで用いる月面地形はクレータ形状の統計的性質(Pike 1974)とクレータ径分布の統計的性質(Baldwin 1964)を考慮したモンテ・カルロ法によって数値的に生成している。この時、球面である月表面平均面の曲率は無視し、月平均面は平面であるとした。シミュレーション空間の広がりは 200km×200km とし、観測レンジにシミュレーション月面の境界端がかからないようにしている。モデル
月面の最小クレータ半径は 10m とし、クレータ分布密度は 100[km-2]とした。因みに、月面クレータの飽和分布密度は月面の最小クレータ半径を 10m とすると200[km-2]である（Kobayashi, 2000）。
2.3 シミュレーション条件
主なシミュレーション条件は以下の通り。このシミュレーション条件に基づいて 10 の平行する軌道に沿った観測を行なう。
隣り合う軌道の間隔は 500m あるいは 1000m である。
LRS･……………… 出力：800W
周波数：４−６MHz 直線掃引
周波数掃引時間：200μsec
パルス繰り返し周波数：60Hz
SELENE･………… 高度：月面平均面上空 100km
速度：1600 m/sec
観測飛行距離：40km
月面･………………上層月物質誘電率： 4.0 + i 0.02
クレータ分布密度：100[km-2]
平均（全球）面曲率：考慮せず

３． SELENE LRS 合成開口レーダ（SAR）解析
LRS SAR 解析の原理は、すでに実用化されて久しいマイクロ波 SAR の原理と変わるところはない。しかし、LRS SAR ではダイポールアンテナを送受信アンテナとして用いるため観測によって決定されたターゲット位置に、SELENE 速度方向に関する左右の任意性が生じてしまう。これが、マイクロ波 SAR と LRS SAR の決定的な違いである。この現象は図１にはっきりと現れている。図１は図中黄色の矢印で示される軌道にそって行なわれたLRS 観測データから得られる SAR 画像である。軌道の左右に対称な SAR 画像が作られていることがよくわかる。
ターゲット位置の左右の任意性を取り除くためには複数の軌道による SAR 画像の相関をとる（小林他 2000）。そのため、まず月面平均面上に固定した座標でターゲットの位置を定義し、第ｎ軌道の観測による月面反射源分布を Tn(x,y)とする。この Tn(x,y)にはSELENE 進行方向に関する反射源位置の左右の任意性がまだ残されている。このような
SAR 画像データ Tn(x,y)が各軌道の観測で得られる。次に、これら Tn(x,y)の積相関を
1

 N
N
T ( x, y ) = ∏ Tn ( x, y )
 n =1

に従ってとり、月面反射源の分布を決定する。積相関をとるのは効率よく偽ターゲット信号を抑圧して真ターゲット信号を残すためである。この処理を行なうことで SAR 画像から軌道に関するターゲット位置の左右の任意性が除去される。

図 1 １軌道の観測で得られる LRS SAR 画像。
黄色の矢印が SELENE の軌道を示す。
画像が軌道に関して対称的であることに注意されたい。

図 2 LRS SAR 観測領域のモデル月面の光学画像（上）と LRS SAR 画像（下）。
上図の明度が高く表示された矩形領域は下図の LRS SAR 観測領域に対応する。
下図の SAR 画像は９軌道の観測で得られた SAR 画像の平均である。SAR 解析は原画素の大きさを 40m×40m として行ない、得られた原 SAR 画像に対し 5×5 画素（200m×200m）の平均をとった。図１で見られた対称性が認められないことに注意されたい。

10 の軌道による観測から得られた SAR 画像の積相関処理の結果を図２に示す。比較のため観測領域を含むモデル月面の光学画像を並べて示す。両者の比較を容易にするため、光学画像中、SAR 観測が行なわれる領域は明度を上げて表示している。SAR 解析は原画素の大きさを 40m×40m として行ない、得られた原 SAR 画像に対しスペックルノイズを低減させるため 5×5 画素（200m×200m）の平均をとった。合成開口長は 10km としたので、SELENE 軌道方向の空間分解能は約 300m である。
一見して明らかなように、地球観測衛星の観測で得られている SAR 画像から連想される、光学画像を彷彿とさせるような画像は LRS SAR では得られない。光学画像のクレータ地形との対応を調べると直径がおよそ 10km よりも大きなクレータについては、SELENE 軌道に正対するクレータ内壁部分が明るく見えていることがわかる。白く抜けてしまっている部分は SELENE 軌道直下領域であり、この領域からの表面反射信号の強度は十二分にあるが、信号位相遅れの変化が SAR 画像を構成できるほど十分ではないため画像は判然としない。その他の領域には明るさの違う大小さまざまな光点が散在している。SAR 画像を拡大して、光学画像と比較するとこれらの光点がすべて小クレータの内壁部分に対応していることがわかる。図３に図２の一部拡大図を示す。図３で比較すると、SAR 画像に現れている全ての光点が小クレータに対応し、図では最小直径が数百メートルのクレータとの対応までが認められる。また、カラーコードから明らかなように、大きなクレータほど SAR画像で明るく見えることもわかる。このことから、LRS SAR 画像の明るさを解析することでクレータの分布を推定できることが示される。
SAR 画像の明るさから地表面粗度を推定することは地球観測分野では広く行なわれているが、地球では地中の水分の量が SAR 画像の明るさに大きく影響を与えるため推定精度を上げることは難しい。一方、月面には水分はないので、LRS SAR ではこのような問題を心配する必要はない。このことから、LRS SAR 観測によって、月面の地質年代決定に重要な情報であるクレータ分布が比較的容易に得られることが期待される。

図３ 光学画像（左）と LRS SAR 画像（右）の比較。
図２の一部を拡大図して表示している。SAR 画像に見られる光点すべてがクレータに対応していることに注目されたい。

４．SELENE LRS 干渉計合成開口レーダ（InSAR）解析
複数の軌道による SAR 解析複素画像データが得られると InSAR 解析が可能となる。
InSAR 解析は異なる軌道の観測で得られる同一地点からの反射波の位相遅れの差がその地点の高度と１対１対応することを利用して、表面地形の高度分布を算出する解析法である（Graham, 1974）。
本稿では、まず独立した 10 の軌道から適当に 13 組の軌道の対（軌道間隔 500m、1000m、1500m）を選び、それぞれの組について上述した SAR 解析同様 40m×40m を原画素の大きさとして InSAR 解析を行なった。次いで、各軌道組の InSAR 解析結果に対し 65×65画素（2600m×2600m）の移動平均をとった後、９つの InSAR 解析結果の平均をとった。
このように平均操作を繰り返したのは、以下の理由による。
まず、組となる各軌道の観測で得られる SAR 画像要素複素データには、着目する地形表面からの反射波だけではなく、同レンジの他地点からの反射波データが重畳している。これは、LRS の送受信アンテナがダイポールアンテナであるためである。このため、InSAR解析で用いられる（第ｎ組の）２軌道間の観測データ位相差φ n には、期待される位相差φ 0に重畳波成分による位相差 ∆φ n が加わる。すなわち、
φ n = φ 0,n + ∆φ n 。
軌道の選び方と観測される表面地形との間には何らの相関も存在しないので重畳波成分の位相差 ∆φ n は一様ランダムな分布を示し、さらに、位相差φを − π ≤ φ ≤ π で定義すれば、∆φ n の平均値の期待値は０となる。さて、この第ｎ組軌道の観測位相差から求まる表面地形高さ hn には φ 0,n と ∆φ n に対応して、真の高さ h0 と高さ誤差 ∆hn が含まれることになる。
∆φ n と ∆hn とは線形関係にはないため ∆hn の分布は一様分布とはならないが、近似的には一様ランダム分布を示すものとして扱ってよい。したがって、観測数 N が十分大きければ ∆hnの平均値の期待値は０となるはずである。すなわち、
1 N
1 Nhn = h0 + ∑ ∆hn → h0∑
N n =1
N n =1となり、着目地点の高さ h0 が求まるのである
図４に InSAR 解析の結果をモデル月面の高さ分布図とともに示す。解析領域は前述のSAR 解析を行なった領域と同じである。
さて、両図の比較から明らかなように、LRS InSAR 解析の結果では数 km スケール以上の大きな地形がよく再現されている。しかし、一方で一部に不可解な地形が出現している。
まず、図４（下図）中央部に存在する深く長い巨大渓谷は SELENE 軌道直下領域に相当し、前述したように、受信信号の十分な位相変化が観測されない領域である。このような領域では InSAR 解析を行なうことも不可能となるために、この図に見られるような実在しない地形が現れる。さらに、直径が 10km を超えるいくつかのクレータで不自然な地形高度分布が見られるが、これは phase wrapping と呼ばれる現象によるもので、地形高度と位相差は 1 対 1 対応するが、実際に観測される位相は − π ≤ φ ≤ π に限定されて 2πの整数倍差を認識できないことによる。この問題の解はすでに得られており phase unwrappingの方法としてマイクロ波 InSAR 観測で実用化されている（Massonnet and Feigl, 1998）。
次に前項の SAR 解析結果で見た図３と同じ領域の InSAR 解析結果（図５）を見てみる。

図４ モデル月面標高分布図（上）と LRS InSAR 解析結果（下）。
図２と同様、上図の明るい矩形領域は下図の InSAR 解析を行なった領域に対応している。下図は 40m×40m を原画素の大きさとして行なわれた InSAR 解析結果それぞれについて 65×65 画素（2600m×2600m）の移動平均をとり、９つの InSAR 解析結果の平均をとったものである。

図では InSAR 解析で得られた地形標高分布にモデル月面の光学画像を重ね実地形との対応を見易くした。この InSAR 解析では、13 の軌道対のデータから得られた原 InSAR 画像のそれぞれに対し 17×17 画素（680m×680m）の移動平均をとり、13 の InSAR 画像の平均をとった。図には比較のため、モデル月面の地形高度分布を同様にして示したが、両図の比較から明らかなよう LRS InSAR 解析では確かに数 km オーダー以上の地形がよく再現されている。実際、両図のずれの大きさを調べてみると、地形高度分布の誤差の標準偏差は 112m であった。
以上から、LRS InSAR 解析は LRS 送信波波長の数倍程度の誤差、km オーダーの空間分解能で月面の地形高度分布を再現できる能力のあることが結論される。
ところで、SELENE 計画では月面地形の精密な測量を目的としたステレオカメラによる地形計測を行なうミッションも搭載される。地形測量精度のみを問題とする場合、光学カメラと LRS InSAR の能力の差は歴然としており、LRS InSAR の積極的意義はその実験的試み以外にはないように考えられるかもしれない。しかし、HF 帯電波の地中透過力を考えると LRS 観測によって InSAR 解析を行なう意義が明らかとなる。
表面粗さが小さい領域に斜め入射する電波の表面反射波のエネルギーは大部分が鏡面反射成分で、LRS で受信される後方散乱成分は非常に小さい。一方、表面で屈折して地下へ伝播する成分は地下境界面で正面反射条件を満足して表面の後方散乱波に比べ十分な強さをもって受信される可能性がある。そのような場合には、地下反射波成分が表面反射波成分に卓越して、InSAR 解析で決定される地形高度は地下境界面地形の高度となる。このような条件が満足される場所は 海 として月面に広がっている。すなわち、LRS InSAR 解析では月地下境界面地形を直接 見る ことができる可能性を持っており、ここに LRSInSAR 解析を行なう意義が存在するのである。

図５ モデル地形標高図（左）と LRS InSAR 解析結果の比較。
図３と同じ領域を拡大して示している。右図は原 InSAR 画像に対し 17×17 画素（680m×680m）の移動平均をとり、さらに 13 の InSAR 画像の平均をとったものである。

５．月極域における２次元合成開口レーダ（２D-SAR）観測SELENE 衛星は極軌道を運航するため極域上空では軌道の密度、ひいては LRS 観測数密度が非常に大きくなる。観測点の空間密度、観測数ともに十分そろえば SAR 解析アルゴ
リズムを拡張することにより、データ空間上で２次元開口アンテナを合成（２D-SAR）して、表面地形の撮像と地形高度分布の計測が可能となる。この拡張は、SAR 解析アルゴリ
ズムで言うところのアジマス圧縮処理が、実は受信信号位相遅れの理論値と観測値の相関処理と同値であることから極めて自然に行なわれる。
今、月表面上のターゲットの位置を R 0 とし、 n 番目の観測における R 0 からの反射波の位相遅れの理論値を ϕ theo,n (R 0 ) 、同じく n 番目の観測の複素 A-scope データにおいて R 0 が属するレンジビンの観測値を An exp{ iϕ obs,n (R 0 )}とすると、
I (R 0 ) =∑ A exp[iϕnobs , n(R 0 )]exp[− iϕ theo,n (R 0 )]nは、SAR 画像におけるターゲットの明るさとなる。この式は観測点相互の幾何学的配置には何の制限も要求しない。仮に、観測点が一直線上に並べば、この式は SAR 解析におけるアジマス圧縮処理そのものを意味する。
次に R 0 を
R 0 ≡ R 0 ( z ) = ( x0 , y 0 , z )として、 x0 、 y 0 を固定して z （ターゲット標高）を変化させると、 I (R 0 ) （ターゲットの明るさ）が変化するが、式の形から明らかなように、 I (R 0 ) は∀n ϕ obs,n (R 0 ( z ) ) = ϕ theo,n (R 0 ( z ) )のとき最大値を示す。つまり、z0 ≡ { z0 | I ( x0 , y0 , z0 ) = max[I ( x0 , y0 , z )] }なる z 0 が実ターゲットの標高となる。
この２D-SAR アルゴリズムにしたがってデータを処理したシミュレーション結果を図６と図７に示す。観測点は − 10 ≤ x ≤ +6 km、 − 16 ≤ y ≤ 0 km の領域上空 100km で x 方向、y 方向ともに 40m 間隔に設定し、合計 401 × 401 = 160801 点である。図７の標高分布を得るに当たっては、 z の関数 I (x0 , y 0 , z ) に対してウエーブレット分解アルゴリズムを適用してI (x0 , y 0 , z ) の平滑化を行なってから z 0 を求めた。用いたウエーブレットは４階のスプラインウエーブレットで、分解レベルは−2 である。
図６では前述した LRS SAR 画像と違い、大きなクレータの丸い輪郭がくっきりと浮かび上がり直径が数 km 以上のクレータを数えることもできる。図７では２D-SAR 解析でモデル月面の標高分布がよく再現されていることがわかる。図中、右上の領域では標高分布がモデルとかけ離れた値を示しているが、これは解析領域が２D-SAR アンテナ直下点から遠く離れて相関強度が小さくなり真値を見出すことができなくなったためである。
この２D-SAR シミュレーションでは 16 万点を超える LRS 観測のシミュレーションを行なわなければならない。クレータ地形計算コードが完成した当初は、コードの並列化がなされておらず、本稿で使われた月面モデルで LRS の観測シミュレーションを行なうと１パルスの観測シミュレーションに約１時間の CPU 時間がかかった。したがって 16 万点のLRS 観測シミュレーションを行なうためには 18 年間の CPU 時間が必要となる計算になる。
今回のシミュレーションでは、コードの並列化を進めるとともに冗長な計算ループを削除

図６

図７モデル月面光学画像（左）と２D-SAR 解析画像（右）。

モデル月面標高図（左）と２D-SAR 解析によって得られた標高図。

するなどしてコードの最適化を行なった。そして、利用者が極端に少なくなる年度末の時期に 32CPU 並列計算のジョブを複数同時に走らせることで、１パルスの観測シミュレーションに要する時間を実時間で約 10 秒とするまでに短縮し、１ヶ月余りですべてのシミュレーションを終えることができた。18 年と 1 ヶ月の差は大きい。並列計算のできるスーパーコンピュータのありがたみが強く実感させられた。現在はこの 16 万パルスを越えるシミュレーション結果の詳しい解析を進めているところである。
まとめ
2004 年に打ち上げを予定する月探査計画 SELENE に搭載される HF 帯月地下探査レーダ LRS について、データ解析手法の確立を目的として観測からデータ解析までを含むシミュレーションを行なった。その結果、これまでに前例のないダイポールアンテナを用いた
HF 帯 SAR 解析、InSAR 解析が可能であることを実証しさらに、極域における表面地形撮像・計測が２D-SAR アルゴリズムにより可能であることを示した。
謝辞
本研究は、東北大学大型計算機センターの平成 12 年度共同研究、「HF 帯レーダによる月探査における月面表層での電波の反射透過特性に関する計算機シミュレーション」として行なわれた。シミュレーションコードの並列化・最適化にあたっては、大型計算機センター（現：情報シナジーセンター）から、適切かつ有益な助言を頂いた。この場を借りて感謝したい。
参考文献
Baldwin, R. A., Lunar crater counts The Astronomical Journal vol. 69 no. 5377-392, 1964
Beckman and Spizzichino., The Scattering of Electromagnetic Waves from RoughSurfaces, Oxford, Pergamon, 1963
Graham, L. C., Synthetic interferometer radar for topographic mapping, Proc. IEEE,62, 763-768, 1974
Kobayashi, T., Computer simulation on investigation of lunar subsurface structure byradar sounders ---- studies related to the SELENE project, Doctoral dissertation,
Tohoku University, 2000
Massonnet, D. and Kurt L. Feigl, Radar interferometry and its application to changes
in the earth’s surface, Rev. Geophys., 36, 441-500, 1998
Pike, R. J., Depth/diameter relations of fresh lunar craters: revision from spacecraft
data, Geophysical Research Letters vol. 1, no. 7, 1974
小林敬生、小野高幸、大家寛 SELENE 計画 月レーダサウンダによる月地下構造探査の計算機シミュレーション、SENAC,Vol.32-2、1999
小林敬生、小野高幸、大家寛 SELENE 計画 月レーダサウンダによる月高地領域地下構造探査と合成開口レーダ解析の計算機シミュレーション、SENAC,Vol.33-3、2000