線形地震逆解析

逆モデリングは、与えられた地震記録を生成した地球領域の地下の物理的特性を決定することを目的とする数学的手法である。CookeとSchneider（1983）^{[1]は、これを観測された}地震データから地球の構造と物理的パラメータを計算することと定義した。この手法の根底にある前提は、収集された地震データが、逆解析アルゴリズム^[2]から計算された断面と一致する地球構造からのものであるということである。逆解析の対象となる一般的な地球特性には、音速、地層と流体の密度、音響インピーダンス、ポアソン比、地層の圧縮率、せん断剛性、多孔度、流体飽和度などがある。

この手法は地球物理学者にとって長年有用であり、大きく分けて2つのタイプに分類できます。^[3] 決定論的逆解析と確率論的逆解析です。決定論的逆解析は、地球モデルの出力と観測されたフィールドデータを比較し、地球モデルパラメータを継続的に更新することで、通常、モデル出力とフィールド観測値との間の何らかの差異を最小化する関数を最小化する手法です。線形逆解析もこの逆解析手法に含まれますが、この逆解析手法は最小化問題として提示され、受け入れられる地球モデルとは、収集されたフィールド地震データと最もよく一致する数値地震記録を生成するための目的関数を最小化するモデルパラメータの集合です。

一方、確率論的インバージョン法は、クリギングなどの地統計学的ツールを用いて、貯留層流動シミュレーションで使用される制約モデルを生成するために使用されます。単一のモデルパラメータセットを生成する決定論的インバージョン法とは対照的に、確率論的インバージョン法は、モデルの制約に従う一連の代替地球モデルパラメータを生成します。しかし、決定論的モデルの結果は、確率論的インバージョン法の可能な非一意解の平均であるため、この2つの手法は関連しています。^[3]地震線形インバージョンは決定論的インバージョン法であるため、確率論的インバージョン法についてはこれ以上説明しません。

図1：線形地震インバージョンフローチャート

線形反転

線形逆変換の決定論的性質により、逆変換する地震変数を地球モデルパラメータでモデル化する関数関係が必要になります。この機能関係は物理学の基本法則から導かれる数学モデルであり、多くの場合、順方向モデルと呼ばれます。この手法の目的は、順方向モデルとソースウェーブレットの畳み込みと現場で収集された地震トレースの差に依存する関数を最小化することです。最適化の分野と同様に、この最小化される関数は目的関数と呼ばれ、従来の逆モデリングでは、畳み込まれた順方向モデルと地震トレースの差に過ぎません。前述のように、さまざまな種類の変数を逆変換できますが、わかりやすくするために、これらの変数は地球モデルのインピーダンス級数と呼ばれます。次のサブセクションでは、最小化問題としての線形逆変換のコンテキストで、地震データを逆変換するために必要なさまざまなコンポーネントを詳細に説明します。

フォワードモデル

地震線形インバージョンの中核となるのは、収集された実験データの生成をモデル化する順方向モデルです。^{[1] Wiggins (1972)、 [4]}によると、このモデルは、モデルパラメータと観測されたトレースの計算値との間の機能的 (計算的) 関係を提供します。収集された地震データに応じて、このモデルは、岩石または流体中の音波伝播の粒子変位または流体圧力を予測する古典的な波動方程式から、これらの古典的な方程式のいくつかの変形まで変化する可能性があります。たとえば、Tarantola (1984) ^{[5]の順方向モデルは、}地震波伝播中の液体媒体における圧力変化の波動方程式であり、Kanasewich and Chiu (1985) ^[6]は、平面インターフェイスを持つ一定速度層を仮定して、パスに沿った光線の移動時間に対して John Bernoulli の brachistotrone モデルを使用しました。Cooke and Schneider (1983) ^[1]では、モデルは式で表現される合成トレース生成アルゴリズムです。 3、ここでR(t)はZ領域で再帰式によって生成される。順モデルがどのような形で現れるにせよ、収集されたフィールドデータを予測するだけでなく、データがどのように生成されるかをモデル化することが重要である。したがって、CookeとSchneider (1983) ^[1]による順モデルは、平面波源に対する 横方向に均質な地球の応答を模倣することで、拡散損失がないと常に仮定しているため、CMPデータの逆解析にのみ使用できる。

