郭振波 孫鵬遠(yuǎn) 李培明 任曉喬 錢忠平 唐博文

(東方地球物理公司物探技術(shù)研究中心，河北涿州 072751)

0 引言

對于油氣地震勘探，需建立較精確的近地表模型以消除地表起伏及近地表速度異常對反射信號(hào)的影響，特別是對于山地等復(fù)雜探區(qū)，近地表模型的精度直接決定了最終的成像效果[1-2]。由于初至旅行時(shí)的穩(wěn)健性及旅行時(shí)層析的高效性，目前利用初至旅行時(shí)層析反演進(jìn)行近地表建模是應(yīng)用最廣泛的方法。近地表模型通常用于靜校正[3]、深度域建模[4]等處理流程。

經(jīng)典的初至層析反演方法基于射線追蹤方程，在反演中需要顯式計(jì)算炮點(diǎn)到檢波點(diǎn)的旅行時(shí)及射線路徑[5]，本文將這類經(jīng)典的初至層析反演方法稱為基于射線追蹤方程的初至層析反演方法。反演中通常采用兩種方式求取旅行時(shí)及射線路徑：①采用兩點(diǎn)射線追蹤[6]或類似方法同時(shí)求取旅行時(shí)及射線路徑；②采用波前追蹤類方法，如基于程函方程的有限差分類方法[7-8]、基于線性旅行時(shí)假設(shè)的線性旅行時(shí)插值方法[9]、基于圖理論的最短路徑方法[10]等，該類方法需在求解得到初至?xí)r間場的基礎(chǔ)上從檢波點(diǎn)開始反向追蹤得到射線路徑[9，11-12]。反演通常采用反向投影[13]、最小二乘奇異值分解(LSQR)[11]、同時(shí)迭代重構(gòu)(SIRT)[14]等方法?；谏渚€追蹤方程的層析反演方法已在實(shí)際生產(chǎn)中廣泛應(yīng)用，目前主要研究方向是如何與波形反演等精度較高的反演方法相結(jié)合[15-16]。

除了上述經(jīng)典算法外，還存在一類層析反演算法，本文將其稱為基于程函方程的層析反演方法。該類方法僅需要計(jì)算初至?xí)r間場而不需要計(jì)算射線路徑，借助伴隨狀態(tài)方法[17]求取模型更新量。該類方法最早由Sei等[18]于1994年提出，近年來得益于波形反演的興起，伴隨狀態(tài)方法為地球物理學(xué)家所熟知，進(jìn)而開展了廣泛的研究。Leung等[19]、Taillandier等[20]、謝春等[21]采用快速掃描方法進(jìn)行初至旅行時(shí)層析； Huang等[22]利用該方法進(jìn)行透射及反射旅行時(shí)層析； Benaichouche等[23]、李勇德等[24]討論了該類層析方法中預(yù)處理方法對層析的影響； Waheed等[25]將該類方法用于各向異性介質(zhì)的初至旅行時(shí)層析。

雖然兩類方法已得到了廣泛研究，但目前還沒有從理論及數(shù)值測試上進(jìn)行系統(tǒng)的對比分析。本文首先基于統(tǒng)一的反演框架從方法原理上對兩種方法進(jìn)行了詳細(xì)的對比分析，并做了定性的評價(jià)；然后通過一個(gè)數(shù)值測試實(shí)例，對兩種方法的反演精度、計(jì)算效率、內(nèi)存占用等情況進(jìn)行了定量的對比分析。

1 方法原理及對比分析

1.1 層析反演基本理論

旅行時(shí)層析通過迭代反演尋找正演旅行時(shí)與觀測旅行時(shí)最優(yōu)擬合的模型，若不考慮正則化項(xiàng)，L2模的目標(biāo)函數(shù)可寫為

(1)

將目標(biāo)函數(shù)在m0處展開且只保留至二階項(xiàng)，可整理為

(2)

通常將目標(biāo)函數(shù)相對于介質(zhì)參數(shù)的一階導(dǎo)數(shù)稱為梯度，即

(3)

將目標(biāo)函數(shù)相對于介質(zhì)參數(shù)的二階導(dǎo)數(shù)稱為Hessian矩陣，即

(4)

由于求取目標(biāo)函數(shù)最小，需要求取目標(biāo)函數(shù)的駐點(diǎn)，即g=0的點(diǎn)。對式(2)求導(dǎo)并將其設(shè)為零，可得

