JP2023526713A - Ｉｎｓａｒ画像を処理して地盤変動信号を抽出するための方法 - Google Patents

Abstract

本発明は、同一エリアのノイズが多い画像時系列を処理するための方法に関し、この方法は、入力画像時系列の各画像の各画素（ＰＸ［ｉ，ｊ］）を入力画像時系列の画像および隣接画像における選択された近隣画素と第１の線形結合によって組み合わせることにより、入力画像時系列（ＲＩ［ｔ］）から画像時系列のセットを生成することと、各々が、セットの各時系列の各画像の各画素をセットの各時系列内の画像および隣接画像における選択された近隣画素と第２の線形結合によって組み合わせるフィルタリング演算（ＦＰｓ、ｓ＝１～３、５、７～１１）をカスケード方式でセットに適用することと、セットの各時系列を単一の画像に低減するために画像組合せ演算（ＳＬｓ、ｓ＝４、６）を行うことと、エリアのモデル画像（ＥＭ）をフィルタ画像としてセット内に導入することと、セット内の各画像を第３の線形結合によって出力画像（ＰＩ）に組み合わせることとを備える。
【選択図】図１

Description

本発明は、たとえば地球観測衛星からの合成開口レーダ（ＳＡＲ）画像などの画像、特に干渉ＳＡＲ（ＩｎＳＡＲ）画像を処理するための方法およびデバイスに関する。本発明は、そのような画像から、構造過程および火山活動を含む自然現象に起因する地盤変動、または帯水層の枯渇、（排水の注入および地熱還元を含む）注水、および（石油、ガス、鉱石を含む）採収活動を含む人間の活動に起因する地盤変動に特に関連する有用な情報を抽出するために適用される。具体的には、本発明は、断層に沿ったスロースリップ事象の検出、活断層に沿った先行滑りの検出に適用される。

構造プレート間の相対運動は、原則としてＩｎＳＡＲによる測定が可能な活断層の滑りによって調整される。大半の断層は、ほとんどの時間動いていないというのが従来のパラダイムである。大きな地震の際、蓄積している弾性応力が解放され得るが、一部の断層または断層の一部における滑りは主として非地震性であり得るため、地震災害を引き起こすことはほとんどない。予想外の滑り事象の発生を明らかにしたこの２０年間の研究により、この考えは重要な点で安易であることが示されている。実際、昨今では、低周波地震、スロースリップ事象およびそれに伴う構造微震、ならびに連続的な非地震性滑りを含む様々な種類の滑り現象が識別されている。

地震は、人口密度の高い地域にとって大きな脅威であるが、個々の地震の発生に関する短期または長期予測を提供することは未だ不可能である。地震の前に、構造断層がゆっくりと滑り始めて急速な断裂の前に加速するという核生成段階が先行することが理論的に示されている。そのような核生成段階の検出は、大きな地震の短期予測を可能にする。しかし、一部の断層は、大きな地震を発生させることなく、また十分に密度の高い観測セットを伴わず、偶発的にゆっくりと滑るので、破壊的な事象の前兆である核生成段階と、無害なスロースリップ事象とを効率的に区別することは難しい。また、非地震性滑りと地震性滑りとは密接に関連しており、断層クリープは地震の前およびしばしば後に生じると考えられるので、一般に、断層滑りの機械的および物理的記述が大きく進歩することの必要性が示される。

１９９０年代初めに開発された合成開口レーダ干渉法（ＩｎＳＡＲ）は現在、水循環過程、火山過程、および構造過程に起因する地盤変動を測定するために日常的に用いられている。現在は、軌道衛星が地球表面の日々の画像を提供し、レーダシステムを有する干渉法の概念を用いて、地上と衛星システムとの間の視距離の展開のマップが導出され得る。これらのマップは、そのような距離の連続的な監視を提供するために組み立てられ得る。この視距離は、経時的な地盤運動に直接関連するので、今日、大きな地震の前兆であろうとなかろうと、地球上のあらゆる場所でのスロースリップ事象を検出することが可能である。Ｓｅｎｔｉｎｅｌ １ ＩｎＳＡＲ衛星コンステレーションは、全世界の全ての変動領域の系統的マッピングが可能な能力を有しており、断層変動から大気信号を分離するためのロバストかつ自動的な方法が開発可能であれば、構造断層の全ての滑りモードの背景にある機構を明らかにする。Ｓｅｎｔｉｎｅｌ １コンステレーションの前、ＳＡＲ画像の典型的な範囲は、最高水準のＣバンドおよびＬバンド衛星（ＥＲＳ、Ｅｎｖｉｓａｔ、ＡＬＯＳ・・・）に関して、（レーダおよびアンテナ構成の波長に依存して）約数十メートルの画素サイズで１００×１００ｋｍであった。Ｓｅｎｔｉｎｅｌ １コンステレーションの場合、画像のサイズは平均３倍の大きさであり、取得の頻度は最低でも５倍である。Ｓｅｎｔｉｎｅｌ １衛星コンステレーションの発射は、６日間の再現期間で世界中の全ての活動性変動領域の系統的レーダマッピングを提供するという点で変革的である。今後の衛星発射（ＮＩＳＡＲなど）により、解像度および周波数が更に高いＩｎＳＡＲデータが提供され、自動検出ツールの必要性が更に高まる。

ＩｎＳＡＲ技術は、地震に起因する大きな変位、氷床運動、帯水層の枯渇に関連する小さな変位、地震間の変動、ゆっくりと動く地滑り、およびスロースリップ事象を監視するための適用に成功している。たとえば地震時の変位場や火山構造性事象などの急激な大振幅変動信号もまたＩｎＳＡＲによって日常的に分析される。同様に、比較的長期間にわたる、ゆっくりではあるが着実な変動の蓄積もまた、結果的に生じる変動が十分に大きいことにより、ＩｎＳＡＲを用いて容易に検出される。ただし、これらの測定は一般に、専門家による変動データセットの入念な手動調査がなければ成功しないものであり、そのため大規模研究が阻まれる。

２つのＳＡＲ画像を組み合わせてインターフェログラムにすること（ＩｎＳＡＲ）により、広範な面積にわたる地盤運動の測定が可能である。インターフェログラムのフェーズは、レーダキャリアと地面との間のレーダ波の双方向移動時間に影響を及ぼす任意の物理現象の関数として変化する。これらの影響は、取得間の軌道の変化ならびに大気および地盤変動の変化を含む。軌道は周知であるため、変動の測定値を不明瞭にする最も重要なノイズ源は、温度、圧力、および大気中の水蒸気の空間および時間的変化である。対流圏遅延および電離層遅延によりＩｎＳＡＲ処理に大きなエラーが生じ、しばしば１０～２０ｃｍのエラーとなり、これは関心信号よりも著しく大きい場合がある。レーダ信号における遅延は、主に、対流圏に存在する水蒸気に起因する。これらの遅延は、１）対流圏における乱気流、および２）局所地形に密接に関わる構成要素という主に２つの原因が識別されている。この問題を軽減し、ノイズを低減するために、本文書において２つの主なアプローチが展開される。第１のアプローチは、大気ノイズが確率論的な性質を持ち、平均に近付き得ることを前提として、時間的に平滑化またはフィルタリングすることから成る。これらのアプローチは、フィルタリングおよび相関法に基づいてインターフェログラム時系列分析を積み重ねることを含む。しかし、大気ノイズは実際には時間および空間の両方において非確率論的であり、平滑化は、地面からの信号の解像度の低下を招く。言い換えると、これらの方法は、関心信号を部分的にマスクアウトすることが示されている。

