有限差分的時間和空間頻散及其對逆時偏移和全波形反演的影響

任志明， 戴雪， 包乾宗， 蔡曉慧， 劉洋

1 長安大學地質工程與測繪學院， 西安 710064 2 南京工業大學巖土工程研究所， 南京 210009 3 西安郵電大學自動化學院， 西安 710121

0 引言

有限差分法計算量小、容易實現，是地球物理學和地震學領域最常用的波動方程數值解法，也被廣泛應用于逆時偏移和全波形反演中進行波場正反向延拓.但由于時間和空間離散，有限差分法會產生數值頻散，不可避免地降低正演、成像和反演的精度.為改善其精度，不同學者先后提出規則網格和交錯網格(Virieux，1986；Luo and Schuster，1990；董良國等，2000)、顯式和隱式(Liu and Sen，2009a；Kosloff et al.，2010；Chu and Stoffa，2012)、空間域和時空域(Liu and Sen，2019b；Ren and Liu，2015；蔡曉慧等，2015)、非規則網格和無網格(裴正林，2004；孫輝和張劍鋒，2019；Li et al.，2017；Liu et al.，2019；劉立彬等，2020；姜占東等，2021)有限差分.有限差分法在勘探地球物理和計算地震學的詳細進展請參考Carcione等(2002)、劉洋(2014)和Moczo 等(2014)，這里不再贅述.

有限差分法的模擬精度由時間和空間導數的離散方式決定.提高空間精度的方法主要包括：增大空間差分算子的長度；優化空間差分系數(Zhang and Yao，2013；Liu，2014；Yang et al.，2017)；偽譜法(Kosloff and Baysal，1982；Fornberg，1987；Liu and Wei，2005；劉財等，2013；李青陽等，2018).采用長的差分算子逼近空間偏導數會增加計算量，且存在飽和效應(即算子長度繼續增大對模擬精度的改善非常小).常規有限差分法的差分系數是基于泰勒級數展開得到的，只能保證低波數區域的模擬精度.優化有限差分法采用優化算法(如最小二乘、模擬退火等)逼近頻散關系求取差分系數，可以大幅提高中-高波數區域的模擬精度.偽譜法通過正反傅里葉變換求取空間偏導數，理論上具有無窮階/譜精度，但需要更多的計算時間.提高有限差分時間精度的方法包括：Lax-Wendroff法(Dablain，1986；Chen，2009；Amundsen and Pedersen，2017)；新差分模板法(Liu and Sen，2013；Wang et al.，2016；Tang and Huang，2014；張保慶等，2016；Ren et al.，2017)；時間頻散校正法(Stork，2013；Wang and Xu，2015；Koene et al.，2018).Lax-Wendroff法通過空間偏導數替換高階時間偏導數來改善時間精度.新差分模板法采用非“十字”或“正交”模板逼近空間偏導數.Liu和Sen(2013)給出一種“菱形”差分模板，新方法可以獲得任意偶數階時間精度，但計算量較大.為平衡精度與效率，Wang等(2016)提出基于“十字+菱形”模板的有限差分法.Tang和Huang(2014)和Ren等(2017)發展了時間4～6階和任意偶數階交錯網格有限差分法.新差分模板法的差分系數隨速度變化，相對比較復雜.時間頻散校正法通過正反時間頻散變換向觀測記錄中增加或消除時間頻散.Wang和Xu(2015)介紹了有限差分的正反時間頻散變換對，并將其應用于逆時偏移中消除時間頻散對成像的影響.Koene等(2018)指出采用共軛算子代替逆算子的不足，推導了更精確的正時間頻散算子，使得輸入信號經過正反變換后能夠被精確重建.與其他方法相比，時間頻散校正法的效率更高，且可以徹底消除波場中的時間頻散誤差.

波動方程正演是逆時偏移和全波形反演的基本單元，成像和反演的精度依賴于所采用的數值模擬算法.到目前為止，單獨研究有限差分的時間和空間頻散特性及其對逆時偏移和全波形反演的影響的文獻相對較少.本文采用時間有限差分+空間偽譜法和時間頻散校正+空間有限差分分離時間頻散和空間頻散；詳細分析了有限差分的時間和空間頻散特性及其對成像和反演精度的影響；發展了抗時間頻散、抗空間頻散、抗時間+空間頻散的逆時偏移和全波形反演方法；通過理論數據和實際資料對提出的成像和反演方法進行了測試.