1. ${\displaystyle t=\sum _{i=1}^{n}{\frac {{\big [}(x_{i}-x_{i-1})^{2}+(y_{i}-y_{i-1})^{2}+(z_{i}-z_{i-1})^{2}{\big ]}^{\frac {1}{2}}}{v_{i}}}}$
ここで、t は光線の移動時間、x、y、z は深さ座標、vi は界面 i − 1 と i の間の一定速度です。
2. ${\displaystyle \left[{\frac {1}{K({\vec {r}})}}{\frac {\partial ^{2}}{\partial t^{2}}}-\nabla \cdot {\big (}{\frac {1}{\rho ({\vec {r}}{\big )}}}\nabla )\right]U({\vec {r}},t)=s({\vec {r}},t)}$
ここで、体積弾性率、密度、音波の発生源、および 圧力変化を表します。 ${\displaystyle K({\vec {r}})}$ ${\displaystyle \rho ({\vec {r}})}$ ${\displaystyle s({\vec {r}},t)}$ ${\displaystyle U({\vec {r}},t)}$
3. ${\displaystyle s(t)=w(t)*R(t)}$

ここで、s ( t ) = 合成トレース、w ( t ) = ソースウェーブレット、R ( t ) = 反射率関数です。

目的関数

逆モデリングにおける重要な数値処理は、目的関数を最小化することです。目的関数とは、現場で収集された地震データと数値計算された地震データの差によって定義される関数です。古典的な目的関数としては、最小二乗法における実験データと数値データ間の偏差の二乗和、現場データと数値データ間の差の大きさの合計、あるいはこれらの定義の変形などが挙げられます。どのような定義を用いるかに関わらず、逆問題の数値解は、目的関数を最小化する地球モデルとして得られます。

目的関数に加えて、既知のモデルパラメータや地球の一部の地域における既知の層界面といった制約条件も、逆モデリング手順に組み込まれます。Francis 2006 ^[3]によれば、これらの制約条件は、逆解析データには含まれていない事前情報を提供することで、逆解析解の非一意性を低減するのに役立ちます。一方、Cooke and Schneider (1983) ^[1]は、ノイズの制御や地球物理学的によく知られている地域での作業において、これらの制約条件が有用であると報告しています。

一般化線形逆変換手順の数学的解析

逆モデリングの数学的解析の目的は、前節で述べたすべての要素、すなわち順モデル、目的関数などを考慮して、一般化線形逆問題を単純な行列代数に落とし込むことです。一般的に、数値的に生成された地震データは、地球モデルパラメータの非線形関数です。この非線形性を除去し、線形代数の概念を適用するための基盤を構築するために、順モデルはテイラー級数を用いた展開によって線形化されます。詳細については、Wiggins (1972)、^[4]、Cooke and Schneider (1983) ^{[1]を参照してください。}

地震観測点（ ）と、逆変換対象となる地球モデルパラメータ（ ）の集合を考えます。観測点は または で表すことができます。ここで、 とはモデルパラメータのベクトル表現であり、地球パラメータの関数としての観測点です。同様に、モデルパラメータの推定値を表す は、1.3節の順モデルを用いて数値的に計算された地震データのベクトルです。 のテイラー級数展開は、以下で与えられます。 ${\displaystyle m}$ ${\displaystyle F_{j}}$ ${\displaystyle j=1,\ldots ,m}$ ${\displaystyle n}$ ${\displaystyle p_{i}}$ ${\displaystyle i=1,\ldots ,n}$ ${\displaystyle {\vec {F}}\,({\vec {p}})}$ ${\displaystyle F_{j}\,(p_{i})}$ ${\displaystyle {\vec {p}}}$ ${\displaystyle {\vec {F}}\,({\vec {p}})}$ ${\displaystyle q_{i}}$ ${\displaystyle {\vec {F}}\,({\vec {q}})}$ ${\displaystyle {\vec {F}}\,({\vec {p}})}$ ${\displaystyle {\vec {q}}}$