第２のアプローチは、気象データ、分光放射計および分光計のデータから得た水蒸気マップ、（グローバルナビゲーション衛星システムからの）地理的位置データ、および天候モデルおよびシミュレーションによる推定値を含む追加のデータセットに依拠して、遅延を推定し、地盤変動推定値におけるノイズを軽減する。これらの様々なデータセットは全て、それぞれの利点および欠点を有する。大気モデルは、系統的または局所的エラーに苛まれる。水蒸気マップは、高い空間分解能を有するが、時間分解能が比較的乏しく、雲が存在する場合には制限される。対照的に、地理的位置データは、最も高い時間分解能の大気遅延推定値を提供するが、地表の局を必要とし、データは、利用可能な局に空間的に限定される（局間の補間は明瞭ではなく、その品質を評価することは難しい）。現時点で、ＩｎＳＡＲ地盤運動推定値への大気に起因するエラーの影響を軽減するために普遍的に用いられる信頼性の高いアプローチは存在しない。

したがって、地震事象間の歪み蓄積や地震事象後の運動による変動場などの低振幅で長波長の変動場の検出は、干渉計の非相関性、不正確な衛星軌道、および大気が衛星信号を低速化させると地面と衛星との間の距離が増加して見え、地面が動いたようになるという大気伝搬遅延によって、未だに課題を残している。この目標を達成するためには、ゆっくりとした変動のＩｎＳＡＲ測定に大きな進歩がなくてはならない。

したがって、ＩｎＳＡＲ画像においてデータが利用可能であるが、構造断層に沿ったこれらの滑り事象の系統的検出は、未だ可能ではない。実際、既存の解決策は、信号が大気ノイズと構造信号との混合である地盤変動のＩｎＳＡＲ画像時系列の目視検査による、人手に基づく検出に依拠している。

よって、考えられる最も小さい様々な変動信号を強調するために、たとえばＩｎＳＡＲ画像時系列などの大規模な画像セットの時系列を自動的に処理し、大気による摂動を除去または減衰することが必要である。また、無害なスロースリップ事象と、地震の前兆に関連する小さな事象とを確信的に区別するために、様々なスケールでの低振幅および／または長波長の地盤変動を系統的かつ全体的に検出および特徴化することも必要である。また、専門家の解釈を必要とせずに断層変動を識別することも望まれ得る。

ノイズに作用される同一エリアの画像時系列を処理するために方法が説明される。この方法は、入力画像時系列を受信する第１のフィルタリングステップと、出力画像を提供する最後のフィルタリングステップと、中間フィルタリングステップとを備える複数の連続的なフィルタリングステップを実行することであって、第１のフィルタリングステップおよび中間フィルタリングステップの各々は、最初に入力画像時系列を含むフィルタ画像時系列セットを変換し、このセットは、セットの各時系列の各画像の各画素をセットの画像時系列の画像および隣接画像における選択された近隣画素と線形結合によって組み合わせることによって、複数のフィルタ画像時系列を生成することによって変換され、各線形結合は、それぞれの加重係数セットを用いて、生成されたフィルタ画像のうち１つの画素となることと、１つ以上の組合せ演算を実行することであって、各組合せ演算は、セットのフィルタ画像時系列の各々における隣接画像サブセットの画像を組み合わせることによってセットのフィルタ画像時系列の各々における画像の数を低減するために、２つの連続する中間フィルタリングステップの間に行われ、組合せ演算のうち最後の１つは、セットの各フィルタ画像時系列を単一のフィルタ画像に低減することと、を含んでよい。

実施形態によると、方法は、最後の組合せ演算の後、セット内のエリアのモデル画像を追加のフィルタ画像として導入することを更に備え、モデル画像の導入の後に、フィルタリングステップのうち１つ以上が後続する。

実施形態によると、線形結合の各々は、線形結合の結果に加算されるバイアス係数を含む。

実施形態によると、加重係数の各々は、乗算される画素値の正または負の符号に依存する値を有し、または、線形結合の各々の結果は、正規化線形ユニット関数によって変換される。

実施形態によると、方法は、ランダムに選択されたパラメータを用いて地盤変動モデルから第１の画像時系列を生成することと、様々なノイズ信号のモデルを用いて、第１の画像時系列から訓練画像時系列を生成することと、加重係数の値を調整するために、生成された訓練時系列を用いることとを更に備える。

実施形態によると、加重係数の値は、訓練画像時系列およびモデルコスト関数を用いてモデルコスト関数の結果を最小にするように反復勾配に基づく降下最小法を適用することによって反復的に調整される。

実施形態によると、入力画像時系列からフィルタ画像時系列のセットを生成することは、式

を用いて行われ、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ］は、フィルタ画像時系列ｆのフィルタ画像のうち１つの１画素であり、ＰＸ［ｉ，ｊ，ｔ］は、入力画像時系列の１つの画像の１画素であり、Ｗ［ｌ，ｍ，ｕ，ｆ］は、フィルタ画像時系列ｆに関する第１の係数のうち１つであり、Ｂ［ｆ］は、時系列ｆに関するバイアス係数であり、ＬＲ（）は、正規化線形ユニット関数である。

実施形態によると、フィルタリング演算の第１のフィルタリング演算の各線形結合は、式

を適用し、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ’］は、フィルタ画像時系列ｆ’のフィルタ画像のうち１つの１画素であり、Ｗ［ｓ，ｌ，ｍ，ｕ，ｆ，ｆ’］は、フィルタリング演算ｓに関するフィルタ画像時系列ｆ’に関する第２の加重係数のうち１つであり、Ｂ［ｓ，ｆ’］は、時系列ｆ’およびフィルタリング演算ｓに関するバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である。

実施形態によると、フィルタリング演算の第２のフィルタリング演算の各線形結合は、式

を適用し、式中、ＰＸ［ｉ，ｊ，ｆ’］は、フィルタ画像ｆ’のうち１つの１画素であり、Ｗ［ｓ，ｌ，ｍ，ｆ，ｆ’］は、フィルタリング演算ｓに関するフィルタ画像ｆ’に関する第２の加重係数のうち１つであり、Ｂ［ｓ，ｆ’］は、フィルタ画像ｆ’ およびフィルタリング演算ｓに関するバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である。

実施形態によると、出力画像の各画素は、式