H|m0δm=-g|m0

(5)

δm=-(H|m0)-1g|m0

(6)

式(6)為在當(dāng)前模型m0下對應(yīng)模型更新量的顯式表示。利用δm更新當(dāng)前模型參數(shù)m0并將其作為下次迭代的初始模型，可實(shí)現(xiàn)非線性反演。當(dāng)目標(biāo)函數(shù)小于給定的閾值時(shí)，停止迭代，得到最終的反演模型。

由于Hessian矩陣的逆很難直接求取，通常通過迭代反演的方式求取式(5)所示的線性方程，此時(shí)反演方法通常稱為高斯—牛頓法。對于尺度較大的問題，Hessian矩陣H的求取非常耗時(shí)且存儲(chǔ)困難，通常需要對其進(jìn)行近似處理。若將Hessian矩陣近似為一單位矩陣，則退化為最速梯度方法；若將Hessian矩陣近似為對角或帶限矩陣，則退化為預(yù)處理的最速梯度方法。

1.2 基于射線追蹤方程的初至層析反演

對基于射線追蹤方程的初至層析反演，需要求取每個(gè)炮檢對的射線路徑，在得到射線路徑之后，可將初至?xí)r間表示為

(7)

(8)

假設(shè)每次迭代的模型更新量不大，不足以引起射線路徑的劇烈變化，即?Ls(m)/?m=0。利用該假設(shè)，將式(8)代入式(1)，利用式(3)可求得一階梯度為

(9)

利用式(4)可求得二階Hessian矩陣為

(10)

若利用對角元素近似Hessian矩陣，并將其代入式(6)，通過整理可得到反向投影算法的基本公式

(11)

式中δts，r為第s炮的第r個(gè)檢波點(diǎn)處的初至旅行時(shí)殘差。

將式(9)和式(10)代入式(5)，整理可得

(12)

1.3 基于程函方程的初至層析反演

基于程函方程的初至層析反演方法直接采用程函方程而不是射線追蹤方法作為反演中的正演方法。各向同性介質(zhì)中初至旅行時(shí)的傳播滿足程函方程，利用迎風(fēng)有限差分進(jìn)行離散之后的形式為

(Dxts)(Dxts)+(Dyts)(Dyts)+(Dzts)(Dzts)

=diag(m)m

(13)

滿足初始條件

ts(is)=0

(14)

式中：ts為第s炮對應(yīng)的所有網(wǎng)格點(diǎn)上的初至?xí)r間；Dx、Dy、Dz分別為x、y、z方向的迎風(fēng)差分矩陣；diag(m)為將慢度向量的元素放在對角位置形成的對角矩陣；is為第s炮震源網(wǎng)格序號(hào)。本文利用基于快速匹配方法的迎風(fēng)差分方法[26]求解式(13)。

(15)

將式(15)代入式(1)，利用式(3)可求得一階梯度為

(16)

diag(Dzts)Dz]-1diag(m)

(17)

將式(17)代入式(16),可得一階梯度為

(18)

由于很難直接求解上式中矩陣的逆，對于單炮梯度可通過有限差分方法求解如下方程得到

(19)

式(19)是Taillandier等[20]的連續(xù)形式偏微分方程式(13)的離散形式。本文對式(19)采用迎風(fēng)差分方法求解。

利用式(4)可求取Hessian矩陣為

(20)

由于式(17)中存在矩陣的逆，通常無法直接求解。將式(20)代入式(5)，通過整理可獲取其等價(jià)的形式為

(21)

若在每次非線性迭代內(nèi)部迭代求解式(21)即為高斯—牛頓法。若只保留Hessian矩陣對角元素，則為預(yù)處理的最速下降方法。

1.4 對比分析

上文基于相同的目標(biāo)函數(shù)，在統(tǒng)一的反演框架下分別推導(dǎo)了基于射線追蹤方程與程函方程初至旅行時(shí)層析迭代反演的基礎(chǔ)公式。詳細(xì)對比兩者的理論推導(dǎo)過程，兩類方法具有如下共性。

(1)具有相同的理論基礎(chǔ)。理論初至旅行時(shí)的計(jì)算均采用基于高頻近似的射線理論進(jìn)行；初至旅行時(shí)反演采用統(tǒng)一的目標(biāo)函數(shù)，且均可在統(tǒng)一的反演框架下推導(dǎo)迭代反演的基本公式。