1 有限差分的時間和空間頻散

二階常密度聲波方程形式如下(Dablain，1986；Liu and Sen，2009b)：

(1)

式中，L=?2/?x2+?2/?z2，v為傳播速度，p為壓強，s為震源項.常規有限差分時間和空間偏導數離散格式為：

(2)

(3)

其中，cn為差分系數，M為空間算子長度參數，Δt為時間步長，Δx和Δz為空間采樣間隔.方程(2)和(3)兩邊進行傅里葉變換并化解可得：

(4)

(5)

其中，ω為角頻率，k為波數，kx和kz分別為水平與垂直方向的波數.

(6)

(7)

圖1 常規有限差分法離散前后的頻率(a)和波數(b)變化曲線Fig.1 Variation curves of frequency (a) and wavenumber (b) before and after the discretization in the conventional finite difference method

(8)

有限差分法的空間精度隨M的增加而增大.當M趨向于無窮大時，有限差分法變為偽譜法，空間偏導數通過式(9)求解(Kosloff and Baysal，1982)：

Lp=F-1(F(Lp)),

(9)

其中，F和F-1分別表示二維正反傅里葉變換.與空間頻散不同，時間頻散是可預測的，與傳播路徑和速度等因素無關(Koene et al.，2018).有限差分的時間頻散是由于數值離散引起的頻率移動(從ω到ωFD)而產生的.將具有頻散的數據從ωFD校回到真實頻率ω處即可實現時間頻散校正，得到原始的無頻散數據.由方程(6)可得：

(10)

方程(10)定量表示離散前后的頻率變化關系.時間頻散校正步驟為(Wang and Xu，2015)：

(1)通過方程(10)計算不同ω對應的ωFD；

常規有限差分具有時間2階和空間2M階精度，同時產生時間和空間頻散.偽譜法的空間模擬精度為無窮階.采用時間有限差分和偽譜法(方程(2)和(9))求解波動方程，可以得到只包含時間頻散的數據.時間頻散校正法能夠消除常規時間二階有限差分引入的時間頻散.采用時間頻散校正和空間有限差分(方程(3)和(10))可以得到只包含空間頻散的數據.將時間頻散校正和偽譜法(方程(9)和(10))相結合可獲取不含時間和空間頻散的數據.為了書寫方便，下文中將有限差分、偽譜法和時間頻散校正分別記為FD、PS和TDC.因此，FD、PS、FD+TDC和PS+TDC分別表示全頻散、只有時間頻散、只有空間頻散和無頻散數據或者對應的方法.

這里采用一個均勻介質模型研究有限差分的時間和空間頻散特性.介質速度為2000 m·s-1，網格點數為301×301，空間網格大小為10 m×10 m，時間步長為 2.0 ms，最大記錄時間為2.0 s.震源為30 Hz的雷克子波，位于點(1000 m,1000 m)處.圖2給出不同方法得到的點(3000 m,2000 m)處的波形圖及頻譜.其中，Semi-analytical為由Cagniard-De Hoop方法(De Hoop，1960)得到的半解析波場.FD表示時間2階空間10階精度的有限差分.由圖2可見：時間頻散(PS)會引起波場相移，空間頻散(FD+TDC)導致波形震蕩；時間和空間頻散分別出現在正常波前的前部和尾部；PS+TDC可以有效壓制時間和空間頻散，獲得與解析解相近的模擬結果；空間頻散引起頻譜畸變，時間頻散對頻譜的影響相對較小.

圖2 均勻介質模型(3000 m, 2000 m)處不同方法的波形圖(a)及頻譜(b)Fig.2 Waveforms (a) and spectra (b) at the point (3000 m, 2000 m) by different methods on a homogeneous model

2 抗頻散逆時偏移方法