によって計算され、式中、ＰＸ［ｉ，ｊ］は、出力画像の１画素であり、ＰＸ［ｉ，ｊ，ｆ］は、フィルタ画像ｆの１画素であり、Ｗ［ｌ，ｍ，ｆ］は、第３の係数のうち１つであり、Ｂはバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である。

実施形態によると、漏洩正規化線形ユニット関数ＬＲは、ｘ≧０の場合、ＬＲ（ｘ）＝ｘであり、ｘ＜０の場合、ＬＲ（ｘ）＝０．５ｘとなる。

実施形態によると、画像組合せ演算の各々は、以下の式：

を適用し、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ］は、フィルタ画像時系列ｆのフィルタ画像のうち１つの１画素であり、ＭＡＸ（ＰＸ［ｉ，ｊ，ｔ＋ｕ，ｆ］）は、ｕ＝－１、０、および１を有する画素値ＰＸ［ｉ，ｊ，ｔ＋ｕ，ｆ］の中から最大値を提供する関数である。

また実施形態は、上述の方法を実行するように構成されたプロセッサと、メモリとを備えるコンピュータにも関連し得る。

実施形態によると、コンピュータは、プロセッサによって制御されるグラフィックプロセッサを更に備え、プロセッサは、上述の方法の演算のうちいくつかを実行するようにグラフィックプロセッサを構成するように構成される。

また実施形態は、コンピュータによって実行されるときに、上述の方法を実行するようにコンピュータを構成するコード部分を備える、コンピュータメモリにロード可能なコンピュータプログラム製品にも関連し得る。

方法および／またはデバイスは、以下の図面および説明を参照し、より適切に理解され得る。以下の図面と共に非限定的かつ非包括的な説明が記述される。図面において、特に指定されない限り、異なる図面を通して同様の参照符号が同様の部品を示し得る。

実施形態に係る、画像を処理するための方法を実施するコンピュータの処理ステージのブロック図表である。 実施形態に係る、プロセッサの第１の処理ステージのブロック図表である。 実施形態に係る、１つの画素値を計算するために第１の処理ステージによって行われる動作を示すブロック図表である。 実施形態に係る、プロセッサの第２のステージのブロック図表である。 実施形態に係る、１つの画素値を計算するために第２の処理ステージによって行われる動作を示すブロック図表である。 実施形態に係る、プロセッサの第３の処理ステージのブロック図表である。 実施形態に係る、１つの画素値を計算するために第３の処理ステージによって行われる動作を示すブロック図表である。 実施形態に係る、プロセッサの第４の処理ステージのブロック図表である。 実施形態に係る、１つの画素値を計算するために第４の処理ステージによって行われる動作を示すブロック図表である。 実施形態に係る、プロセッサの第５の処理ステージのブロック図表である。 実施形態に係る、１つの画素値を計算するために第５の処理ステージによって行われる動作を示すブロック図表である。 画像処理方法を実施するように構成され得る従来のパーソナルコンピュータのブロック図表である。

たとえばＩｎＳＡＲ画像などの画像を処理するための方法が開示される。この方法は、
１）位相に関するコヒーレント画像の再構成、
２）衛星位置誤差を考慮に入れるための画像の補正、
３）大気ノイズの除去、および
４）変動時系列の再構成
を備える。

コンピュータ上のソフトウェアの実施形態は、以下で、コンピュータにロードされ実行されるときに、データ処理動作の実行を制御するために協働するコンピュータ可読プログラム構成要素のセットとして説明される。本発明の個々の方法ステップがコンピュータプログラムコードで実施され得ること、および、本明細書で説明する方法を実施するために様々なプログラミング言語および符号化の実施が用いられ得ることが当業者には明らかである。また、ソフトウェアに含まれるコンピュータプログラムは、本明細書で説明する特定の制御フローに限定されることは意図されておらず、コンピュータプログラムのステップのうち１つ以上は、並行または連続して行われてよい。コンピュータプログラム制御型実装の文脈で説明される動作のうち１つ以上は、代替的にハードウェア電子部品として実装されてもよい。

図１は、実施形態に従ってコンピュータによって行われる、衛星位置誤差が補正された画像から大気ノイズを除去する（ステップ３）ための方法の動作を示す。図１は、ノイズ、特に大気ノイズを除去するために画像の時系列ＲＩを連続的に処理するように構成された畳み込みニューラルネットワークＣＮＮの直列接続された処理ステージＦＰ０、ＦＰ１、ＦＰ２、ＦＰ３、ＳＬ４、ＦＰ５、ＳＬ６、ＦＰ７、ＦＰ８、ＦＰ９、ＦＰ１０、ＦＰ１１、ＦＰ１２を示す。時系列ＲＩの各画像ＲＩ［ｔ］は、画素セットＰＸ［ｉ，ｊ］を備え、ｉ＝１、・・・Ｘおよびｊ＝１、・・・Ｙである。画像時系列ＲＩは、異なる時間に撮影された、または所与の時間間隔で取得された同一シーンエリア（たとえば同じ地面領域）のＴ個の画像ＲＩ［１］、・・・ＲＩ［Ｔ］を備える。画像時系列ＲＩは、第１の処理ステージによって処理され、出力画像ＰＩを提供する最終処理または出力ステージＦＰ１２まで、ステージの各々によって計算された画像が次のステージに提供され、出力画像ＰＩでは、大気ノイズが存在せず、または強度に減衰されていながら、画像の入力時系列ＲＩに存在する地盤運動信号が存在している。

図２は、実施形態に従って、第１の処理または入力ステージＦＰ０を示す。処理ステージＦＰ０は、Ｆ個のフィルタリング成分ＦＴ［０，ｆ］を備え、ｆは０～Ｆに等しく、各フィルタリング成分ＦＴ［０，ｆ］が入力時系列ＲＩを処理する。ＦＴ［０，ｆ］の表記中の「０」は、第１の処理ステージＦＰ０に割り当てられた数を指す。各フィルタリング成分ＦＴ［０，ｆ］は、Ｔ個のフィルタ画像ＦＩ［０，ｆ，ｔ］の時系列ＦＩ［０，ｆ］を生成する。したがって、処理ステージＦＰ０は、フィルタ画像ＦＩ［０，ｆ，ｔ］のＦ個の時系列を生成する。

図３は、実施形態に従って、第１の処理ステージＦＰ０の各フィルタリング成分ＦＴ［０，ｆ］によって行われる動作を示す。各フィルタ画像ＦＩ［０，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］は、対応するフィルタリング成分ＦＴ［０，ｆ］によって、入力時系列ＲＩの全画素を含む３次元画素テーブルＰＸにおける画素ＰＸ［ｉ，ｊ，ｔ］の近傍で選択された近隣画素（図３の画像ＲＩ［ｔ］およびＲＩ［ｔ＋１］内の黒色画素）の線形結合を行うことによって計算される。線形結合の結果にバイアス係数が加算されてよく、加算結果は、正規化線形ユニット（ＲｅＬＵ）関数によって変換されてよく、この関数は「漏洩」型であってよい。ＲｅＬＵ関数は、ステージからの合計加重入力をステージのアクティブ化またはその入力に関する出力に変換する役割を担う活性化関数である。一般に、ＲｅＬＵ関数は、入力が正であれば入力を直接出力し、正でなければゼロを出力する区分的線形関数である。この目的は、ディープネットワークが効果的に学習することを妨げ得る「勾配消失問題」として知られる問題を解消することである。他の活性化関数が利用できる。漏洩ＲｅＬＵ関数は、線形減衰係数を有する負の値を維持する。