(2)具體反演算法之間具有對應(yīng)性?；谏渚€追蹤方程方法中的反向投影算法與基于程函方程方法的預(yù)處理最速下降方法是對應(yīng)的；基于射線追蹤方程方法中的LSQR、SIRT等與基于程函方程方法中的高斯—牛頓算法是對應(yīng)的。對應(yīng)方法之間可產(chǎn)生相近的效果。

由于兩者在理論上的相似性，兩者具有相似的分辨率以及處理復(fù)雜問題的能力，如兩者均無法解決速度反轉(zhuǎn)等復(fù)雜探區(qū)的問題。

兩種方法在理論上具有如下的主要區(qū)別。

(1)二者所基于的正演算子不同。前者基于射線追蹤方程，即旅行時(shí)為沿射線路徑對慢度的積分；后者基于程函方程，即利用有限差分等數(shù)值解法通過求解程函方程得到檢波點(diǎn)處的初至?xí)r間。

(2)二者求取或構(gòu)造目標(biāo)函數(shù)相對于介質(zhì)參數(shù)一階梯度及二階Hessian矩陣的方式不同。對于前者，可利用不同網(wǎng)格內(nèi)射線段長度顯式地構(gòu)建一階及二階導(dǎo)數(shù)所對應(yīng)的矩陣(式(9)、式(10))，直觀且物理意義明確；對于后者，只能通過隱式數(shù)值求解偏微分方程(式(19))得到單炮所對應(yīng)的梯度，對于Hessian矩陣則很難顯式地進(jìn)行構(gòu)建，只能通過間接的方式考慮Hessian矩陣的影響，物理意義相對不明確。

基于理論上的異同，結(jié)合具體實(shí)現(xiàn)細(xì)節(jié)(如圖1所示兩種方法的反演流程對比)，在實(shí)際應(yīng)用中，基于射線追蹤方程的初至層析反演方法具有如下優(yōu)勢。

(1)反演過程中針對性處理更為靈活?？衫靡阎纳渚€路徑計(jì)算射線角度、射線密度等信息，可利用這些信息在反演的過程中進(jìn)行針對性處理。

(2)豐富的質(zhì)控方法。將射線路徑、射線角度、射線密度等利用圖形化的方式進(jìn)行顯示，可對最終的反演方法進(jìn)行有效的質(zhì)控，如可檢查射線是否到達(dá)模型邊界、根據(jù)射線密度判斷不同部位反演模型的可靠程度等。

基于程函方程的初至層析反演方法具有如下優(yōu)勢。

(1)避免了射線路徑的計(jì)算以及由于射線路徑計(jì)算引入的誤差。由于初至旅行時(shí)只考慮最小時(shí)間且由于地表起伏可能引起的散射，實(shí)際處理中很少直接進(jìn)行顯式的射線追蹤，而是在利用程函方程計(jì)算得到初至?xí)r間場的基礎(chǔ)上通過從檢波點(diǎn)處反向追蹤得到射線路徑。不管是采用梯度方法[7]還是線性旅行插值方法[9]，均基于局部近似，求取射線路徑存在一定誤差，特別是在地表起伏等復(fù)雜近地表?xiàng)l件下。若采用基于射線追蹤方程的層析反演，需顯式進(jìn)行射線路徑的求取，在增加計(jì)算量的同時(shí)也在反演過程中也引入了部分誤差。由于不需要進(jìn)行射線路徑的求取，因此基于程函方程的層析反演有效地避免了這個(gè)問題。

圖1 兩種層析反演方法反演流程對比

圖2 兩種反演方法核函數(shù)對比背景圖片為基于程函方程層析反演的核函數(shù)； 紅色

(2)計(jì)算效率及占用內(nèi)存大小主要受模型大小的影響，與檢波點(diǎn)個(gè)數(shù)基本無關(guān)。在基于程函方程層析反演內(nèi)部，主要利用有限差分求解式(13)、式(19)和式(21)，這些方程的計(jì)算效率及所需內(nèi)存的大小僅與所用模型的大小有關(guān)，與檢波點(diǎn)的個(gè)數(shù)無關(guān)。