4. ${\displaystyle {\vec {F}}\,({\vec {p}})={\vec {F}}\,({\vec {q}})+({\vec {p}}-{\vec {q}}){\frac {\partial {\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}}}+({\vec {p}}-{\vec {q}})^{2}{\frac {\partial ^{2}{\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}^{2}}}+O({\vec {p}}-{\vec {q}})^{3}}$
非線形項（(p⃗ − ⃗q)の次数が2以上の項）を除いた線形化により、式は次のようになる。
5. ${\displaystyle {\vec {F}}\,({\vec {p}})-{\vec {F}}\,({\vec {q}})=({\vec {p}}-{\vec {q}}){\frac {\partial {\vec {F}}\,({\vec {q}})}{\partial {\vec {p}}}}}$
には要素があり、には要素がある ことを考慮すると、式 5 の離散形式は、以下に示す行列形式の変数 の線形方程式のシステムになります。${\displaystyle {\vec {F}}}$ ${\displaystyle m}$ ${\displaystyle {\vec {p}}}$ ${\displaystyle {\vec {q}}}$ ${\displaystyle n}$ ${\displaystyle m}$ ${\displaystyle n}$
6. ${\displaystyle \Delta {\vec {F}}=\mathbf {A} \,\Delta {\vec {p}}}$
7. ${\displaystyle \Delta {\vec {F}}={\begin{bmatrix}F_{1}({\vec {p}})-F_{1}({\vec {q}})\\\vdots \\F_{m}({\vec {p}})-F_{m}({\vec {q}})\end{bmatrix}}}$
8. ${\displaystyle \Delta {\vec {p}}={\vec {p}}-{\vec {q}}={\begin{bmatrix}p_{1}-q_{1}\\\vdots \\p_{n}-q_{n}\end{bmatrix}}}$
9. ${\displaystyle \mathbf {A} ={\begin{bmatrix}{\frac {\partial F_{1}({\vec {q}})}{\partial p_{1}}}&{\frac {\partial F_{1}({\vec {q}})}{\partial p_{2}}}&\cdots &{\frac {\partial F_{1}({\vec {q}})}{\partial p_{n}}}\\{\frac {\partial F_{2}({\vec {q}})}{\partial p_{1}}}&\cdots &{\frac {\partial F_{2}({\vec {q}})}{\partial p_{n-1}}}&{\frac {\partial F_{2}({\vec {q}})}{\partial p_{n}}}\\\vdots &{\frac {\partial F_{j}({\vec {q}})}{\partial p_{i}}}&\vdots &\vdots \\{\frac {\partial F_{m}({\vec {q}})}{\partial p_{1}}}&{\frac {\partial F_{m}({\vec {q}})}{\partial p_{2}}}&\cdots &{\frac {\partial F_{m}({\vec {q}})}{\partial p_{n}}}\\\end{bmatrix}}}$

${\displaystyle \Delta {\vec {F}}}$ はCookeとSchneider（1983）において差分ベクトルと呼ばれている。 ^[1]サイズは で、その成分は観測されたトレースと数値計算された地震データの差である。はサイズ の補正ベクトルであり、は感度行列と呼ばれる。 サイズは で、各列は未知の地球モデルパラメータの1つに関する順関数の成分の偏微分である、という注釈がある。同様に、各行は数値計算された地震トレースの成分の全ての未知のモデルパラメータに関する偏微分である。 ${\displaystyle m\times 1}$ ${\displaystyle \Delta {\vec {p}}}$ ${\displaystyle n\times 1}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle m\times n}$

解決アルゴリズム

${\displaystyle {\vec {F}}\,({\vec {q}})}$ は順モデルから計算され、 は実験データです。したがって、は既知の特性です。一方、は未知であり、式10を解くことで得られます。この式は、が逆行列である場合、つまり観測値の数が未知の地球パラメータの数と等しい正方行列である場合にのみ、理論的に解くことができます。この場合、未知の補正ベクトル は、以下に示すように、一連の線形方程式を解くための古典的な直接法または反復法のいずれかを使用して解かれます。 ${\displaystyle {\vec {F}}\,({\vec {p}})}$ ${\displaystyle \Delta {\vec {F}}}$ ${\displaystyle \Delta {\vec {p}}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle m}$ ${\displaystyle n}$ ${\displaystyle \Delta {\vec {p}}}$

10. ${\displaystyle \Delta {\vec {p}}=\mathbf {A} ^{-1}\,\Delta {\vec {F}}}$