実施形態によると、各フィルタ画像ＦＩ［０，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］（ｉ＝１、・・・Ｘ、ｊ＝１、・・・Ｙ、ｔ＝１、・・・Ｔ、ｆ＝１、・・・Ｆ）は、以下の式：

に従って、対応するフィルタリング成分ＦＴ［０，ｆ］によって計算され、式中、ＬＲは、漏洩ＲｅＬＵ関数であり、Ｂ［０］は、各々がフィルタリング成分ＦＴ［０，ｆ］によって用いられる一定値を表すバイアス係数Ｂ［０，ｆ］の１次元テーブルであり、Ｗ［０］は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｕ、１・・・Ｆ）を有する加重係数Ｗ［０，ｌ，ｍ，ｕ，ｆ］の４次元テーブルであり、加重係数Ｗ［０，ｌ，ｍ，ｕ，ｆ］は実数の値を有しており、ａ、ｂ、ｃは、線形結合に関与する近隣画素をテーブルＰＸ内で選択するためにパラメータＬ、Ｍ、およびＵで調整され得る係数である。画像ＲＩ［ｔ］外の画素、すなわちｉ＋ａ．ｌ＞Ｘおよび／またはｊ＋ｂ．ｍ＞Ｙおよび／またはｔ＋ｃ．ｕ＞Ｔである画素は、０に設定される。

例によると、ｘ≧０の場合、ＬＲ（ｘ）＝ｘであり、ｘ＜０の場合、ＬＲ（ｘ）＝０．５ｘである。図３の例において、ＬおよびＭは３に設定され、Ｕは２に設定され、ａおよびｂは３に設定され、ｃは１に設定される。したがって、各フィルタ画像ＦＩ［０，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］は、画像時系列ＲＩ［１・・・Ｔ］の画像ＲＩ［ｔ］およびＲＩ［ｔ＋１］の各々における９個の近隣画素から計算される。

図４は、実施形態に従って、処理ステージＦＰ１、ＦＰ２、ＦＰ３、およびＦＰ５を示す。処理ステージＦＰ１～ＦＰ３、ＦＰ５の各々は、Ｆ個のフィルタリング成分ＦＴ［ｓ，ｆ］を備え、ｓは１、２、３、または５に等しく、ｆは１～Ｆに等しく、各フィルタリング成分ＦＴ［ｓ，ｆ］がフィルタ画像時系列ＦＩ［ｓ－１，ｔ］を処理する。各フィルタリング成分ＦＴ［ｓ，ｆ］は、Ｔ個のフィルタ画像ＦＩ［ｓ，ｔ，ｆ］の時系列を生成する。したがって、各処理ステージＦＰｓ（ｓ：１、２、３、または５）は、Ｔ個のフィルタ画像から成るＦ個の時系列を生成する。

図５は、実施形態に従って、処理ステージＦＰｓ（ｓ：１、２、３、または５）の各フィルタリング成分ＦＴ［ｓ，ｆ］によって行われる動作を示す。各フィルタ画像ＦＩ［ｓ，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］は、対応するフィルタリング成分ＦＴ［ｓ，ｆ］によって、フィルタ画像ＦＩ［ｓ，ｔ，ｆ］のＦ個の時系列の全てにおいて、３次元画素テーブルＰＸ［ｉ，ｊ，ｔ，ｆ］内の画素ＰＸ［ｉ，ｊ，ｔ，ｆ］の近傍で選択された近隣画素（図５の画像ＦＩ［ｓ－１，ｔ，ｆ］およびＦＩ［ｓ－１，ｔ＋１，ｆ］内の黒色画素）の線形結合を行うことによって計算される。ここで再び、線形結合の結果にバイアス係数が加算されてよく、加算結果は、同じ正規化線形ユニット（ＲｅＬＵ）関数によって変換され得る。

実施形態によると、各フィルタ画像ＦＩ［ｓ，ｔ，ｆ’］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ’］（ｆ’は１～Ｆに等しい）は、対応するフィルタリング成分ＦＴ［ｓ，ｆ］によって、以下の式：

に従って計算され、式中、Ｂ［ｓ，ｆ’］は、各々がフィルタリング成分ＦＴ［ｓ，ｆ’］に起因する一定値を有するバイアス係数Ｂ［ｓ，ｆ’］の２次元テーブルであり、Ｗ［ｓ］（ｓ＝１、２、３、または５）は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｕ、１・・・Ｆ、１・・・Ｆ）を有する加重係数Ｗ［ｓ，ｌ，ｍ，ｕ，ｆ，ｆ’］の５次元テーブルであり、加重係数Ｗ［ｓ，ｌ，ｍ，ｕ，ｆ，ｆ’］は、実数の値を有しており、ａ、ｂ、ｃは、線形結合に関与する近隣画素をテーブルＰＸ内で選択するためにパラメータＬ、Ｍ、およびＵで調整され得る係数である。

例によると、ｘ≧０の場合、ＬＲ（ｘ）＝ｘであり、ｘ＜０の場合、ＬＲ（ｘ）＝０．５ｘである。図５の例において、ＬおよびＭは３に設定され、Ｕは２に設定され、ａおよびｂは３に設定され、ｃは１に設定される。したがって、各フィルタ画像ＦＩ［ｓ，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］は、画像セットＦＩ［ｓ－１，ｔ，ｆ］およびＦＩ［ｓ－１，ｔ＋１，ｆ］の全てにおける９個の（画素ＰＸ［ｉ，ｊ］の）近隣画素から計算され、ｆ＝１・・・Ｆであり、すなわち、各フィルタ画像ＦＩ［ｓ，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］は、３×３×２×Ｆの画素から計算される。