(3)如圖2所示，基于程函方程的層析反演方法核函數(shù)具有一定的寬度，在復(fù)雜區(qū)域反演效果略好。核函數(shù)代表空間任一位置一個(gè)散射點(diǎn)對走時(shí)變化的影響[27-28]，決定了梯度的形態(tài)。主要受差分算子的影響，基于程函方程的層析反演其核函數(shù)具有一定的寬度，類似于有限頻效應(yīng)，在復(fù)雜區(qū)域理論上具有更好的反演效果。

2 數(shù)值對比分析

為了對比兩種方法在反演精度、計(jì)算效率、內(nèi)存占用等方面的差異，選用Amoco靜校正基準(zhǔn)測試模型(圖3)進(jìn)行數(shù)值測試。該模型包含了大部分常見的近地表地質(zhì)構(gòu)造，如高速層出露、局部高速、低速異常體、淺層低速層、近地表復(fù)雜構(gòu)造及極淺層低速體等，可進(jìn)行較全面的對比分析。由于常規(guī)近地表建模通常的反演深度大約在1km左右，圖3所示模型是對原有模型在深度方向進(jìn)行壓縮之后的結(jié)果，在保留原有模型復(fù)雜性的同時(shí)使其深度與常規(guī)近地表問題一致。

圖3 Amoco靜校正基準(zhǔn)測試模型

2.1 反演精度對比

兩種反演方法內(nèi)部的正演算子均通過數(shù)值求解程函方程實(shí)現(xiàn)，具體采用基于快速匹配追蹤的迎風(fēng)差分算法[26]。對于基于射線追蹤方程的層析反演，在通過求解程函方程獲取時(shí)間場之后由檢波點(diǎn)位置開始沿初至?xí)r間場的負(fù)梯度方向追蹤射線路徑，最后計(jì)算每條射線所經(jīng)過網(wǎng)格內(nèi)射線段的長度來構(gòu)建式(7)中對應(yīng)的正演算子[13]?；诔毯匠痰膶游龇囱莶捎妙A(yù)處理的最速梯度方法，預(yù)處理算子采用與Benaichouche等[23]、李勇德等[24]相同的形式，式(19)對應(yīng)的方程采用迎風(fēng)有限差分方法進(jìn)行求解；基于射線追蹤方程的層析反演方法采用與最速梯度法對應(yīng)的反向投影算法。除以上差別外，兩種反演方法均采用相同的處理方式，數(shù)值試驗(yàn)時(shí)也選用相同的參數(shù)。

以25m為間隔在地表均勻布置共1998炮，每炮最大炮檢距為7500m，檢波點(diǎn)間距為20m。采用數(shù)值正演方法為層析反演生成合成數(shù)據(jù)。選用圖4所示線性遞增模型為初始模型，地表速度為4500m/s，速度隨深度變化梯度為0.5，在模型底部速度約為5500 m/s。圖5為基于射線追蹤方程的層析反演結(jié)果，圖6為基于程函方程的層析反演結(jié)果。對比圖5與圖6可知，從整體形態(tài)上來說，兩種反演結(jié)果基本一致，高速層出露、局部高速、低速異常體、淺層低速層、極淺層低速體等近地表地質(zhì)構(gòu)造均得到較好的反演。由于射線層析分辨率的限制，模型右側(cè)復(fù)雜構(gòu)造區(qū)域均未取得理想的效果。從反演細(xì)節(jié)上來說，由于基于程函方程的層析反演其核函數(shù)是帶限的，即具有類似有限頻的特性，反演結(jié)果更加穩(wěn)定，特別是在復(fù)雜構(gòu)造區(qū)域，如圖5黑色箭頭所示區(qū)域，基于射線追蹤方程的層析反演結(jié)果存在局部的高速異常，但在基于程函方程的層析反演結(jié)果中不存在這類高速異常。從反演精度上來說，兩類方法具有相似的反演精度，在局部細(xì)節(jié)上基于程函方程的層析反演方法略優(yōu)。由圖7的兩種方法殘差曲線對比可知，兩者具有相似的收斂效率。

圖4 初始模型

圖5 基于射線追蹤方程的層析反演結(jié)果

圖6 基于程函方程的層析反演結(jié)果

圖8為基于射線追蹤方程層析反演最后一次迭代的射線密度，利用射線密度可分析不同區(qū)域速度反演的可靠性。由于基于程函方程的層析反演沒有射線密度的概念，因此缺少類似的質(zhì)控圖件。

2.2 內(nèi)存占用及計(jì)算效率對比

圖7 兩種方法反演殘差曲線對比