ほとんどの地震インバージョンアプリケーションでは、インバージョン対象となる地球パラメータの数よりも観測点の数の方が多いため、数学的に過剰決定された連立方程式が導かれます。その結果、式10は理論的に解くことができず、正確な解を得ることができません。^[6]補正ベクトルの推定値は、誤差の二乗和である を最小化する補正ベクトルを求める最小二乗法を用いて得られます。^[6]${\displaystyle m>n}$ ${\displaystyle \Delta {\vec {p}}}$ ${\displaystyle {\vec {e}}\,^{T}{\vec {e}}}$ ${\displaystyle {\vec {e}}}$

誤差は次のように表される。 ${\displaystyle {\vec {e}}}$

11. ${\displaystyle {\vec {e}}=\Delta {\vec {F}}-\mathbf {A} \,\Delta {\vec {p}}}$

最小二乗法では、最小化する補正ベクトルは以下のように得られます。 ${\displaystyle {\vec {e}}\,^{T}{\vec {e}}}$

12. ${\displaystyle {\begin{aligned}\mathbf {A} \,\Delta {\vec {p}}&=\Delta {\vec {F}}\\\mathbf {A} ^{T}\mathbf {A} \,\Delta {\vec {p}}&=\mathbf {A} ^{T}\Delta {\vec {F}}\end{aligned}}}$

したがって、

13. ${\displaystyle \Delta {\vec {p}}=(\mathbf {A} ^{T}\mathbf {A} )^{-1}\,\mathbf {A} ^{T}\Delta {\vec {F}}}$

上記の議論から、目的関数は、 または で与えられた の または ノルム、あるいは または で与えられた の のいずれかとし て定義されます。 ${\displaystyle L_{1}}$ ${\displaystyle L_{2}}$ ${\displaystyle \Delta {\vec {p}}}$ ${\displaystyle \sum _{j=0}^{n}|\Delta p_{j}|}$ ${\displaystyle \sum _{j=0}^{n}|\Delta p_{j}|^{2}}$ ${\displaystyle \Delta {\vec {F}}}$ ${\displaystyle \sum _{i=0}^{m}|\Delta F_{i}|}$ ${\displaystyle \sum _{i=0}^{m}|\Delta F_{i}|^{2}}$

上で説明した逆モデリングの数学的理論を使用して、 またはの任意の実験地震データを逆変換する一般的な手順を図 1 に示し、次のように説明します。${\displaystyle m=n}$ ${\displaystyle m>n}$

インバージョン処理を開始するために、モデルインピーダンスの初期推定値が提供されます。順方向モデルはこの初期推定値を用いて合成地震データを計算し、これを観測地震データから差し引いて差分ベクトルを計算します。

1. 逆変換プロセスを開始するために、モデルインピーダンスの初期推定が提供されます。${\displaystyle {\vec {q}}}$
2. 上記のモデルインピーダンスを使用して、順方向モデルによって合成地震データが計算されます。${\displaystyle {\vec {F}}({\vec {q}})}$
3. 差分ベクトルは、実験地震データと合成地震データの差として計算されます。${\displaystyle {\vec {F}}({\vec {p}})-{\vec {F}}({\vec {q}})}$
4. 感度マトリックスは、インピーダンス プロファイルのこの値で計算されます。${\displaystyle \mathbf {A} }$
5. と上記3の差ベクトルを用いて補正ベクトルを計算する。新しいインピーダンスプロファイルは次のように得られる。${\displaystyle \mathbf {A} }$ ${\displaystyle \Delta {\vec {p}}}$
   14. ${\displaystyle {\vec {p}}={\vec {q}}+\Delta {\vec {p}}}$
6. 計算された補正ベクトルのノルムまたはノルムは、与えられた許容値と比較されます。計算されたノルムが許容値より小さい場合、数値計算手順は終了し、式14から、地球領域の逆インピーダンスプロファイルが次のように与えられます。一方、ノルムが許容値より大きい場合、ステップ2～6の反復が繰り返されますが、式14から計算されたインピーダンスプロファイルが更新されます。図2 ^{[7] は、連続反復プロセスにおけるインピーダンスプロファイル更新の典型的な例を示しています。Cooke and Schneider (1983) [1]}によれば、式14からの補正推定値を反復処理中の新しい初期推定値として使用することで、誤差が減少します。${\displaystyle L_{1}}$ ${\displaystyle L_{2}}$ ${\displaystyle {\vec {p}}}$

地球モデル空間のパラメータ化