図６および図７は、実施形態に従って、プロセッサの処理ステージＳＬ４およびＳＬ６を示す。処理ステージＳＬ４およびＳＬ６の各々は、Ｆ個の組合せ成分ＣＢ［ｓ，ｆ］（ｓは４または６に等しく、ｆは１～Ｆに等しい）を備え、各組合せ成分ＣＢ［ｓ，ｆ］は、フィルタ画像時系列ＦＩ［ｓ－１，ｔ］を処理し、Ｔ／ｎ個のフィルタ画像ＦＩ［ｓ，ｔ，ｆ］から成る時系列を提供し、Ｔはｎの倍数である。組合せ成分ＣＢ［ｓ，ｆ］によって行われる組合せ演算は、フィルタ画像のＦ個の時系列ＦＩ［ｓ－１，ｆ］の各々における画素グループＰＸ［ｉ，ｊ，ｔ］、ＰＸ［ｉ，ｊ，ｔ＋１］、・・・、［ｉ，ｊ，ｔ＋ｋ］の全てをそれぞれ単一の合成画素に組み合わせ、フィルタ画像ＦＩ［ｓ－１，ｆ］の各々における各画素は、時系列ＦＩ［ｓ，ｆ］を計算するために一度だけ考慮され、ｋは、結合成分ＣＢ［ｓ，ｆ］によって組み合わせられたフィルタ画像ＦＩ［ｓ，ｔ，ｆ］の数を表す。

実施形態によると、各フィルタ画像ＦＩ［ｓ，ｔ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｔ，ｆ］（ｆは１～Ｆに等しい）は、３つのフィルタ画像ＦＩ［ｓ－１，ｔ］ＦＩ［ｓ－１，ｔ＋１］、およびＦＩ［ｓ－１，ｔ＋２］から、対応する組合せ成分ＣＢ［ｓ，ｆ］によって、以下の式：

に従って計算され、式中、ＭＡＸは、３つの画素の各グループＰＸ［ｉ，ｊ，３（ｔ－１）＋１］、ＰＸ［ｉ，ｊ，３（ｔ－１）＋２］、およびＰＸ［ｉ，ｊ，３（ｔ－１）＋３］内の最大値を提供する関数であり、ｔは１、２、・・・Ｔ／３に等しい。その結果、処理ステージＳＬ４およびＳＬ６によって提供されたフィルタ画像ＦＩ（ｓ，ｆ）から成る時系列の各々は、Ｔ／３個の画像を備え、Ｔは３の倍数である。

処理ステージＳＬ６の後、フィルタ画像ＦＩ（ｓ，ｆ，ｔ）から成る時系列の各々には１つの画像ＦＩ（ｓ，ｆ）のみが残る。

図８は、実施形態に従って、プロセッサの処理ステージＦＰ７、ＦＰ８、ＦＰ９、ＦＰ１０を示す。処理ステージＦＰ７～１０の各々は、Ｆ個のフィルタリング成分ＦＴ［ｓ，ｆ］を備え、ｓは７、８、９、または１０に等しく、ｆは１～Ｆに等しく、各フィルタリング成分ＦＴ［ｓ，ｆ］は、フィルタ画像ＦＩ［ｓ－１，ｆ］を処理し、１つのフィルタ画像ＦＩ［ｓ，ｆ］を生成する。

図９は、実施形態に従って、処理ステージＦＰ［ｓ］（ｓ＝７、８、９、１０、または１１）の各フィルタリング成分ＦＴ［ｓ，ｆ］によって行われる動作を示す。各フィルタ画像ＦＩ［ｓ，ｆ］の各画素ＰＸ［ｉ，ｊ，ｆ］は、対応するフィルタリング成分ＦＴ［ｓ，ｆ］によって、Ｆ個のフィルタ画像ＦＩ［ｓ，ｆ］の全てにおいて、３次元画素テーブルＰＸ内の画素ＰＸ［ｉ，ｊ，ｆ］の近傍で選択された隣接画素（図９のＦＩ［ｓ－１，ｆ］における黒色画素）の線形結合を行うことによって計算される。ここで再び、線形結合の結果にバイアス係数が加算されてよく、加算結果は、同じ正規化線形ユニット（ＲｅＬＵ）関数によって変換され得る。

実施形態によると、各フィルタ画像ＦＩ［ｓ，ｆ’］（ｆ’は１～Ｆに等しい）の各画素ＰＸ［ｉ，ｊ，ｆ’］は、対応するフィルタリング成分ＦＴ［ｓ，ｆ］によって、以下の式：

に従って計算され、式中、Ｂ［ｓ］は、各々がフィルタリング成分ＦＴ［ｓ，ｆ’］に起因する一定値を有するバイアス係数Ｂ［ｓ，ｆ’］の１次元テーブルであり、Ｗ［ｓ］は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｆ、１・・・Ｆ）を有する、ｓ＝７、８、９、１０、または１１の加重係数Ｗ［ｓ，ｌ，ｍ，ｆ，ｆ’］の４次元テーブルであり、加重係数Ｗ［ｓ，ｌ，ｍ，ｆ，ｆ’］は、実数の値を有しており、ａおよびｂは、線形結合に関与する近隣画素をテーブルＰＸ内で選択するためにパラメータＬおよびＭで調整され得る係数である。図９の例において、ＬおよびＭは３に設定され、ａおよびｂは３に設定される。したがって、フィルタ画像ＦＩ［ｓ，ｆ］の各々の各画素ＰＸ［ｉ，ｊ］は、画像ＦＩ［ｓ－１，ｆ］（ｆ＝１～Ｆ）の全てにおける９つの（位置［ｉ，ｊ］と）近隣画素から計算される。

画像時系列ＲＩ［ｔ］によって撮像された地面領域の標高モデルの画像ＥＭが、インデックスＦ＋１を有するフィルタ画像として処理ステージＦＰ７に導入される。したがって、処理ステージＦＰ７で行われる計算において、１とＦ＋１との間で式（４）におけるインデックスｆの集計が行われ、４次元テーブルＷ［７］は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｆ＋１、１・・・Ｆ）を有する。処理ステージＦＰ８～１１において、４次元テーブルＷ［ｓ］は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｆ、１・・・Ｆ）を有する。

図１０は、実施形態に従って、プロセッサの最終処理または出力ステージＦＰ１２を示す。出力ステージＦＰ１２は、フィルタ画像ＦＩ［１１，ｆ］を処理して出力画像ＰＩを生成する１つのフィルタリング成分ＦＴ［１２］を備える。

図１１は、実施形態に従って、出力ステージＦＰ１２のフィルタリング成分ＦＴ［１２］によって行われる動作を示す。出力画像ＰＩの各画素ＰＩ［ｉ，ｊ］は、フィルタリング成分ＦＴ［１２］によって、Ｆ個のフィルタ画像ＦＩ［１１，ｆ］の全ての画素を備える３次元画素テーブルＰＸ［ｉ，ｊ，ｆ］において計算された画素の近傍で選択された近隣画素の線形結合を行うことによって計算される。ここで再び、線形結合の結果にバイアス係数が加算されてよく、加算結果は、同じ正規化線形ユニット（ＲｅＬＵ）関数によって変換され得る。

実施形態によると、出力画像ＰＩの各画素ＰＸ［ｉ，ｊ］は、フィルタリング成分Ｆ［１１］によって、以下の式：