逆時偏移包括震源波場正向傳播、檢波點波場反向傳播和正反向波場互相關成像.將FD、PS、FD+TDC和PS+TDC應用于波場正傳和反傳中，即可研究時間和空間頻散對逆時偏移的影響.但直接進行時間頻散校正，需要存儲所有網格點所有時刻的正反向波場，并進行大量的正反傅里葉變換，計算和存儲需求巨大.為了提高計算效率，Koene等(2018)介紹了一種實用的抗時間頻散逆時偏移(FD+TDC)方法：震源s中加入時間頻散，源波場正向傳播；觀測波場中加入時間頻散，檢波點波場反向傳播；正反向波場互相關成像.向無頻散數據中添加時間頻散就是實現從ω到ωFD的過程，具體步驟為(Koene et al.，2018)：

(1)通過方程(6)計算不同ωFD對應的ω；

震源波場中只包含源-反射點路徑上的時間頻散.觀測波場中加入時間頻散后包含源-反射點-檢波點整個路徑的時間頻散.由于波場正向和反向傳播沿不同的時間方向積累時間頻散，檢波點波場反向傳播會抵消掉反射點-檢波點路徑上積累的時間頻散.因此，震源波場和檢波點波場具有相同的時間頻散(相移)，仍可以實現零延遲互相關成像.

偽譜法具有空間無窮階精度，采用偽譜法進行波場正反向延拓得到抗空間頻散逆時偏移(PS)方法.將時間頻散校正和偽譜法結合得到抗時間+空間頻散逆時偏移(PS+TDC)方法：震源s中加入時間頻散，偽譜法源波場正向傳播；觀測波場中加入時間頻散，偽譜法檢波點波場反向傳播；正反向波場互相關成像.

3 抗頻散全波形反演方法

全波形反演包括震源波場正向傳播、殘差波場反向傳播、模型參數梯度計算和迭代更新.將FD、PS、FD+TDC和PS+TDC應用于波場正傳和反傳中，同樣可以研究時間和空間頻散對全波形反演的影響.抗時間頻散全波形反演(FD+TDC)方法為：震源s中加入時間頻散，源波場正向傳播；觀測波場中加入時間頻散，殘差波場反向傳播；模型參數梯度計算；反演迭代(如共軛梯度法)直至滿足收斂條件.

由于震源波場和伴隨波場中殘留有源-反射點路徑的時間頻散，導致得到的模型參數梯度與無頻散情況仍有差異.模擬波場和觀測波場包含源-反射點-檢波點路徑的時間頻散引起波形震蕩、頻譜改變，也會增加反演的多解性.因此，抗空間頻散(PS)和抗時間+空間頻散(PS+TDC)全波形反演方法仍無法徹底消除時間頻散的影響.與波場延拓相比，觀測波場中添加時間頻散的計算量可忽略不計(且可在偏移和反演之前完成).抗時間+空間頻散(PS+TDC)逆時偏移和全波形反演方法的計算效率與空間頻散校正(PS)方法相近.

4 數值算例

采用salt模型、Sigsbee2A模型、實際資料對抗頻散逆時偏移和全波形反演方法進行測試.通過對比抗空間頻散和抗時間+空間頻散方法的結果研究時間頻散對成像和反演的影響；對比抗時間頻散和抗時間+空間頻散方法的結果研究空間頻散對成像和反演的影響.基于不同方法與參考解的接近程度驗證抗時間+空間頻散方法的有效性.

4.1 逆時偏移