インバージョン対象となる変数に関わらず、地球のインピーダンスは深度（地震データの場合は時間）の連続関数であり、数値線形インバージョン手法をこの連続物理モデルに適用するには、連続特性を離散化するか、地球モデルの深度に沿って離散間隔でサンプリングする必要があります。したがって、モデル特性を決定する対象となる全深度は、離散化の出発点として不可欠です。一般的に、図3に示すように、この特性はこの深度にわたって狭い離散間隔でサンプリングされ、地球深度に沿ったインピーダンス変化の高解像度を確保します。アルゴリズムからインバージョンされたインピーダンス値は、離散間隔における平均値を表します。

逆モデリング問題は、特性をサンプリングするための離散区間の数が逆変換対象となるトレースの観測数と等しい場合にのみ理論的に解けることを考慮すると、高解像度のサンプリングは行列が大きくなり、逆変換に非常に多くのコストがかかります。さらに、従属方程式の場合、行列は特異になる可能性があり、ノイズが存在すると逆変換が不安定になる可能性があり、逆変換の対象となる主変数以外のパラメータが必要な場合は、システムが制約不足になる可能性があります。インピーダンス以外の必要なパラメータに関しては、CookeとSchneider (1983) ^[1]は、ソースウェーブレットとスケールファクタを含めることを示しました。

最後に、制約をいくつかのレイヤーまたは離散間隔内の既知のインピーダンス値として扱うことで、解決すべき未知のインピーダンス値の数が減り、逆変換アルゴリズムの結果の精度が向上します。

図8：振幅ログ

図9a: 振幅から反転したインピーダンス対数

図9b: インピーダンスウェルログ

反転の例

Marescot（2010年）の気温逆転

出典: ^[8]

まず、地球のある領域における温度深度分布から地球パラメータ値を逆変換する例から始めます。この例は伝播音波を伴わないため地震波逆変換とは直接関係ありませんが、地震波への応用に移る前に、逆変換技術の実用的な応用を分かりやすく紹介します。この例では、温度センサーを目標深度に設置することで、坑井内の離散的な位置における地球温度を測定します。温度が深度に対して線形分布するという順方向モデルを仮定し、温度深度測定値から2つのパラメータを逆変換します。

順方向モデルは次のように与えられる。

15. ${\displaystyle {\vec {F}}({\vec {q}})={\vec {T}}=a+bz}$

ここで。したがって、 の次元は2 であり、つまり について反転されるパラメータの数は 2 です。 ${\displaystyle {\vec {q}}=[a,b]}$ ${\displaystyle {\vec {q}}}$

この逆変換アルゴリズムの目的は、観測された温度分布と式15の順方向モデルを使用して得られた温度分布との差を最小にする の値 を見つけることです。順方向モデルの次元または温度観測の数を とすると、順方向モデルの要素は次のように表されます。 ${\displaystyle {\vec {p}}}$ ${\displaystyle [a,b]}$ ${\displaystyle n}$

16. ${\displaystyle {\begin{aligned}T_{1}&=a+bz_{1}\\T_{2}&=a+bz_{2}\\\vdots \\T_{n-1}&=a+bz_{n-1}\\T_{n}&=a+bz_{n}\\\end{aligned}}}$
となることによって${\displaystyle {\vec {F}}({\vec {q}})=T}$
17. ${\displaystyle \mathbf {A} ={\begin{bmatrix}1&z_{1}\\1&z_{2}\\\vdots &\vdots \\1&z_{n-1}\\1&z_{n}\\\end{bmatrix}}}$

^{Marescot (2010) [8]}の結果から、深度における観測温度が 、 であった場合の結果を示す。これらの実験データは、地球パラメータ値およびを取得するために逆変換された。より一般的な、多数の温度観測値を伴うケースとして、図4 は、およびの逆変換値を使用して得られた最終的な線形順方向モデルを示している。この図は、実験データと数値データが良好に一致していることを示す。 ${\displaystyle n=2}$ ${\displaystyle T_{1}=19^{\circ }C}$ ${\displaystyle z=2m}$ ${\displaystyle T_{2}=22^{\circ }C}$ ${\displaystyle z=8m}$ ${\displaystyle a=0.5}$ ${\displaystyle b=18^{\circ }C}$ ${\displaystyle a}$ ${\displaystyle b}$

Marescot (2010) による波動伝播時間逆解析

出典: ^[8]