に従って計算され、式中、Ｗ［１２］は、次元（１・・・Ｌ、１・・・Ｍ、１・・・Ｆ）を有する加重係数Ｗ［１２，ｌ，ｍ，ｆ］の３次元テーブルであり、加重係数Ｗ［１２，ｌ，ｍ，ｆ］は、実数の値を有しており、ａおよびｂは、線形結合に関与する近隣画素をテーブルＰＸ内で選択するためにパラメータＬおよびＭで調整され得る係数である。図１１の例において、ＬおよびＭは３に設定され、ａおよびｂは３に設定される。したがって、出力画像ＰＩの各画素ＰＸ［ｉ，ｊ］は、画像ＦＩ［１１，ｆ］（ｆ＝１～Ｆ）の各々における９個の近隣画素から計算される。

ステージＦＰｓ（ｓ＝１～３、５、７～１１）におけるフィルタリング成分ＦＴ［ｓ，ｆ］の数Ｆの選択は、入力画像からノイズを除去するためのニューラルネットワークの効率と計算時間とのトレードオフの結果である。一般に、数Ｆは、１０～１００から選択される。上記例において、Ｆは６４に設定される。したがって、
Ｗ［０］は、２×３×３×６４（＝１１６２）の加重係数を含み、Ｂ［０，ｆ］は、６４のバイアス係数を含み、
ｓ＝１、２、３、および５であるＷ［ｓ］およびＢ［ｓ，ｆ］は、それぞれ２×３×３×６４×６４（＝７３７２８）の加重係数および６４のバイアス係数を含み、
Ｗ［７］は、３×３×６４×６５（＝３７４４０）の加重係数を含み、Ｂ［７，ｆ］は、６４のバイアス係数を含み、
ｓ＝８、９、１０、および１１であるＷ［ｓ］およびＢ［ｓ，ｆ］は、それぞれ３×３×６４×６４（＝３６８６４）の加重係数および６４のバイアス係数を含み、
Ｗ［１２］は、３×３×６４（＝５７６）の加重係数を含み、Ｂ［１２］は、１つのバイアス係数を表す。

実施形態によると、全ての係数ＷおよびＢは、訓練データを用いて学習フェーズによって決定される。実施形態によると、訓練データは、断層運動の単純な弾性モデルを用いて地盤運動のシミュレーションから生成された数百万または数十億の画像時系列を備える。画像時系列は、たとえば、ランダムな緯度、経度、深度、断層走行角度、地平線俯角、および幅で徐々に不活動化する断層をランダムに生成することによって生成される。その後、地盤変動時系列は、様々なノイズ信号で損なわれる。たとえば、空間相関ガウスノイズは、様々な長さスケールの大気乱流によってＩｎＳＡＲ画像に生じた遅延を模擬する。また、様々な追加のスプリアス信号が、地盤変動時系列を更に損なうように用いられ得る。たとえば、地盤運動信号のように見えるが経時的な相関を有さない非常に急峻な天候関連信号の最悪の場合のシナリオを模擬するために、過渡地盤運動信号が追加され得る。また、画素の相関除去を模擬しエラーをアンラップする誤った画素または画素パッチも画像時系列内に生じ得る。

各訓練画像時系列は、ニューラルネットワークＣＮＮによって提供される最終画像に関連付けられる。ニューラルネットワークの訓練は、係数ＷおよびＢが未知数である方程式に対する解を求めることを備え、画像時系列ＲＩおよび結果画像ＰＩの画素ＰＸ［ｉ，ｊ］は、訓練事例から知られている。実施形態によると、解は、モデルコスト関数を用いる反復勾配に基づく降下最小化法を用いて計算され、係数ＷおよびＢの初期値は、たとえば均一分布を有するように－１～１からランダムに選択される。

実施形態によると、用いられるモデルコスト関数ＭＣＦは、変動再構成のエラーノルム、たとえば各画素の再構成エラーの絶対値の和であり、

であり、式中、ＣＰＸ［ｉ，ｊ］は、訓練画像時系列からニューラルネットワークＣＮＮによって提供された画像ＰＩ内の計算された画素を表し、ＭＰＸ［ｉ，ｊ］は、訓練画像時系列ＰＩに関連付けられた、提供される最終画像内の対応する画素を表す。係数ＷおよびＢは、訓練データのグループにおいてモデルコスト関数ＭＣＦの結果を最小にするために反復的に調整される。凸状ではないコスト関数を用いる場合、大域最小値を満たすことはなく、初期値に依存する極小値のみを満たすことがほぼ確実である。したがって、訓練データに依存して係数ＷおよびＢの唯一のセットに到達することが確実である。

実施形態によると、係数ＷおよびＢは、以下の式：

に従ってＡＤＡＭ規則に則りバッチデータから計算されるコスト関数の勾配に基づいて反復ごとに更新され、式中、
ｗは、係数ＷおよびＢのうち１つであり、ｗｓｔｐおよびεは、それぞれ０．００１および１０^－６に設定された一定値であり、

であり、
ｍ＝β_１ｍ＋ｇ（１－β_１）、ｖ＝β_２ｖ＋ｇ^２（１－β_２）であり、ｍおよびｖは最初に０に設定され、β_１およびβ_２は、それぞれ０．９および０．９９９に設定された一定値であり、ｔは反復回数であり、
ｇ＝勾配（ｘ，ｙ）であり、ｘは、１度の反復で用いられる時系列のセットであり、ｙは、対応する真の出力変動であり、ｇは、モデルパラメータの関数としてコスト関数の勾配の平均値に設定可能である。

モデルの係数ＷおよびＢを調整するために他の勾配降下アルゴリズムが同様に用いられ得ることが当業者には明らかである。

完全な畳み込みニューラルネットワークの利点は、処理される画像のサイズが訓練画像のサイズに依存しない点であり、上記の式（１）～（５）は全て単一の画素を計算し、画像の各画素を計算するために用いられるので、これらの画像の画素数に依存しない。

実施形態によると、ニューラルネットワークは、モデル化地盤変動およびモデル化ノイズから生成された４０×４０の画素から成る９つの画像を備える合成訓練時系列を用いて訓練される。

訓練された畳み込みネットワークＣＮＮの効率を試験するために、専門家によって既に「手動で」分析されている２つの領域、トルコの北アナトリアン断層（ＣＯＳＭＯ－ＳｋｙＭｅｄのデータ）、およびアフガニスタンとパキスタンの国境にあるチャマン断層（Ｓｅｎｔｉｎｅｌ １のデータ）から得たＩｎＳＡＲ画像時系列を処理することが用いられる。チャマンから得たＩｎＳＡＲ画像は、おおよそ１８００００ｋｍ^２の面積をカバーする７０２４×２４８８の画素から成る９個の画像の時系列を備える。トルコおよびチャマンの両方において、従来の大気補正の後でも、断層運動の専門家によって識別された信号よりも強い信号がデータに残っている。両方の事例において、断層変動時系列を生成するためには、断層位置の先験的知識およびデータの手動検査が必要である。合成訓練データを用いて訓練され、実際のデータに追加の微調整を加えていないニューラルネットワークＣＮＮは、専門家の分析でも信号が発見されたトルコおよびチャマンにおけるクリーンな変動信号を自動的に分離および復元する。北アナトリアン断層の時系列において、ニューラルネットワークＣＮＮは１．５ｃｍの視界滑り線を発見し、人手を介さず、また断層位置の知識を有することなく、実質的に断層以外のどこにもノイズ信号は残っていなかった。