SEG/EAGEsalt模型如圖3所示.空間網格大小為20 m×20 m，時間步長為2.0 ms，最大記錄時間為4.0 s.震源采用主頻為30 Hz的Ricker子波.31個震源和601個檢波點均勻分布于地表.采用偽譜法和小時間步長(Δt=0.5 ms)模擬無頻散觀測記錄.FD表示時間2階空間16階精度的有限差分.常規(FD)和抗空間頻散(PS)逆時偏移中直接對無頻散觀測記錄進行反傳，而抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法中先在觀測記錄中加入時間頻散后再反傳.圖4為第16炮(7800 m, 2160 m)處不同偏移方法的正反向傳播波場和頻譜.為便于比較，采用偽譜法和小時間步長(Δt=0.5 ms)進行波場正反向延拓并成像作為參考解(Ref).由圖4可見，正傳和反傳波場的時間頻散分別引起波形向負時間和正時間方向移動，正傳和反傳波場的空間頻散分別出現在正常波前時間增大和減小的方向；常規(FD)和抗時間頻散(FD+TDC)方法的正傳波場同時包含時間和空間頻散(相移和震蕩)，而抗空間頻散(PS)和抗時間+空間頻散(PS+TDC)方法的正傳波場中只積累了源-反射點路徑的時間頻散(相移)；常規(FD)和抗空間頻散(PS)方法的反傳波場包含檢波點-反射點路徑的時間頻散，而抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法的反傳波場中積累了源-反射點路徑的時間頻散，常規(FD)和抗時間頻散(FD+TDC)方法的反傳波場空間頻散較強；常規(FD)和抗時間頻散(FD+TDC)方法的正傳波場頻譜變形嚴重，抗空間頻散(PS)和抗時間+空間頻散(PS+TDC)方法的正反傳波場頻譜與參考解的差異較小.圖5給出不同方法的正反向傳播波場匹配圖.常規(FD)和抗空間頻散(PS)方法正反傳波場中分別積累了源-反射點和檢波點-反射點路徑的時間頻散，波形沿相反的方向移動，導致正反傳波場無法匹配；抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法的正反傳波場中包含相同的源-反射點路徑時間頻散，波形都沿時間負方向移動，正反傳波場相位相同，仍然能進行互相關條件成像；采用偽譜法(PS)進行正反向波場延拓，可以有效壓制空間頻散，減小波形震蕩對成像精度的影響.

圖3 鹽丘模型Fig.3 Salt dome model

圖4 鹽丘模型(7800 m, 2160 m)處不同方法的正反向波場及頻譜(a) 正向波場； (b) 反向波場； (c) 正向波場頻譜； (d) 反向波場頻譜.Fig.4 Wavefields and spectra at the point (7800 m, 2160 m) by different methods on salt dome model(a) Forward wavefield； (b) Backward wavefield； (c) Spectrum of (a)； (d) Spectrum of (b).