この例では、記録された地震波伝播時間から地層速度を逆解析します。図5は初期速度推定値と現場で記録された伝播時間を示しており、図6aは逆解析アルゴリズムを30回繰り返して得られた解である、逆解析された不均質速度モデルを示しています。図6bに示すように、逆解析された速度を用いた順解析モデルから得られた最終的な伝播時間と現場記録の伝播時間との間には良好な比較が見られます。これらの解を用いて波線経路を再構築したところ、図7に示すように、地球モデル内を非常に曲がりくねった経路であることが示されました。

CookeとSchneider（1983）による地震痕跡逆解析

Cooke and Schneider (1983) ^[1]から引用したこの例は、 CMP地震波形を地球モデルのインピーダンス（密度と速度の積）プロファイルに逆変換した例です。逆変換された地震波形を図8に示し、図9aは逆変換アルゴリズムに使用した初期インピーダンスを入力した逆変換後のインピーダンスプロファイルを示しています。また、図9bに示すように、地震波形と並んで地球領域のインピーダンスログも記録されています。これらの図は、記録されたインピーダンスログと地震波形から逆変換された数値インピーダンスの良好な比較を示しています。

参考文献

1. Cooke, DA; Schneider WA (1983年6月). 「反射地震データの一般化線形インバージョン」.地球物理学. 48 (6): 665– 676. Bibcode :1983Geop...48..665C. doi :10.1190/1.1441497.
2. Pica, A.; Diet JP; Tarantola A. (1990年3月). 「横方向不変媒質における地震反射データの非線形逆解析」.地球物理学. 55 (3): 284– 292. Bibcode :1990Geop...55..284P. doi :10.1190/1.1442836.
3. Francis, AM (2006年11月). 「確率論的および地震学的インバージョンを理解する」. First Break . 24 (11). doi :10.3997/1365-2397.2006026.
4. Wiggins, Ralph (1972年2月). 「一般線形逆問題：表面波と自由振動の地球構造への影響」Reviews of Geophysics . 10 (1): 251– 285. Bibcode :1972RvGSP..10..251W. doi :10.1029/RG010i001p00251.
5. Tarantola, A. (1984). 「地震反射データの線形化と逆解析」.地球物理探査. 32 (6): 908– 1015. doi :10.1111/j.1365-2478.1984.tb00751.x.
6. Kanasewich, ER; Chiu SK L (1985年6月). 「空間地震屈折データの最小二乗法逆変換」アメリカ地震学会誌. 75 (3): 865–880 .
7. Cooke, D; Cant J. (2010年4月). 「モデルベース地震インバージョン：決定論的アプローチと確率論的アプローチの比較」. CSEG Recorder .
8. Marescot, Laurent. 「地球物理学におけるインバージョン入門」 . 2013年5月3日閲覧。

さらに読む

• Backus, G. 1970.「不十分かつ不正確なデータからの推論」米国科学アカデミー紀要65巻1号。
• Backus, G., F. Gilbert. 1968.「地球全体のデータの分解能」王立天文学会地球物理学誌16(2): 169–205.
• Backus, GE, JF Gilbert. 1967.「地球物理学的逆問題に対する形式主義の数値的応用」王立天文学会地球物理学誌. 13 (1–3): 247.
• Bamberger, A., G. Chavent, C. Hemon, P. Lailly. 1982.「垂直入射地震記録の逆解析」地球物理学47(5):757–770.
• Clayton, RW, and RH Stolt. 1981.「音響反射データのためのBorn-WKBJ逆変換法」. Geophysics 46 (11): 1559–1567.
• フランクリン, JN 1970.「不適切線形問題の適切確率的拡張」数学解析応用ジャーナル31(3):682.
• パーカー, RL 1977.「逆理論の理解」地球惑星科学年次レビュー5:35–64。
• ローリンソン、N. 2000.「地殻構造の層状化のための地震データの逆解析」モナッシュ大学博士論文
• Wang, B., LW Braile. 1996.「反射法と屈折法の地震探査データの同時インバージョンとリオグランデ北部リフトの現場データへの適用」Geophysical Journal International 125 (2): 443–458.
• Weglein, AB , HY Zhang, AC Ramirez, F. Liu, JEM Lira. 2009.「地震データの近似線形逆解析の根底にある根本的な意味の解明」地球物理学74(6): 6WCD1–WCD13.