チャマンから得た画像時系列から地盤変動を抽出するためにニューラルネットワークＣＮＮを実行する時に行われる計算は全て、従来のパーソナルコンピュータのグラフィックプロセッサ上で行うことが可能であり、最小限の計算時間、たとえば７０２４×２４８８の画素から成る９個の取得物の時系列から地盤変動を抽出するために、単一のＮｖｉｄｉａ（登録商標）ＲＴＸ６０００グラフィックプロセッサにおいて約４分を要する。

図１２は、従来のパーソナルコンピュータの構成要素を示す。コンピュータＰＣは、少なくとも１つのプロセッサＰＲＣと、プロセッサに動作可能に結合された、ユーザインタフェースＤＳＰ、ＣＭ、メモリＭＥＭ、グラフィックプロセッサＧＰ、および通信回路ＮＩＴとを備える。メモリＭＥＭは、オペレーティングシステムおよびアプリケーションを格納する。ユーザインタフェースＤＳＰ、ＣＭは、限定はされないがたとえばキーボードＣＭ、およびグラフィックカードＧＰに接続されたコンピュータディスプレイスクリーンＤＳＰを含んでよい。通信回路ＮＩＴは、プロセッサＰＲＣが、たとえばインターネットなどの電子通信ネットワークＮＴに動作可能に結合されることを可能にする。メモリＭＥＭは、ニューラルネットワークＣＮＮを実装するコンピュータ可読プログラムコード構成要素も含む。たとえば、これらのコンピュータ可読プログラムコード構成要素がプロセッサＰＲＣによって処理されると、プログラムコード構成要素は、上述のような画像時系列を処理するための方法を実行させるように構成される。

本発明の様々な実施形態に関する上記説明は、当業者への説明を目的として提供される。上記説明は、包括的であること、または本発明を単一の開示された実施形態に限定することを意図されているものではない。上記教示の当業者には、本発明に対する数々の代替例および変形例が明らかである。したがって、いくつかの代替実施形態が具体的に説明されているが、当業者には他の実施形態が明らかであり、または比較的容易に開発される。

この点に関して、図１に開示する処理ステージによって行われる動作を他の順序で行い、処理ステージの数が変更され得ることが当業者には明らかである。特に、処理ステージＳＬ４およびＳＬ６は、他のフィルタリングステージＦＰの間に配置され得る。処理ステージＳＬ４、ＳＬ６によって行われる組合せ演算ＣＢ［ｓ，ｆ］は、たとえば加算機能や平均化機能を用いて、最大画素値の選択とは異なる多数の方法で行われてよい。処理ステージＳＬの数および組合せ成分ＣＢ［ｓ，ｆ］の各々によって処理される画像の数は、ニューラルネットワークＣＮＮにおける最後の処理ステージＳＬの組合せ成分ＣＢ［ｓ，ｆ］の各々が時系列ごとに単一の画像ＦＩ［ｓ，ｆ］を提供することを前提として、入力時系列ＲＩ内の画像の数の関数として設定される。加えて、時間成分（インデックスｔ）がもう存在しない他のステージにモデル画像ＥＭが挿入されてよく、モデル画像ＥＭの導入の前および後のフィルタリングステージＦＰの数もまた、処理する画像の種類、特に処理する画像に影響を及ぼすノイズ信号の種類に依存して変更され得る。

加えて、線形結合を計算するために考慮される時系列内の近隣画素の数およびそれぞれの位置および隣接画像の数もまた、処理する画像の種類および処理する画像に影響を及ぼすノイズ信号の種類に依存して変更され得る。線形結合ＦＴ［ｓ，ｆ］に関与する近隣画素を選択するパラメータＬ、Ｍ、およびＵおよび係数ａ、ｂ、ｃの調整は、処理される画像に影響を及ぼす摂動の振幅の関数として行われ得る。一般に、これらのパラメータおよび係数は、１～４から選択される。

フィルタリング演算ＦＰにおけるバイアス係数Ｂ［ｓ，ｆ］の追加も任意選択であるが、ニューラルネットワークＣＮＮの設計に他の自由度をもたらす。このバイアス係数の除去は、膨大な数の加重係数Ｗを導入することにより多数の自由度をもたらす１つ以上のフィルタリング演算ＦＰを追加することによって、自由度に関して容易に補償され得る。

正規化線形ユニット関数ＲＬの適用も任意選択であり、処理される画像の種類に依存する。上記例において、特に処理される画像が大気乱流によって損なわれたＩｎＳＡＲ画像である場合、画素の負の値よりも正の値をより重視することが望ましい。他の応用において、他の正規化線形ユニット関数が適当であり得る。加えて、そのような関数を用いる代わりに、様々な加重係数Ｗが、加重係数と乗算される画素値の関数として定義され得る。

加えて、ニューラルネットワークＣＮＮの係数（加重係数Ｗおよびバイアス係数Ｂ）を計算するための訓練フェーズは、これらの係数が、処理する画像および信号の種類のみに依存することにより、各真正画像の処理より前に行われる必要はない。したがって、１種類の画像を処理するように設計されたニューラルネットワークを実装するコンピュータは、訓練方法および訓練データの生成方法を実施する必要はない。

また、上述の方法は、各観測衛星によって取得された地上のＩｎＳＡＲ画像のみに適用されるものではない。この方法は、時系列の画像間の運動信号を抽出するために、ノイズモデルが既知であるノイズの多い画像時系列を処理するために容易に適合され得る。

上記説明は、本明細書で説明されている本発明の全ての代替、修正、および変形例と、上記説明の主旨および範囲に収まる他の実施形態とを包括することが意図されている。特許請求の範囲に記載の範囲は、特許請求の範囲で用いられる言語に基づいて幅広く解釈すべきであり、そのような範囲は、本明細書で説明する特定の例に限定されてはならない。

Claims (15)