圖8 基于射線追蹤方程層析反演最后迭代的射線密度

為了在內(nèi)存占用、計(jì)算效率等方面對兩種方法進(jìn)行更全面的對比，設(shè)計(jì)了四種不同的觀測系統(tǒng)進(jìn)行測試。由前面的理論分析可發(fā)現(xiàn)，基于射線追蹤方程層析反演的一個(gè)主要特點(diǎn)是其內(nèi)存占用、計(jì)算效率等均是與檢波點(diǎn)個(gè)數(shù)相關(guān)的?；谶@個(gè)主要差別，在固定總炮數(shù)、最大炮檢距的基礎(chǔ)上通過設(shè)置不同的檢波點(diǎn)間距(5、10、20、40m)模擬不同的觀測系統(tǒng)。測試環(huán)境為Linux集群，操作系統(tǒng)為企業(yè)版Red Hat 6.4，CPU為Intel(R) Xeon(R) CPU E5-2670(2.60GHz)，編譯器為Intel icpc。程序源代碼采用C++編寫，內(nèi)部采用MPI并行處理。數(shù)值測試共用7個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)有30個(gè)線程，每次層析反演迭代30次。

在不同觀測系統(tǒng)下內(nèi)存占用(單進(jìn)程對應(yīng)內(nèi)存占用)及運(yùn)行時(shí)間對比如圖9所示，可以看出：①基于程函方程的層析反演內(nèi)存占用及計(jì)算時(shí)間與檢波點(diǎn)個(gè)數(shù)無關(guān)，只與當(dāng)前模型的大小有關(guān)；②基于射線追蹤方程的層析反演占用內(nèi)存及計(jì)算時(shí)間均是檢波點(diǎn)個(gè)數(shù)的函數(shù)，檢波點(diǎn)個(gè)數(shù)越多，占用內(nèi)存越大、計(jì)算時(shí)間越長；③當(dāng)檢波點(diǎn)較密時(shí)，如寬方位、高密度等觀測數(shù)據(jù)，基于程函方程的層析反演方法在內(nèi)存占用、計(jì)算效率等具有較大優(yōu)勢；當(dāng)檢波點(diǎn)稀疏時(shí)，基于射線追蹤方程的層析反演方法具有更大的優(yōu)勢。

圖9 兩種方法在不同檢波點(diǎn)間距下內(nèi)存占用(a)和運(yùn)行時(shí)間(b)對比

3 結(jié)論

本文在統(tǒng)一的反演框架下推導(dǎo)了兩種層析方法的基本公式，在理論上對兩種方法的主要差別進(jìn)行了定性的分析；通過理論數(shù)值測試，對兩種方法在反演精度、計(jì)算效率等方面進(jìn)行了定量的分析?；谝陨系姆治觯梢缘贸鋈缦陆Y(jié)論：

(1)兩種方法在理論上具有對應(yīng)性，可在統(tǒng)一的反演框架下導(dǎo)出，且不同迭代反演方法之間均具有對應(yīng)關(guān)系，主要的差別是由于反演中正演算子的不同引起；

(2)當(dāng)基于射線追蹤方程的方法采用反向投影反演方法、基于程函方程的方法采用與之對應(yīng)的最速梯度方法反演方法時(shí)，兩者具有相近的反演精度，基于程函方程的反演方法核函數(shù)具有類似有限頻的特性，在復(fù)雜構(gòu)造中的反演更加穩(wěn)定；

(3)計(jì)算效率及內(nèi)存占用是兩者的主要差別。對于基于射線追蹤方程的層析反演方法，計(jì)算效率等依賴于檢波點(diǎn)個(gè)數(shù)，在檢波點(diǎn)特別稀疏時(shí)具有優(yōu)勢；對于基于程函方程的層析反演方法，其主要依賴于模型大小，在檢波點(diǎn)間距較小時(shí)具有優(yōu)勢；

(4)基于射線追蹤方程的層析反演方法具有射線密度等質(zhì)控手段，是一個(gè)主要的優(yōu)勢。

鑒于目前寬方位、高密度數(shù)據(jù)的大量采集，建議在實(shí)際中以應(yīng)用基于程函方程的層析反演方法為主；在層析反演過程中選用部分炮進(jìn)行射線追蹤以指導(dǎo)、控制反演過程并生成射線密度等質(zhì)控圖件，從而充分利用不同方法之間的優(yōu)勢。