合成開口レーダーシステムによって取得された入力信号を処理する方法
【課題】本方法は、アジマスアンビギュイティ（azimuth ambiguity：アジマス曖昧性）を制御し、画像モデルを組み込んで飽和したＳＡＲデータからの再構成を可能にする。
【解決手段】合成開口レーダー（ＳＡＲ）システムによって取得された飽和入力信号が、該入力信号を生成した再構成を推定し、推定された再構成から入力信号を再生して再生信号を生成し、該再生信号を入力信号と比較し、推定された再構成を該比較に基づいて調整し、終了条件に達するまで再生するステップから反復することによって処理される。
【解決手段】合成開口レーダー（ＳＡＲ）システムによって取得された飽和入力信号が、該入力信号を生成した再構成を推定し、推定された再構成から入力信号を再生して再生信号を生成し、該再生信号を入力信号と比較し、推定された再構成を該比較に基づいて調整し、終了条件に達するまで再生するステップから反復することによって処理される。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、包括的には合成開口レーダーに関し、より詳細には飽和データから高品質画像を形成することに関する。
【背景技術】
【０００２】
逆問題
逆問題の研究は、信号処理システムおよび合成開口レーダー（ＳＡＲ）システムを含む種々の分野において、いくつかの用途を有する十分確立された学問分野である。逆問題の理論は、関数を逆にして、該関数の出力から入力を推測することができる場合、および方法を調べる。そのような問題は、信号取得システムおよび信号処理システムにおいて、豊富に存在する。これらの問題を逆問題として定式化することによって、画像のノイズ除去、ぶれ除去、および超解像の分野において大きな成果が生み出されてきた。
【０００３】
信号処理における一般的な逆問題は、測定値ｙの集合からの信号ｘの復元である。
【数１】
ここで、関数Ａは、取得システムをモデリングしており、ｎは、ノイズを表す。再構成プロセスの目標は、信号ｘに近い信号推定値
【数２】
を求めることである。
【０００４】
取得関数Ａが可逆であり、かつノイズが無視できる場合、関数Ａの逆関数を用いてｘを求めることが自明の選択肢である。しかしながら、この方法は、Ａが単射でない場合、すなわち、データを複数の信号で説明することができる場合、またはかなりのノイズが存在する場合に、失敗する可能性がある。より一般的な手法は、以下のコスト関数を最小化することによって、ｘを推定することである。
【数３】
ここで、ｆ（・，・）は、取得システムおよびノイズの特性に従ってデータ忠実度を測定するコスト関数であり、ｇ（ｘ）は、対象の信号の特性に関する知識を組み込んで望ましくない解にペナルティーを課す正則化項であり、λは、コスト関数内の２つの項間のトレードオフを制御する。式（２）における定式化によって、取得シナリオおよび信号モデルのレンジに適応する際の柔軟性がもたらされる。
【０００５】
ｆ（ｙ，Ａ（ｘ））＝‖ｙ−Ａ（ｘ）‖２、ｇ（ｘ）＝０であり、関数Ａが線形であり、かつシステムが過剰決定である特殊な場合、式（２）の解法は、疑似逆関数Ａ†を用いる。疑似逆関数は、システムが劣決定であり、かつｇ（ｘ）＝‖ｘ‖２である、すなわち、最小エネルギー解が所望されるときにも解となる。
【０００６】
ＳＡＲ画像形成
ＳＡＲ画像の形成は、逆問題として定式化することができる。特に、理想的なＳＡＲ取得プロセスは、線形系、すなわち、式（１）の実例とみなすことができる。ここで、ｘは、表面反射率の２次元（２−Ｄ）画像を表し、ｙは、受信信号を表し、関数Ａは、線形である。
【０００７】
これ以降、Ａを用いて線形取得関数を表す。受信信号ｙは、通常、ＳＡＲ画像に類似した２−Ｄ形式で構成される。各反射パルスのサンプルは、行ベクトルを形成するとみなすことができ、異なる反射パルスが一緒にスタックされて、データの行列を形成する。関数Ａが線形であることを強調するために、表記Ａｘを用いてグラウンド信号（画像）ｘへの関数の適用Ａ（・）を表す。ｘを２Ｄ画像の代わりに１つの単一列ベクトルに再構成すると、Ａｘは、行列乗算とみなすことができる。
【０００８】
データ行列ｙの各行は、取得プラットフォームの経路に沿った該取得プラットフォームの異なる位置に対応する。経路に沿った次元は、アジマスと呼ばれる。各列は、パルスの送信からの遅延に対応する。遅延は、パルスが進んだ距離に比例するので、この第２の次元は、レンジと呼ばれる。
【０００９】
逆問題の観点から、いくつかの既知の画像形成手順は、疑似逆関数Ａ†の近似を求めることとして解釈することができる。近似によって、ＳＡＲ画像ｘの推定を効率的に実施することが可能になる。例えば、多くの手順が、パルス圧縮に依存する。パルス圧縮は、送信パルスとの相関を通じた受信信号の近似逆畳み込みを指す。相関は、高速フーリエ変換（ＦＦＴ）を用いることによって非常に効率的にすることができる。パルス圧縮は、受信パルスがナイキストレートでサンプリングされている場合、少なくとも対象の帯域幅内で厳密な逆関数に対応し、受信パルスが僅かにオーバーサンプリングされている場合（これは、実際のＳＡＲシステムにおいて通例である）、疑似逆関数に密に対応する。
【００１０】
最も広く用いられている画像形成手順のうちの１つは、チャープスケーリングアルゴリズム（ＣＳＡ）である（参照により本明細書に援用される、Runge他「A novel high precision SAR focusing algorithm based on chirp scaling」（International Geoscience and Remote Sensing Symposium（IGARSS），May 1992，pp. 372-375））。
【００１１】
図１に示すように、入力信号１０１にアジマスＦＦＴ１１０が適用され、続いて、チャープスケーリング１２０、レンジＦＦＴ１３０、レンジ圧縮１４０、レンジＩＦＦＴ１５０、チャープスケーリング１６０、アジマス圧縮１７０、およびアジマス逆ＦＦＴ（ＩＦＦＴ）１８０が適用され、出力信号１０９が生成される。
【００１２】
チャープスケーリングは、線形チャープ信号の時間変動遅延を、２つの他のチャープ信号と乗算することによって概算する。ＣＳＡでは、取得プラットフォームの動きによって引き起こされたレンジマイグレーション、所与のターゲットへの距離の変動、したがって関連する遅延を補正するのにチャープスケーリングが用いられる。ＣＳＡは、その組成に起因して、ＦＦＴおよび乗算に関してのみ、非常に効率的である。この理由により、ＣＳＡは、後述するように、取得関数およびその随伴関数を求める本方法の基礎を成す。
【００１３】
飽和および逆問題
データ振幅における大きな変動および精度の低い量子化器の使用に起因して、データ値の飽和は、ＳＡＲ信号取得において大きな問題である。ＣＳＡ等の従来の画像形成方法は、通常、データ飽和に明確に対処しない。
【００１４】
信号飽和は、アナログ／デジタル（Ａ／Ｄ）変換システムにおいて遭遇する、非常に一般的な問題である。なぜなら、電子部品は、有限の電圧レンジを有し、信号振幅は、それを超えて変動してはならないためである。加えて、量子化器は、有限数の量子化レベルを用いる。取得した信号振幅が飽和閾値±Ｔを超えて変動する場合、該取得した振幅は、±Ｔに飽和する。
【００１５】
飽和自体は、望ましくないが、厳しい量子化の存在下では、飽和は、有益となり得る。特に、信号の利得を増大させて飽和を促すことによって、非飽和測定値における信号対量子化雑音比も増大する。飽和を適切に考慮に入れた場合、全体再構成性能を増大させることができる。
【発明の概要】
【発明が解決しようとする課題】
【００１６】
飽和に対処する最も一般的な手法、すなわち、飽和が発生しなかったかのように再構成において飽和値をそのまま用いることによって、重大なアーティファクトが生じる。再構成誤差は、一貫した再構成手法、すなわち、再構成された信号推定値が再取得された場合に、同じ飽和を生じることを確保することを用いて低減することができる。一貫した再構成は、量子化され飽和した測定値とともに用いられた場合に、大幅に再構成誤差を改善する。
【課題を解決するための手段】
【００１７】
本発明の実施の形態によれば、合成開口レーダー（ＳＡＲ）画像の形成が逆問題として定式化される。本発明の実施の形態は、データの飽和が存在するとき、特に一方向二次コスト関数を最適化して受信データとの一貫性をもたらすことによって、ＳＡＲ画像の品質を増大させる。ＳＡＲ取得プロセスは、線形関数を用いてモデリングされる。関数は、実施するのに効率的であり、反復最適化プロセスにおいて、該関数の随伴関数とともに用いることができる。本方法は、アジマスアンビギュイティを制御し、画像モデルを組み込んで飽和したＳＡＲデータからの再構成を可能にする。
【発明の効果】
【００１８】
本発明の実施の形態は、合成開口レーダー（ＳＡＲ）画像の再構成を逆問題として定式化する。ＳＡＲは、リモートセンシングにおいて、比較的小さなアンテナを用いて高解像度で広域のレーダー画像を得るのに広く用いられている。高解像度は、通常、平面アンテナまたは衛星アンテナであるアンテナを搭載するプラットフォームの動きを利用して、はるかに大きな開口を合成することによって可能になる。ＳＡＲにおける取得プロセスは、プラットフォームが動くときに定期的な間隔でパルスを送信し、反射を記録して下方の表面の反射率に関する情報を得ることを伴う。
【図面の簡単な説明】
【００１９】
【図１】従来のチャープスケーリングアルゴリズムのブロック図である。
【図２】本発明の実施の形態による合成開口レーダー画像形成に用いられる取得関数のブロック図である。
【図３】本発明の実施の形態によるＳＡＲ画像を再構成する方法の流れ図である。
【発明を実施するための形態】
【００２０】
有限の動的レンジを有する線形系を
【数４】
によって表す。ここで、Ａは、システム取得関数であり、Ｓ（・）は、要素単位で入力信号の全ての成分に適用された非線形スカラー飽和関数
【数５】
である。
【００２１】
従来の最小二乗忠実度基準の代わりに、飽和測定値に対し一方的な二次ペナルティを強制し、一貫した画像再構成を確保する。より具体的には、二次ペナルティについて、コスト関数
【数６】
を用いる。ここで、（ｙ）−＝ｍｉｎ（ｙ，０）であり、（ｙ）＋＝ｍａｘ（ｙ，０）であり、ｉは、データの実数部および虚数部を別個にインデックス付けする。利点として、コスト関数は、一貫性のない解にペナルティを課し、最適化が単純であり、その二次特性に起因して飽和が発生する前に加法性雑音を低減する。
【００２２】
式（５）において、飽和にロバストなコスト関数を最小化することによって、多くの場合に受入可能な解がもたらされるが、信号モデルも用いて、データの飽和および量子化に固有のアンビギュイティをさらに解消することができる。
【００２３】
ＳＡＲ画像形成
ＡおよびＡＨの実施
式（５）のコスト関数は、種々の反復最適化手順を用いて最小化することができる。ほとんどの手順は、取得関数Ａおよびその随伴関数ＡＨを繰返し計算することを必要とする。これらの関数の効率的な実施は、画像の高次元性に起因して、ＳＡＲ画像再構成において、特に重要である。
【００２４】
効率性を達成するために、本実施の形態は、上述したチャープスケーリングアルゴリズム（ＣＳＡ）に大まかに基づく。ＣＳＡは、疑似逆関数Ａ†の近似であることから、未加工のＳＡＲデータ２０９を再生するのに用いることができる取得関数Ａは、図２に示すブロック図を有する。これは、本質的に従来のＣＳＡの関数の反対である。
【００２５】
入力信号２０１にアジマス高速フーリエ変換（ＦＦＴ）２１０が適用され、続いてアジマス周波数応答Ｈａ（Ｒ０，ｆη）の乗算２２０、第１のチャープスケーリング２３０、レンジＦＦＴ２４０、レンジ周波数応答Ｐ２ｄｆ（ｆτ，ｆη）の乗算２５０、レンジ逆ＦＦＴ（ＩＦＦＴ）２６０、第２のチャープスケーリング２７０、およびアジマスＩＦＦＴ２８０が適用され、未加工ＳＡＲデータである出力信号２０９が生成される。
【００２６】
上記の関数は、当該技術分野において既知のメモリおよび入力／出力インターフェースに接続されたプロセッサ２００、またはディスクリートハードウェア構成要素において実施することができることを理解されたい。
【００２７】
ＳＡＲ画像を生成する随伴関数を実施するには、関数を合成したものである随伴関数が、個々の随伴関数を逆の順序で合成したものであるということを利用する。ＦＦＴの随伴関数は、スケール係数次第でＩＦＦＴであり、関数の乗算の随伴関数は、同じ関数の複素共役の乗算である。このため、随伴関数は、アジマス圧縮関数がアジマス周波数応答
【数７】
と置き換えられた、ＣＳＡと同じブロック図を有する。ここで、＊は、複素共役を示し、ａは、アジマス応答を意味し、Ｒ０は、最接近のレンジを示し、ｆηは、アジマス周波数である。
【００２８】
このとき、レンジ圧縮関数は、フーリエ変換
【数８】
であり、２ｄｆは、２次元周波数領域を意味し、ｆτは、レンジ周波数である。
【００２９】
次に、本発明者らの取得関数Ａの実施を、より詳細に説明する。画像が点目標物の長方形格子を含むと仮定する。単位振幅の単一点目標物による応答に焦点を当てる。関数Ａは、アジマスにおいてシフト不変であるので、レンジドップラー（レンジ時間τ、アジマス周波数ｆη）領域における応答を、等しく特定することができる。最接近のレンジＲ０における単位点目標物の場合、レンジドップラー応答は、
【数９】
によって近似的に与えられる。ここで、ｒｄは、レンジドップラー領域を意味し、Ｈａ（Ｒ０，ｆη）は、アジマス周波数応答であり、ｐ（τ，ｆη）は、送信パルスであり、Ｄ（ｆη）は、レンジマイグレーション係数であり、ｃは、光速である。パルスｐ（τ，ｆη）は、線形チャープであると仮定され、レンジドップラー領域において、ｆηの関数である。
【００３０】
式（７）における応答は、レンジマイグレーション係数Ｄ（ｆη）の存在に起因して、レンジがシフト不変でない。しかしながら、線形チャープパルスの場合、変更されたチャープを用いた畳み込み、その前後におけるチャープスケーリング関数との乗算によって、チャープスケーリングを用いて、式（７）におけるレンジ応答を近似することができる。図２に示すように、畳み込みは、ＦＦＴを用いてフーリエ変換Ｐ２ｄｆ（ｆτ，ｆη）との乗算として実施される。
【００３１】
図２は、取得関数Ａを計算する複雑度が図１のＣＳＡを計算する複雑度と同じであることを示している。より高速な計算のために、信号と乗算された図２の関数の全てを予め求め、格納することができる。離散時間型の実施では、Ｈａ（Ｒ０，ｆη）およびＰ２ｄｆ（ｆτ，ｆη）を最小のエイリアシングで表すのに十分高いサンプリングレートを用いるように注意すべきである。
【００３２】
基本最適化
図３に示すように、一次勾配降下プロセスを用いて、式（５）のコスト関数を最小にする。より洗練された二次手順も検討することができる。この最適化は、初期構成ｘを推定すること（３１０）から始まる。初期再構成ｘは、例えば、ＣＳＡを用いて計算することもできるし、０に設定することもできるし、任意に設定することもできる。
【００３３】
各反復３６０中に、
【数１０】
を用いて信号ｙが再生され（３２０）、再生信号は、入力信号と比較され（３３０）、式（６）において残差ベクトルｒが求められ、式（５）において二次ペナルティを用いて総コストが求められる。推定値を調整する（３４０）ために、
【数１１】
を用いて、式（５）に対応する勾配ｇが計算される。ここで、ＡおよびＡＨは、上述したように求められる。次に、サイズτのステップが勾配ｇと反対方向に取られる。すなわち、反復は、終了条件３５０に達するまで
【数１２】
として更新される。
【００３４】
厳密な線探索を通じて、負の勾配−ｇに対して平行な線に沿って、関数ｆ（ｙ，Ａｘ）を最小化するステップサイズτが最適に選択される。この線探索は、１−Ｄ凸区分的二次関数を効率的に最適化する。
【００３５】
いくつかの可能な終了条件が存在する。例えば、再構成された画像推定値ｘの変化が小さい場合、または比較ステップにおける式（５）の二次ペナルティが小さい場合、反復は、終了する場合がある。本方法は、終了条件に達したときに、再構成画像を出力する（３７０）。上記のステップは、当該技術分野において、既知のメモリおよび入力／出力インターフェースに接続されたプロセッサ３００において実行することができる。
【００３６】
アジマスアンビギュイティおよびモデル不整合
アンダーサンプリングに起因したアジマス周波数領域におけるアンビギュイティは、多くのＳＡＲ画像形成手順に影響を与える問題である。パルス繰返し周波数に等しいアジマスサンプリングレートは、通常、アジマス応答Ｈａ（Ｒ０，ｆη）の帯域幅に適応するには不十分であり、画像化領域の幅を減らさずに、増大させることができない。アンダーサンプリングの結果、エイリアシングが生じ、エイリアシングは、ゴースト像が特定のアジマスおよびレンジオフセットにおいて重なる形態で現れる。エイリアシングからの復元は、画像に関する追加情報なしでは困難である。例えば、画像自体が帯域制限されていると仮定することはできない。
【００３７】
上述した基本方法では、真の取得システムと、再構成において用いられるモデルとの間の不整合によって、さらなる難題が提起される。例えば、アジマス応答Ｈａ（Ｒ０，ｆη）は、アンテナの長さおよび重み付けパターン等のパラメータに依拠する。これらのパラメータの値を精密に知ることができず、不正確な結果につながる可能性がある。特に、本発明者らの経験によれば、モデル不整合は、アジマスアンビギュイティによって生じたゴースト像を大きくさせる傾向にあることが示唆される。ゴースト発生を低減する１つの方法は、データを非適応的にフィルタリングするか、またはエイリアシングによって最も影響を受けるアジマス周波数を、空間において適応的に除去することである。しかしながら、その手法は、再構成に利用可能なデータ量、したがって、結果の画像の解像度を低減する可能性がある。
【００３８】
代替的な手法は、上記の基本方法における各反復の終了時に、ゴースト除去を実行することである。ＳＡＲパラメータの知識に基づいて、対応する真の目標物に対するゴーストのロケーションを予測する能力を利用する。本方法の第１の反復において、まず、低解像度でアンビギュイティのない画像を形成し、次に、ゴーストの距離特性で隔てられた領域間の局所電力を比較することによって、主要なゴースト源である可能性が高い画像領域が特定される。後続の反復において、現在の画像を１度に１つのアジマス線ずつ該現在の画像自体と相関させ、次に、推定値を画像から減算することによって、主要なゴースト源によるゴーストが推定される。推定および減算を主要なゴースト源に限定することによって、ゴーストの効率的な除去が可能になる。
【００３９】
逆問題との関連で、本方法は、特に、ある特定のロケーションおよびシフトにおける画像の自己相関において、式（２）の再構成最適化問題に対し、明示的な制約を付加することとみなすことができる。これは、大きな自己相関値は、ゴーストのない画像内に生じるはずがないという前提に基づく。代替の形態は、そのような知識を、信号モデル項ｇ（ｘ）に組み込む。
【００４０】
信号モデル
再構成正確度における大きな改善は、再構成に信号モデルを組み込むことによって達成することができる。上述したゴースト低減モデルに加えて、自然な画像は、大きな構造を呈する。本発明者らが課すモデルは、ＳＡＲ画像の大きさのウェーブレット変換が、スパースであるというものである。本明細書で定義されるとき、スパース信号は、複数の非ゼロ値の数が、ゼロ値の数よりも大幅に小さい場合に１である。
【００４１】
ｌ１最小化の固定点継続手順に類似した方式のソフト閾値処理を用いて、モデルを課す。参照により本明細書に援用されるHale他「Fixed-point continuation for l1-minimization: methodology and convergence」（SIAM J. Optim., vol. 19, no. 3, pp. 1107-1130, Oct. 2008）を参照されたい。
【００４２】
本方法の全ての反復の終了時に、画像またはその大きさを変換し、ソフト閾値処理を実行し、逆変換する。画像の大きさのみが変換された場合、同じ位相を維持する。このモデルの動機付けは、ＳＡＲ画像の大きさが、自然な画像と類似して見え、したがって、ウェーブレット変換を用いて圧縮可能であるということである。一方、画像化されている特徴の高さの変動に依拠し、深刻なラッピングを受ける画像の位相は、再構成において利用することができる構造を提供しないようである。
【技術分野】
【０００１】
本発明は、包括的には合成開口レーダーに関し、より詳細には飽和データから高品質画像を形成することに関する。
【背景技術】
【０００２】
逆問題
逆問題の研究は、信号処理システムおよび合成開口レーダー（ＳＡＲ）システムを含む種々の分野において、いくつかの用途を有する十分確立された学問分野である。逆問題の理論は、関数を逆にして、該関数の出力から入力を推測することができる場合、および方法を調べる。そのような問題は、信号取得システムおよび信号処理システムにおいて、豊富に存在する。これらの問題を逆問題として定式化することによって、画像のノイズ除去、ぶれ除去、および超解像の分野において大きな成果が生み出されてきた。
【０００３】
信号処理における一般的な逆問題は、測定値ｙの集合からの信号ｘの復元である。
【数１】
ここで、関数Ａは、取得システムをモデリングしており、ｎは、ノイズを表す。再構成プロセスの目標は、信号ｘに近い信号推定値
【数２】
を求めることである。
【０００４】
取得関数Ａが可逆であり、かつノイズが無視できる場合、関数Ａの逆関数を用いてｘを求めることが自明の選択肢である。しかしながら、この方法は、Ａが単射でない場合、すなわち、データを複数の信号で説明することができる場合、またはかなりのノイズが存在する場合に、失敗する可能性がある。より一般的な手法は、以下のコスト関数を最小化することによって、ｘを推定することである。
【数３】
ここで、ｆ（・，・）は、取得システムおよびノイズの特性に従ってデータ忠実度を測定するコスト関数であり、ｇ（ｘ）は、対象の信号の特性に関する知識を組み込んで望ましくない解にペナルティーを課す正則化項であり、λは、コスト関数内の２つの項間のトレードオフを制御する。式（２）における定式化によって、取得シナリオおよび信号モデルのレンジに適応する際の柔軟性がもたらされる。
【０００５】
ｆ（ｙ，Ａ（ｘ））＝‖ｙ−Ａ（ｘ）‖２、ｇ（ｘ）＝０であり、関数Ａが線形であり、かつシステムが過剰決定である特殊な場合、式（２）の解法は、疑似逆関数Ａ†を用いる。疑似逆関数は、システムが劣決定であり、かつｇ（ｘ）＝‖ｘ‖２である、すなわち、最小エネルギー解が所望されるときにも解となる。
【０００６】
ＳＡＲ画像形成
ＳＡＲ画像の形成は、逆問題として定式化することができる。特に、理想的なＳＡＲ取得プロセスは、線形系、すなわち、式（１）の実例とみなすことができる。ここで、ｘは、表面反射率の２次元（２−Ｄ）画像を表し、ｙは、受信信号を表し、関数Ａは、線形である。
【０００７】
これ以降、Ａを用いて線形取得関数を表す。受信信号ｙは、通常、ＳＡＲ画像に類似した２−Ｄ形式で構成される。各反射パルスのサンプルは、行ベクトルを形成するとみなすことができ、異なる反射パルスが一緒にスタックされて、データの行列を形成する。関数Ａが線形であることを強調するために、表記Ａｘを用いてグラウンド信号（画像）ｘへの関数の適用Ａ（・）を表す。ｘを２Ｄ画像の代わりに１つの単一列ベクトルに再構成すると、Ａｘは、行列乗算とみなすことができる。
【０００８】
データ行列ｙの各行は、取得プラットフォームの経路に沿った該取得プラットフォームの異なる位置に対応する。経路に沿った次元は、アジマスと呼ばれる。各列は、パルスの送信からの遅延に対応する。遅延は、パルスが進んだ距離に比例するので、この第２の次元は、レンジと呼ばれる。
【０００９】
逆問題の観点から、いくつかの既知の画像形成手順は、疑似逆関数Ａ†の近似を求めることとして解釈することができる。近似によって、ＳＡＲ画像ｘの推定を効率的に実施することが可能になる。例えば、多くの手順が、パルス圧縮に依存する。パルス圧縮は、送信パルスとの相関を通じた受信信号の近似逆畳み込みを指す。相関は、高速フーリエ変換（ＦＦＴ）を用いることによって非常に効率的にすることができる。パルス圧縮は、受信パルスがナイキストレートでサンプリングされている場合、少なくとも対象の帯域幅内で厳密な逆関数に対応し、受信パルスが僅かにオーバーサンプリングされている場合（これは、実際のＳＡＲシステムにおいて通例である）、疑似逆関数に密に対応する。
【００１０】
最も広く用いられている画像形成手順のうちの１つは、チャープスケーリングアルゴリズム（ＣＳＡ）である（参照により本明細書に援用される、Runge他「A novel high precision SAR focusing algorithm based on chirp scaling」（International Geoscience and Remote Sensing Symposium（IGARSS），May 1992，pp. 372-375））。
【００１１】
図１に示すように、入力信号１０１にアジマスＦＦＴ１１０が適用され、続いて、チャープスケーリング１２０、レンジＦＦＴ１３０、レンジ圧縮１４０、レンジＩＦＦＴ１５０、チャープスケーリング１６０、アジマス圧縮１７０、およびアジマス逆ＦＦＴ（ＩＦＦＴ）１８０が適用され、出力信号１０９が生成される。
【００１２】
チャープスケーリングは、線形チャープ信号の時間変動遅延を、２つの他のチャープ信号と乗算することによって概算する。ＣＳＡでは、取得プラットフォームの動きによって引き起こされたレンジマイグレーション、所与のターゲットへの距離の変動、したがって関連する遅延を補正するのにチャープスケーリングが用いられる。ＣＳＡは、その組成に起因して、ＦＦＴおよび乗算に関してのみ、非常に効率的である。この理由により、ＣＳＡは、後述するように、取得関数およびその随伴関数を求める本方法の基礎を成す。
【００１３】
飽和および逆問題
データ振幅における大きな変動および精度の低い量子化器の使用に起因して、データ値の飽和は、ＳＡＲ信号取得において大きな問題である。ＣＳＡ等の従来の画像形成方法は、通常、データ飽和に明確に対処しない。
【００１４】
信号飽和は、アナログ／デジタル（Ａ／Ｄ）変換システムにおいて遭遇する、非常に一般的な問題である。なぜなら、電子部品は、有限の電圧レンジを有し、信号振幅は、それを超えて変動してはならないためである。加えて、量子化器は、有限数の量子化レベルを用いる。取得した信号振幅が飽和閾値±Ｔを超えて変動する場合、該取得した振幅は、±Ｔに飽和する。
【００１５】
飽和自体は、望ましくないが、厳しい量子化の存在下では、飽和は、有益となり得る。特に、信号の利得を増大させて飽和を促すことによって、非飽和測定値における信号対量子化雑音比も増大する。飽和を適切に考慮に入れた場合、全体再構成性能を増大させることができる。
【発明の概要】
【発明が解決しようとする課題】
【００１６】
飽和に対処する最も一般的な手法、すなわち、飽和が発生しなかったかのように再構成において飽和値をそのまま用いることによって、重大なアーティファクトが生じる。再構成誤差は、一貫した再構成手法、すなわち、再構成された信号推定値が再取得された場合に、同じ飽和を生じることを確保することを用いて低減することができる。一貫した再構成は、量子化され飽和した測定値とともに用いられた場合に、大幅に再構成誤差を改善する。
【課題を解決するための手段】
【００１７】
本発明の実施の形態によれば、合成開口レーダー（ＳＡＲ）画像の形成が逆問題として定式化される。本発明の実施の形態は、データの飽和が存在するとき、特に一方向二次コスト関数を最適化して受信データとの一貫性をもたらすことによって、ＳＡＲ画像の品質を増大させる。ＳＡＲ取得プロセスは、線形関数を用いてモデリングされる。関数は、実施するのに効率的であり、反復最適化プロセスにおいて、該関数の随伴関数とともに用いることができる。本方法は、アジマスアンビギュイティを制御し、画像モデルを組み込んで飽和したＳＡＲデータからの再構成を可能にする。
【発明の効果】
【００１８】
本発明の実施の形態は、合成開口レーダー（ＳＡＲ）画像の再構成を逆問題として定式化する。ＳＡＲは、リモートセンシングにおいて、比較的小さなアンテナを用いて高解像度で広域のレーダー画像を得るのに広く用いられている。高解像度は、通常、平面アンテナまたは衛星アンテナであるアンテナを搭載するプラットフォームの動きを利用して、はるかに大きな開口を合成することによって可能になる。ＳＡＲにおける取得プロセスは、プラットフォームが動くときに定期的な間隔でパルスを送信し、反射を記録して下方の表面の反射率に関する情報を得ることを伴う。
【図面の簡単な説明】
【００１９】
【図１】従来のチャープスケーリングアルゴリズムのブロック図である。
【図２】本発明の実施の形態による合成開口レーダー画像形成に用いられる取得関数のブロック図である。
【図３】本発明の実施の形態によるＳＡＲ画像を再構成する方法の流れ図である。
【発明を実施するための形態】
【００２０】
有限の動的レンジを有する線形系を
【数４】
によって表す。ここで、Ａは、システム取得関数であり、Ｓ（・）は、要素単位で入力信号の全ての成分に適用された非線形スカラー飽和関数
【数５】
である。
【００２１】
従来の最小二乗忠実度基準の代わりに、飽和測定値に対し一方的な二次ペナルティを強制し、一貫した画像再構成を確保する。より具体的には、二次ペナルティについて、コスト関数
【数６】
を用いる。ここで、（ｙ）−＝ｍｉｎ（ｙ，０）であり、（ｙ）＋＝ｍａｘ（ｙ，０）であり、ｉは、データの実数部および虚数部を別個にインデックス付けする。利点として、コスト関数は、一貫性のない解にペナルティを課し、最適化が単純であり、その二次特性に起因して飽和が発生する前に加法性雑音を低減する。
【００２２】
式（５）において、飽和にロバストなコスト関数を最小化することによって、多くの場合に受入可能な解がもたらされるが、信号モデルも用いて、データの飽和および量子化に固有のアンビギュイティをさらに解消することができる。
【００２３】
ＳＡＲ画像形成
ＡおよびＡＨの実施
式（５）のコスト関数は、種々の反復最適化手順を用いて最小化することができる。ほとんどの手順は、取得関数Ａおよびその随伴関数ＡＨを繰返し計算することを必要とする。これらの関数の効率的な実施は、画像の高次元性に起因して、ＳＡＲ画像再構成において、特に重要である。
【００２４】
効率性を達成するために、本実施の形態は、上述したチャープスケーリングアルゴリズム（ＣＳＡ）に大まかに基づく。ＣＳＡは、疑似逆関数Ａ†の近似であることから、未加工のＳＡＲデータ２０９を再生するのに用いることができる取得関数Ａは、図２に示すブロック図を有する。これは、本質的に従来のＣＳＡの関数の反対である。
【００２５】
入力信号２０１にアジマス高速フーリエ変換（ＦＦＴ）２１０が適用され、続いてアジマス周波数応答Ｈａ（Ｒ０，ｆη）の乗算２２０、第１のチャープスケーリング２３０、レンジＦＦＴ２４０、レンジ周波数応答Ｐ２ｄｆ（ｆτ，ｆη）の乗算２５０、レンジ逆ＦＦＴ（ＩＦＦＴ）２６０、第２のチャープスケーリング２７０、およびアジマスＩＦＦＴ２８０が適用され、未加工ＳＡＲデータである出力信号２０９が生成される。
【００２６】
上記の関数は、当該技術分野において既知のメモリおよび入力／出力インターフェースに接続されたプロセッサ２００、またはディスクリートハードウェア構成要素において実施することができることを理解されたい。
【００２７】
ＳＡＲ画像を生成する随伴関数を実施するには、関数を合成したものである随伴関数が、個々の随伴関数を逆の順序で合成したものであるということを利用する。ＦＦＴの随伴関数は、スケール係数次第でＩＦＦＴであり、関数の乗算の随伴関数は、同じ関数の複素共役の乗算である。このため、随伴関数は、アジマス圧縮関数がアジマス周波数応答
【数７】
と置き換えられた、ＣＳＡと同じブロック図を有する。ここで、＊は、複素共役を示し、ａは、アジマス応答を意味し、Ｒ０は、最接近のレンジを示し、ｆηは、アジマス周波数である。
【００２８】
このとき、レンジ圧縮関数は、フーリエ変換
【数８】
であり、２ｄｆは、２次元周波数領域を意味し、ｆτは、レンジ周波数である。
【００２９】
次に、本発明者らの取得関数Ａの実施を、より詳細に説明する。画像が点目標物の長方形格子を含むと仮定する。単位振幅の単一点目標物による応答に焦点を当てる。関数Ａは、アジマスにおいてシフト不変であるので、レンジドップラー（レンジ時間τ、アジマス周波数ｆη）領域における応答を、等しく特定することができる。最接近のレンジＲ０における単位点目標物の場合、レンジドップラー応答は、
【数９】
によって近似的に与えられる。ここで、ｒｄは、レンジドップラー領域を意味し、Ｈａ（Ｒ０，ｆη）は、アジマス周波数応答であり、ｐ（τ，ｆη）は、送信パルスであり、Ｄ（ｆη）は、レンジマイグレーション係数であり、ｃは、光速である。パルスｐ（τ，ｆη）は、線形チャープであると仮定され、レンジドップラー領域において、ｆηの関数である。
【００３０】
式（７）における応答は、レンジマイグレーション係数Ｄ（ｆη）の存在に起因して、レンジがシフト不変でない。しかしながら、線形チャープパルスの場合、変更されたチャープを用いた畳み込み、その前後におけるチャープスケーリング関数との乗算によって、チャープスケーリングを用いて、式（７）におけるレンジ応答を近似することができる。図２に示すように、畳み込みは、ＦＦＴを用いてフーリエ変換Ｐ２ｄｆ（ｆτ，ｆη）との乗算として実施される。
【００３１】
図２は、取得関数Ａを計算する複雑度が図１のＣＳＡを計算する複雑度と同じであることを示している。より高速な計算のために、信号と乗算された図２の関数の全てを予め求め、格納することができる。離散時間型の実施では、Ｈａ（Ｒ０，ｆη）およびＰ２ｄｆ（ｆτ，ｆη）を最小のエイリアシングで表すのに十分高いサンプリングレートを用いるように注意すべきである。
【００３２】
基本最適化
図３に示すように、一次勾配降下プロセスを用いて、式（５）のコスト関数を最小にする。より洗練された二次手順も検討することができる。この最適化は、初期構成ｘを推定すること（３１０）から始まる。初期再構成ｘは、例えば、ＣＳＡを用いて計算することもできるし、０に設定することもできるし、任意に設定することもできる。
【００３３】
各反復３６０中に、
【数１０】
を用いて信号ｙが再生され（３２０）、再生信号は、入力信号と比較され（３３０）、式（６）において残差ベクトルｒが求められ、式（５）において二次ペナルティを用いて総コストが求められる。推定値を調整する（３４０）ために、
【数１１】
を用いて、式（５）に対応する勾配ｇが計算される。ここで、ＡおよびＡＨは、上述したように求められる。次に、サイズτのステップが勾配ｇと反対方向に取られる。すなわち、反復は、終了条件３５０に達するまで
【数１２】
として更新される。
【００３４】
厳密な線探索を通じて、負の勾配−ｇに対して平行な線に沿って、関数ｆ（ｙ，Ａｘ）を最小化するステップサイズτが最適に選択される。この線探索は、１−Ｄ凸区分的二次関数を効率的に最適化する。
【００３５】
いくつかの可能な終了条件が存在する。例えば、再構成された画像推定値ｘの変化が小さい場合、または比較ステップにおける式（５）の二次ペナルティが小さい場合、反復は、終了する場合がある。本方法は、終了条件に達したときに、再構成画像を出力する（３７０）。上記のステップは、当該技術分野において、既知のメモリおよび入力／出力インターフェースに接続されたプロセッサ３００において実行することができる。
【００３６】
アジマスアンビギュイティおよびモデル不整合
アンダーサンプリングに起因したアジマス周波数領域におけるアンビギュイティは、多くのＳＡＲ画像形成手順に影響を与える問題である。パルス繰返し周波数に等しいアジマスサンプリングレートは、通常、アジマス応答Ｈａ（Ｒ０，ｆη）の帯域幅に適応するには不十分であり、画像化領域の幅を減らさずに、増大させることができない。アンダーサンプリングの結果、エイリアシングが生じ、エイリアシングは、ゴースト像が特定のアジマスおよびレンジオフセットにおいて重なる形態で現れる。エイリアシングからの復元は、画像に関する追加情報なしでは困難である。例えば、画像自体が帯域制限されていると仮定することはできない。
【００３７】
上述した基本方法では、真の取得システムと、再構成において用いられるモデルとの間の不整合によって、さらなる難題が提起される。例えば、アジマス応答Ｈａ（Ｒ０，ｆη）は、アンテナの長さおよび重み付けパターン等のパラメータに依拠する。これらのパラメータの値を精密に知ることができず、不正確な結果につながる可能性がある。特に、本発明者らの経験によれば、モデル不整合は、アジマスアンビギュイティによって生じたゴースト像を大きくさせる傾向にあることが示唆される。ゴースト発生を低減する１つの方法は、データを非適応的にフィルタリングするか、またはエイリアシングによって最も影響を受けるアジマス周波数を、空間において適応的に除去することである。しかしながら、その手法は、再構成に利用可能なデータ量、したがって、結果の画像の解像度を低減する可能性がある。
【００３８】
代替的な手法は、上記の基本方法における各反復の終了時に、ゴースト除去を実行することである。ＳＡＲパラメータの知識に基づいて、対応する真の目標物に対するゴーストのロケーションを予測する能力を利用する。本方法の第１の反復において、まず、低解像度でアンビギュイティのない画像を形成し、次に、ゴーストの距離特性で隔てられた領域間の局所電力を比較することによって、主要なゴースト源である可能性が高い画像領域が特定される。後続の反復において、現在の画像を１度に１つのアジマス線ずつ該現在の画像自体と相関させ、次に、推定値を画像から減算することによって、主要なゴースト源によるゴーストが推定される。推定および減算を主要なゴースト源に限定することによって、ゴーストの効率的な除去が可能になる。
【００３９】
逆問題との関連で、本方法は、特に、ある特定のロケーションおよびシフトにおける画像の自己相関において、式（２）の再構成最適化問題に対し、明示的な制約を付加することとみなすことができる。これは、大きな自己相関値は、ゴーストのない画像内に生じるはずがないという前提に基づく。代替の形態は、そのような知識を、信号モデル項ｇ（ｘ）に組み込む。
【００４０】
信号モデル
再構成正確度における大きな改善は、再構成に信号モデルを組み込むことによって達成することができる。上述したゴースト低減モデルに加えて、自然な画像は、大きな構造を呈する。本発明者らが課すモデルは、ＳＡＲ画像の大きさのウェーブレット変換が、スパースであるというものである。本明細書で定義されるとき、スパース信号は、複数の非ゼロ値の数が、ゼロ値の数よりも大幅に小さい場合に１である。
【００４１】
ｌ１最小化の固定点継続手順に類似した方式のソフト閾値処理を用いて、モデルを課す。参照により本明細書に援用されるHale他「Fixed-point continuation for l1-minimization: methodology and convergence」（SIAM J. Optim., vol. 19, no. 3, pp. 1107-1130, Oct. 2008）を参照されたい。
【００４２】
本方法の全ての反復の終了時に、画像またはその大きさを変換し、ソフト閾値処理を実行し、逆変換する。画像の大きさのみが変換された場合、同じ位相を維持する。このモデルの動機付けは、ＳＡＲ画像の大きさが、自然な画像と類似して見え、したがって、ウェーブレット変換を用いて圧縮可能であるということである。一方、画像化されている特徴の高さの変動に依拠し、深刻なラッピングを受ける画像の位相は、再構成において利用することができる構造を提供しないようである。
【特許請求の範囲】
【請求項１】
合成開口レーダー（ＳＡＲ）システムによって取得された入力信号を処理する方法であって、前記入力信号は、飽和しており、該方法は、
前記入力信号を生成した初期再構成を推定するステップと、
再生信号を生成するために、前記再構成から前記入力信号を再生するステップと、
前記再生信号を前記入力信号と比較するステップと、
前記推定された再構成を前記比較に基づいて調整するステップと、
再構成されたＳＡＲ画像を形成するために、前記再生するステップから終了条件に達するまで反復するステップと
を含み、これらのステップは、プロセッサにおいて実行される、合成開口レーダー（ＳＡＲ）システムによって取得された入力信号を処理する方法。
【請求項２】
前記推定するステップは、
前記入力信号にアジマス高速フーリエ変換（ＦＦＴ）を適用するステップと、
前記アジマスＦＦＴの出力に第１のチャープスケーリングを適用するステップと、
前記第１のチャープスケーリングの出力にレンジＦＦＴを適用するステップと、
前記レンジＦＦＴの出力にレンジ周波数応答を乗算するステップと、
前記レンジ周波数応答の出力にレンジ逆ＦＦＴ（ＩＦＦＴ）を適用するステップと、
前記レンジＩＦＦＴの出力に第２のチャープスケーリングを適用するステップと、
前記第２のチャープスケーリングの出力にアジマス周波数応答を乗算するステップと、
前記再構成されたＳＡＲ画像を形成するために、前記アジマス周波数応答の出力にアジマスＩＦＦＴを適用するステップと
をさらに含む請求項１に記載の方法。
【請求項３】
前記再生するステップは、
前記入力信号にアジマス高速フーリエ変換（ＦＦＴ）を適用するステップと、
前記アジマスＦＦＴの出力にアジマス周波数応答を乗算するステップと、
前記アジマス周波数応答の出力に第１のチャープスケーリングを適用するステップと、
前記第１のチャープスケーリングの出力にレンジＦＦＴを適用するステップと、
前記レンジＦＦＴの出力にレンジ周波数応答を乗算するステップと、
前記レンジ周波数応答の出力にレンジ逆ＦＦＴ（ＩＦＦＴ）を適用するステップと、
前記レンジＩＦＦＴの出力に第２のチャープスケーリングを適用するステップと、
未加工ＳＡＲデータを形成するために、前記第２のチャープスケーリングの出力にアジマスＩＦＦＴを適用するステップと
をさらに含む請求項１に記載の方法。
【請求項４】
前記調整するステップは、信号モデルに従う請求項１に記載の方法。
【請求項５】
前記反復が終了した後に信号モデルが組み込まれる請求項１に記載の方法。
【請求項６】
前記信号モデルは、前記再構成されたＳＡＲ画像のウェーブレット変換または該再構成されたＳＡＲ画像の大きさのウェーブレット変換がスパースであることを課す請求項４または５に記載の方法。
【請求項７】
変換係数を生成するために、前記再構成されたＳＡＲ画像または前記大きさを変換するステップと、
閾値未満の値を有する前記係数をゼロに設定することによって前記変換係数を閾値処理するステップと、
前記再構成されたＳＡＲ画像を形成するために、前記閾値処理された係数を逆変換するステップと、
をさらに含む請求項６に記載の方法。
【請求項８】
ゴースト像を低減するステップをさらに含む請求項１に記載の方法。
【請求項９】
前記低減するステップは、
前記ゴースト像のロケーションを特定するステップと、
前記ロケーションから前記ゴースト像を除去するステップと
をさらに含む請求項８に記載の方法。
【請求項１０】
前記低減するステップは、前記調整するステップ中に適用される請求項８に記載の方法。
【請求項１１】
前記低減するステップは、前記反復するステップの後に適用される請求項８に記載の方法。
【請求項１２】
前記比較するステップは、飽和一貫性を強制する請求項１に記載の方法。
【請求項１３】
前記調整するステップは、勾配降下を用いて実行される請求項１に記載の方法。
【請求項１４】
前記勾配降下は、線探索を用いる請求項１３に記載の方法。
【請求項１】
合成開口レーダー（ＳＡＲ）システムによって取得された入力信号を処理する方法であって、前記入力信号は、飽和しており、該方法は、
前記入力信号を生成した初期再構成を推定するステップと、
再生信号を生成するために、前記再構成から前記入力信号を再生するステップと、
前記再生信号を前記入力信号と比較するステップと、
前記推定された再構成を前記比較に基づいて調整するステップと、
再構成されたＳＡＲ画像を形成するために、前記再生するステップから終了条件に達するまで反復するステップと
を含み、これらのステップは、プロセッサにおいて実行される、合成開口レーダー（ＳＡＲ）システムによって取得された入力信号を処理する方法。
【請求項２】
前記推定するステップは、
前記入力信号にアジマス高速フーリエ変換（ＦＦＴ）を適用するステップと、
前記アジマスＦＦＴの出力に第１のチャープスケーリングを適用するステップと、
前記第１のチャープスケーリングの出力にレンジＦＦＴを適用するステップと、
前記レンジＦＦＴの出力にレンジ周波数応答を乗算するステップと、
前記レンジ周波数応答の出力にレンジ逆ＦＦＴ（ＩＦＦＴ）を適用するステップと、
前記レンジＩＦＦＴの出力に第２のチャープスケーリングを適用するステップと、
前記第２のチャープスケーリングの出力にアジマス周波数応答を乗算するステップと、
前記再構成されたＳＡＲ画像を形成するために、前記アジマス周波数応答の出力にアジマスＩＦＦＴを適用するステップと
をさらに含む請求項１に記載の方法。
【請求項３】
前記再生するステップは、
前記入力信号にアジマス高速フーリエ変換（ＦＦＴ）を適用するステップと、
前記アジマスＦＦＴの出力にアジマス周波数応答を乗算するステップと、
前記アジマス周波数応答の出力に第１のチャープスケーリングを適用するステップと、
前記第１のチャープスケーリングの出力にレンジＦＦＴを適用するステップと、
前記レンジＦＦＴの出力にレンジ周波数応答を乗算するステップと、
前記レンジ周波数応答の出力にレンジ逆ＦＦＴ（ＩＦＦＴ）を適用するステップと、
前記レンジＩＦＦＴの出力に第２のチャープスケーリングを適用するステップと、
未加工ＳＡＲデータを形成するために、前記第２のチャープスケーリングの出力にアジマスＩＦＦＴを適用するステップと
をさらに含む請求項１に記載の方法。
【請求項４】
前記調整するステップは、信号モデルに従う請求項１に記載の方法。
【請求項５】
前記反復が終了した後に信号モデルが組み込まれる請求項１に記載の方法。
【請求項６】
前記信号モデルは、前記再構成されたＳＡＲ画像のウェーブレット変換または該再構成されたＳＡＲ画像の大きさのウェーブレット変換がスパースであることを課す請求項４または５に記載の方法。
【請求項７】
変換係数を生成するために、前記再構成されたＳＡＲ画像または前記大きさを変換するステップと、
閾値未満の値を有する前記係数をゼロに設定することによって前記変換係数を閾値処理するステップと、
前記再構成されたＳＡＲ画像を形成するために、前記閾値処理された係数を逆変換するステップと、
をさらに含む請求項６に記載の方法。
【請求項８】
ゴースト像を低減するステップをさらに含む請求項１に記載の方法。
【請求項９】
前記低減するステップは、
前記ゴースト像のロケーションを特定するステップと、
前記ロケーションから前記ゴースト像を除去するステップと
をさらに含む請求項８に記載の方法。
【請求項１０】
前記低減するステップは、前記調整するステップ中に適用される請求項８に記載の方法。
【請求項１１】
前記低減するステップは、前記反復するステップの後に適用される請求項８に記載の方法。
【請求項１２】
前記比較するステップは、飽和一貫性を強制する請求項１に記載の方法。
【請求項１３】
前記調整するステップは、勾配降下を用いて実行される請求項１に記載の方法。
【請求項１４】
前記勾配降下は、線探索を用いる請求項１３に記載の方法。
【図１】
【図２】
【図３】
【図２】
【図３】
【公開番号】特開２０１２−１６８１６８（Ｐ２０１２−１６８１６８Ａ）
【公開日】平成２４年９月６日（２０１２．９．６）
【国際特許分類】
【外国語出願】
【出願番号】特願２０１２−８１２５（Ｐ２０１２−８１２５）
【出願日】平成２４年１月１８日（２０１２．１．１８）
【出願人】（５９７０６７５７４）ミツビシ・エレクトリック・リサーチ・ラボラトリーズ・インコーポレイテッド (484)
【住所又は居所原語表記】２０１ ＢＲＯＡＤＷＡＹ， ＣＡＭＢＲＩＤＧＥ， ＭＡＳＳＡＣＨＵＳＥＴＴＳ ０２１３９， Ｕ．Ｓ．Ａ．
【Ｆターム（参考）】
【公開日】平成２４年９月６日（２０１２．９．６）
【国際特許分類】
【出願番号】特願２０１２−８１２５（Ｐ２０１２−８１２５）
【出願日】平成２４年１月１８日（２０１２．１．１８）
【出願人】（５９７０６７５７４）ミツビシ・エレクトリック・リサーチ・ラボラトリーズ・インコーポレイテッド (484)
【住所又は居所原語表記】２０１ ＢＲＯＡＤＷＡＹ， ＣＡＭＢＲＩＤＧＥ， ＭＡＳＳＡＣＨＵＳＥＴＴＳ ０２１３９， Ｕ．Ｓ．Ａ．
【Ｆターム（参考）】
[ Back to top ]