1. ノイズにさらされる、同じ領域の画像の時系列（ＲＩ［ｔ］）を処理するための方法であって、前記方法は、
   入力画像時系列（ＲＩ［ｔ］）を受信する第１のフィルタリングステップ（ＦＰ０）と、出力画像（ＰＩ）を提供する最後のフィルタリングステップ（ＦＰ１２）と、中間フィルタリングステップ（ＦＰｓ、ｓ＝１～３、５、７～１１）とを備える複数の連続的なフィルタリングステップ（ＦＰｓ、ｓ＝０～３、５、７～１２）を実行することであって、前記第１のフィルタリングステップおよび前記中間フィルタリングステップの各々は、最初に前記入力画像時系列を含むフィルタ画像（ＦＩ［ｓ，ｔ，ｆ］）の時系列セット（ＦＩ）を変換し、前記セットは、前記セットの各時系列の各画像の各画素（ＰＸ［ｉ，ｊ，ｔ］）を前記セットの画像時系列の画像および隣接画像における選択された近隣画素と線形結合（ＦＴ［ｓ，ｔ，ｆ］）によって組み合わせることによって、複数（Ｆ、１）のフィルタ画像時系列を生成することによって変換され、各線形結合は、それぞれの加重係数セットを用いて、前記生成されたフィルタ画像のうち１つの画素（ＰＸ［ｉ，ｊ，ｔ，ｆ］）となることと、
   １つ以上の組合せ演算（ＳＬ４、ＳＬ６）を実行することであって、各組合せ演算は、前記セットのフィルタ画像時系列の各々における隣接画像サブセットの画像を組み合わせることによって前記セットのフィルタ画像時系列の各々における画像の数を低減するために、２つの連続する中間フィルタリングステップの間に行われ、前記組合せ演算のうち最後の１つ（ＳＬ６）は、前記セットの各フィルタ画像時系列を単一のフィルタ画像に低減することと、
   を備える、方法。
2. 最後の組合せ演算（ＳＬ６）の後、前記セット（ＦＩ）内のエリアのモデル画像（ＥＭ）を追加のフィルタ画像として導入することを更に備え、前記モデル画像の導入の後に、前記フィルタリングステップ（ＦＰｓ、ｓ＝７～１２）のうち１つ以上が後続する、請求項１に記載の方法。
3. 前記線形結合（ＦＴ［ｓ，ｔ，ｆ］）の各々は、前記線形結合の結果に加算されるバイアス係数を含む、請求項１または２に記載の方法。
4. 前記加重係数の各々は、乗算される画素値の正または負の符号に依存する値を有し、または
   前記線形結合（ＦＴ［ｓ，ｔ，ｆ］）の各々の結果は、正規化線形ユニット関数によって変換される、
   請求項１乃至３の１項に記載の方法。
5. ランダムに選択されたパラメータを用いて地盤変動モデルから第１の画像時系列を生成することと、
   様々なノイズ信号のモデルを用いて、前記第１の画像時系列から訓練画像時系列を生成することと、
   前記加重係数の値を調整するために前記生成された訓練時系列を用いることと、
   を更に備える、請求項１乃至４の１項に記載の方法。
6. 前記加重係数の値は、前記訓練画像時系列およびモデルコスト関数（ＭＥＢ）を用いて前記モデルコスト関数の結果を最小にするように反復勾配に基づく降下最小法を適用することによって反復的に調整される、請求項５に記載の方法。
7. 前記入力画像時系列（ＲＩ［ｔ］）から前記フィルタ画像時系列のセット（ＦＩ［０，ｔ，ｆ］）を生成することは、以下の式：
   

   

   を用いて行われ、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ］は、フィルタ画像時系列ｆのフィルタ画像のうち１つの１画素であり、ＰＸ［ｉ，ｊ，ｔ］は、前記入力画像時系列の１つの画像の１画素であり、Ｗ［ｌ，ｍ，ｕ，ｆ］は、前記フィルタ画像時系列ｆに関する第１の係数のうち１つであり、Ｂ［ｆ］は、前記時系列ｆに関するバイアス係数であり、ＬＲ（）は、正規化線形ユニット関数である、請求項１乃至６の１項に記載の方法。
8. 前記フィルタリング演算の第１のフィルタリング演算（ＦＰｓ、ｓ＝１～３、５）の各線形結合は、以下の式：
   

   

   を適用し、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ’］は、フィルタ画像時系列ｆ’のフィルタ画像のうち１つの１画素であり、Ｗ［ｓ，ｌ，ｍ，ｕ，ｆ，ｆ’］は、前記フィルタリング演算ｓに関するフィルタ画像時系列ｆ’に関する第２の加重係数のうち１つであり、Ｂ［ｓ，ｆ’］は、前記時系列ｆ’および前記フィルタリング演算ｓに関するバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である、請求項１乃至７の１項に記載の方法。
9. 前記フィルタリング演算の第２のフィルタリング演算（ＦＰｓ、ｓ＝７～１１）の各線形結合は、以下の式：
   

   

   を適用し、式中、ＰＸ［ｉ，ｊ，ｆ’］は、フィルタ画像ｆ’のうち１つの１画素であり、Ｗ［ｓ，ｌ，ｍ，ｆ，ｆ’］は、フィルタリング演算ｓに関する前記フィルタ画像ｆ’に関する第２の加重係数のうち１つであり、Ｂ［ｓ，ｆ’］は、前記フィルタ画像ｆ’および前記フィルタリング演算ｓに関するバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である、請求項１乃至８の１項に記載の方法。
10. 前記出力画像（ＰＩ）の各画素（ＰＸ［ｉ，ｊ］）は、以下の式：
    

    

    によって計算され、式中、ＰＸ［ｉ，ｊ］は、前記出力画像の１画素であり、ＰＸ［ｉ，ｊ，ｆ］は、フィルタ画像ｆの１画素であり、Ｗ［ｌ，ｍ，ｆ］は、第３の係数のうち１つであり、Ｂはバイアス係数であり、ＬＲ（）は、漏洩正規化線形ユニット関数である、請求項９に記載の方法。
11. 前記漏洩正規化線形ユニット関数ＬＲは、ｘ≧０の場合、ＬＲ（ｘ）＝ｘであり、ｘ＜０の場合、ＬＲ（ｘ）＝０．５ｘとなる、請求項７乃至１０の１項に記載の方法。
12. 前記画像組合せ演算（ＳＬｓ、ｓ＝４、６）の各々は、以下の式：
    

    

    を適用し、式中、ＰＸ［ｉ，ｊ，ｔ，ｆ］は、前記フィルタ画像時系列ｆのフィルタ画像のうち１つの１画素であり、ＭＡＸ（ＰＸ［ｉ，ｊ，ｔ＋ｕ，ｆ］）は、ｕ＝－１、０、および１を有する画素値ＰＸ［ｉ，ｊ，ｔ＋ｕ，ｆ］の中から最大値を提供する関数である、請求項１乃至１１の１項に記載の方法。
13. 請求項１乃至１２の１項に記載の方法を実行するように構成されたプロセッサ（ＰＲＣ）と、メモリ（ＭＥＭ）とを備えるコンピュータ。
14. 前記プロセッサ（ＰＲＣ）によって制御されるグラフィックプロセッサ（ＧＰ）を更に備え、前記プロセッサは、前記方法の演算のうちいくつかを実行するように前記グラフィックプロセッサを構成するように構成される、請求項１３に記載のコンピュータ。
15. コンピュータによって実行されるときに、請求項１乃至１２の１項に記載の方法を実行するように前記コンピュータを構成するコード部分を備える、コンピュータメモリにロード可能なコンピュータプログラム製品。

Images