圖5 鹽丘模型不同方法的正反向波場匹配圖(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.5 Matching diagrams of forward- and backward-propagated wavefields by different methods for the salt dome model

圖6為不同偏移方法得到的成像剖面.由圖6可知，常規(FD)方法同時受時間和空間頻散影響，無法準確成像；抗空間頻散(PS)方法只與時間頻散有關，出現同相軸不聚焦、分辨率降低、反射界面偏離真實位置等現象；抗時間頻散(FD+TDC)方法只受空間頻散影響，產生了許多高頻噪聲(如黑色箭頭所示)和虛假反射界面(如白色箭頭所示)；抗時間+空間頻散(PS+TDC)方法幾乎不受頻散影響，得到了與參考解相近的結果.圖7給出不同方法像的波數譜.抗空間頻散(PS)方法的譜中包含垂直方向的線條(如黑色箭頭所示)，對應于空間域的周期性虛假水平同相軸；抗時間頻散(FD+TDC)方法的譜中出現了一些高波數成分(如黑色箭頭所示)，對應于空間域的高頻成像噪聲；抗時間頻散(FD+TDC)方法的波數譜與參考解能夠較好地吻合.

圖6 鹽丘模型不同方法的成像結果(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.6 Images by different methods for the salt dome model

圖7 鹽丘模型不同方法成像結果的波數譜(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.7 Spectra of images by different methods for the salt dome model

下面采用海上某工區實際地震資料對不同的偏移方法進行測試.該數據包含240炮，每炮120道，炮間距50 m，道間距25 m.網格大小為25 m×25 m，時間步長為2.0 ms，最大記錄時間為4.0 s.圖8給出其中5炮的地震記錄.圖9為偏移速度模型.震源子波是通過軟件提取的最小相位子波.FD表示時間2階空間10階精度的有限差分.圖10展示了不同偏移方法的成像結果.由圖10可見：常規(FD)和抗時間頻散(FD+TDC)方法偏移剖面中存在較強的成像噪聲；抗空間頻散(PS)方法偏移剖面同相軸的連續性較差；抗時間+空間頻散(FD+TDC)方法不受時間和空間頻散的影響，得到同相軸連續性更好、噪聲少、分辨率更高的像.圖11給出不同方法成像結果的波數譜.通過時間頻散校正和偽譜法壓制時間和空間頻散，可有效恢復出缺失的波數成分并削弱成像噪聲(如虛線框和箭頭所示).但實際應用中速度模型和震源子波未知，成像結果同時受速度和震源精度以及數值頻散的影響，導致不同方法的偏移效果差異不明顯.只有消除其他不利因素(如速度誤差、震源誤差、照明不均等)影響，本文發展的抗頻散逆時偏移方法才更加有效.

圖8 實際地震記錄Fig.8 Real seismic record

圖9 實際資料偏移速度模型Fig.9 Migration velocity model for real data

圖10 實際資料不同方法的成像結果(a) FD; (b) PS; (c FD+TDC; (d) PS+TDC.Fig.10 Images by different methods for real data

圖11 實際資料不同方法的成像結果波數譜(a) FD; (b) PS; (c) FD+TDC; (d) PS+TDC.Fig.11 Spectra of images by different methods for real data

4.2 全波形反演

Sigsbee2A模型如圖12所示.計算區域為5600 m×2960 m，網格大小為16 m×16 m，時間步長為1.5 ms，最大記錄長度為3.0 s.震源采用主頻為30 Hz的Ricker子波.36個震源和351個檢波點均勻分布于地表.初始速度模型為真實模型的平滑結果，如圖13所示.由于初始模型與真實模型相近，我們進行單尺度反演，不同方法迭代10次.采用偽譜法和小時間步長(Δt=0.375 ms)模擬無頻散觀測記錄.抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法中先在觀測記錄中加入時間頻散再計算殘差并進行波場反向延拓.圖14為第18炮(4000 m,2240 m)處不同方法第一次迭代時的正反向傳播波場和頻譜.FD表示時間2階空間12階精度的有限差分.采用偽譜法和小時間步長(Δt=0.375 ms)進行波場正反向延拓并反演作為參考解(Ref).由圖14可見：常規(FD)、抗空間頻散(PS)、抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法的正傳波場相對于參考解向時間負方向發生相移；常規(FD)和抗空間頻散(PS)的反傳波場相對于參考解向時間正方向發生相移，而抗時間頻散(FD+TDC)和抗時間+空間頻散PS+TDC)方法的反傳波場仍向時間負方向移動；正反向波場的空間頻散分別出現在正常波前時間增大和減小的方向，偽譜法(PS)可以有效壓制空間頻散.圖15給出不同反演方法第一次迭代時的正反傳波場匹配圖.時間頻散校正可以保證正反向波場中包含等量的(源-反射點路徑)時間頻散.抗時間頻散(FD+TDC)和抗時間+空間頻散(PS+TDC)方法的正反傳波場的時移量與參考解正反傳波場的時移量大致相同.正反傳波場的不匹配主要是由于速度存在誤差導致的.

圖12 Sigsbee2A模型Fig.12 Sigsbee2A model

圖13 Sigsbee2A模型初始模型Fig.13 Initial model for the Sigsbee2A model

圖14 Sigsbee2A模型(4000 m, 2240 m)處不同方法的正反向波場及頻譜(a) 正向波場； (b) 反向波場； (c) 正向波場頻譜； (d) 反向波場頻譜.Fig.14 Forward and backward wavefields and frequency spectra at the point (4000 m, 2240 m) by different methods for the Sigsbee2A model(a) Forward wavefield； (b) Backward wavefield； (c) Spectrum of (a)； (d)Spectrum of (b).

圖15 Sigsbee2A模型不同方法的正反向波場匹配圖(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.15 Matching diagrams of forward- and backward-propagated wavefields by different methods for the Sigsbee2A model

圖16為不同反演方法得到的速度模型.由圖16可知，常規(FD)和抗時間頻散(FD+TDC)方法受空間頻散影響嚴重，無法得到滿意的結果；抗空間頻散(PS)和抗時間+空間頻散(PS+TDC)方法反演的速度模型有一定改善.為進一步比較，圖17給出Distance = 2720 m處不同方法的反演結果.抗空間頻散(PS)方法受時間頻散影響，反演精度較低；抗時間+空間頻散(PS+TDC)方法可以緩解時間和空間頻散對反演的影響，獲得了與參考解更接近的結果.圖18展示了不同方法的收斂曲線.強的空間頻散導致常規(FD)和抗時間頻散(FD+TDC)方法無法收斂；抗時間+空間頻散(PS+TDC)方法受時間頻散影響較小，比抗空間頻散(PS)方法的收斂速度快.抗時間+空間頻散(PS+TDC)方法中震源和伴隨波場中殘留有源-反射點路徑的時間頻散，模擬和觀測波場也包含源-反射點-檢波點路徑的時間頻散.因此，抗頻散(PS+TDC)全波形反演方法仍無法徹底消除時間頻散的影響，其反演結果和收斂性與參考解仍有一定距離.

圖16 Sigsbee2A模型不同方法的反演結果(a) Ref; (b) FD; (c) PS; (d) FD+TDC; (e) PS+TDC.Fig.16 Inversion results by different methods for the Sigsbee2A model

圖17 Sigsbee2A模型Distance=2720 m處不同方法的反演結果Fig.17 Inversion results at distance 2720 m by different methods for the Sigsbee2A model

圖18 Sigsbee2A模型不同方法的收斂曲線Fig.18 Convergence curves by different methods for the Sigsbee2A model

下面采用一維初始模型(如圖19所示)對不同的反演方法進行測試.由于初始模擬與真實模型相差較遠，我們進行多尺度反演：0～5 Hz、0～10 Hz和0～15 Hz，每個尺度迭代10次.圖20給出Distance=2720 m處不同方法的反演結果.在大尺度低頻反演中，時間和空間頻散相對較弱，對全波形反演的影響非常小；通過時間頻散校正(TDC)減小時間頻散的影響，可以微弱改善反演精度.

圖19 Sigsbee2A模型一維初始模型Fig.19 1D initial model for the Sigsbee2A model

圖20 一維初始模型時Sigsbee2A模型Distance=2720 m處不同方法的反演結果Fig.20 Inversion results at distance 2720 m by different methods with the 1D initial model for the Sigsbee2A model

5 結論

本文詳細研究了有限差分的時間和空間頻散特性及其對逆時偏移和全波形反演的影響.通過時間有限差分+偽譜法和時間頻散校正+空間有限差分分離出時間和空間頻散，進而發展了抗時間頻散、抗空間頻散、抗時間+空間頻散的逆時偏移和全波形反演方法.模型算例和實際資料應用表明：逆時偏移同時受時間和空間頻散影響，時間頻散導致同相軸連續性變差、成像位置偏離，空間頻散產生高頻噪聲和虛假反射界面；全波形反演在低頻大尺度反演中幾乎不受時間和空間頻散影響，高頻精細反演中時間頻散會降低反演精度，空間頻散可能導致不收斂；抗頻散方法可以有效緩解時間和空間頻散影響，獲得滿意的偏移剖面和反演結果.

時間頻散校正和偽譜法不依賴于差分網格和波動方程，本文的抗頻散方法可以直接推廣到復雜介質(如變密度聲波、彈性、各向異性等)逆時偏移和全波形反演中.當模型復雜且存在速度劇烈變化界面時，偽譜法的模擬精度變差.通過對速度進行平滑或者采用譜元法可以緩解上述問題.時間頻散校正+優化空間有限差分也是一種兼顧精度和計算效率的方法.為徹底消除時間頻散對全波形反演的影響，需要進一步研究基于高階時間有限差分的抗頻散方法.

Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elastic wave equation.ChineseJournalofGeophysics(in Chinese), 43(3): 411-419.

Liu X, Liu Y, Ren Z M, et al. 2019. Perfectly matched layer boundary conditions for frequency-domain acoustic wave simulation in the mesh-free discretization.GeophysicalProspecting, 67(7): 1732-1744.

Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation.GeophysicalJournalInternational, 197(2): 1250-1267.