改進(jìn)的伴隨狀態(tài)法初至波走時(shí)層析成像方法

董良國(guó)， 張建明， 韓佩恩

同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

0 引言

反演地下介質(zhì)速度分布是地震學(xué)的核心問(wèn)題之一.鑒于初至波最先到達(dá)、能量強(qiáng)、信息可靠、易于拾取等優(yōu)點(diǎn)，目前被大量研究和應(yīng)用，并產(chǎn)生了利用初至波不同信息的多種反演方法，如初至波走時(shí)反演(Olson, 1989)、初至波包絡(luò)反演(敖瑞德等，2015)、初至波相位反演(Choi and Alkhalifah, 2013)、初至波波形反演(Sheng et al., 2006)以及初至波多信息聯(lián)合反演(Liu and Zhang, 2017; 張建明等，2020)等方法.

從理論上講，利用初至波動(dòng)力學(xué)信息可以取得更高的反演精度，但影響初至波動(dòng)力學(xué)信息的因素比較多，信息的穩(wěn)定性差，導(dǎo)致反演的非線性強(qiáng)，對(duì)初始模型要求更高，也更容易陷入局部極值.而初至波走時(shí)信息基本上只受速度的宏觀分布影響，信息可靠性高，因此初至波走時(shí)反演方法穩(wěn)定性更強(qiáng)，在近地表以及井間速度建模等方面得到了更廣泛的應(yīng)用.

使用不同的地震波正演引擎，產(chǎn)生了不同的初至波走時(shí)反演方法，例如，基于射線追蹤(Bishop et al., 1985)、基于程函方程走時(shí)計(jì)算(Taillandier et al., 2009)以及基于波動(dòng)方程模擬(Luo and Schuster, 1991)的初至波走時(shí)反演方法.在這些方法中，不需要射線追蹤、基于伴隨狀態(tài)法的初至波走時(shí)層析成像方法(Adjoint State Tomography method，以下簡(jiǎn)稱(chēng)AST方法)近年來(lái)得到了很大的關(guān)注.這種方法基于伴隨狀態(tài)法、通過(guò)程函方程計(jì)算的走時(shí)場(chǎng)和伴隨方程計(jì)算的伴隨場(chǎng)來(lái)計(jì)算目標(biāo)函數(shù)的梯度，避免了射線追蹤和Fréchet導(dǎo)數(shù)矩陣的計(jì)算.該方法對(duì)內(nèi)存的要求主要依賴(lài)于地下模型的網(wǎng)格數(shù)，與觀測(cè)的地震數(shù)據(jù)量無(wú)關(guān)，內(nèi)存需求極低，特別適合于并行計(jì)算(Noble et al.，2010；Benaichouche et al.，2015).至于伴隨狀態(tài)法，自從它首次在優(yōu)化控制領(lǐng)域被提出以來(lái)(Lions, 1971)，目前已經(jīng)被廣泛應(yīng)用于包括地震學(xué)在內(nèi)的許多領(lǐng)域中(Plessix，2006；Fichtner and Trampert, 2011).

基于程函方程的AST方法最早由Sei和Symes(1994)所提出，自此之后，幾乎所有的關(guān)于AST方法的理論方法及應(yīng)用的文獻(xiàn)，都是基于面積分來(lái)定義目標(biāo)函數(shù)，由此得到的伴隨方程都依賴(lài)于地表的法向量，本文稱(chēng)這種方法為傳統(tǒng)AST方法.

隨后，這種基于面積分來(lái)定義目標(biāo)函數(shù)的AST方法得到了不斷的發(fā)展和完善.Leung和Qian(2006)重新推導(dǎo)得到了與觀測(cè)系統(tǒng)地表法向量有關(guān)的伴隨狀態(tài)法走時(shí)層析成像方法，構(gòu)造了伴隨方程的差分格式.Taillandier等(2009)以管量場(chǎng)概念對(duì)反傳播方程的計(jì)算格式和意義作了進(jìn)一步說(shuō)明，但從其中的地表形態(tài)不對(duì)稱(chēng)的理論模型試驗(yàn)結(jié)果中，可明顯看到伴隨場(chǎng)的不對(duì)稱(chēng)現(xiàn)象，進(jìn)而導(dǎo)致了梯度形態(tài)畸變的不合理現(xiàn)象.Noble等(2010)指出了伴隨狀態(tài)法層析成像技術(shù)在并行計(jì)算時(shí)所具有的巨大效率優(yōu)勢(shì).Bretaudeau等(2014)進(jìn)一步發(fā)展了基于高斯-牛頓優(yōu)化方法的二階伴隨狀態(tài)法走時(shí)層析方程.為了在走時(shí)反演中有效刻畫(huà)地下模型的界面信息，在AST方法中引入了水平集技術(shù)(Li and Leung, 2013; Li et al., 2014)，可以考慮多傳播路徑問(wèn)題.這種依賴(lài)于地表邊界法向量的傳統(tǒng)AST方法，目前已經(jīng)發(fā)展到透射和反射走時(shí)的同時(shí)反演(Huang and Bellefleur，2012；Li et al., 2014)以及各向異性介質(zhì)模型參數(shù)的反演(Waheed et al., 2014,2016)中.

在國(guó)內(nèi)，謝春等(2014)討論了基于射線理論的伴隨狀態(tài)法初至波走時(shí)層析成像技術(shù)，并擴(kuò)展至基于有限頻的AST走時(shí)層析成像(謝春等, 2015).為優(yōu)化反演效率和提高反演效果，李勇德等(2017)提出了使用近似Hessian矩陣進(jìn)行預(yù)條件的AST走時(shí)反演層析成像方法.

我們已經(jīng)發(fā)現(xiàn)，這種傳統(tǒng)AST方法無(wú)論在理論方法上還是在應(yīng)用過(guò)程中還存在一定的問(wèn)題(Han et al., 2019)，本文的目的是對(duì)這些問(wèn)題進(jìn)行進(jìn)一步的討論.首先對(duì)這種傳統(tǒng)AST方法簡(jiǎn)要介紹之后，從理論模型試驗(yàn)以及地震波傳播物理規(guī)律的角度指出了這些問(wèn)題的存在，并提出了改進(jìn)的方法，即不依賴(lài)于地表法向量的改進(jìn)的AST走時(shí)反演層析成像方法，并通過(guò)兩個(gè)模型反演試驗(yàn)證明了改進(jìn)方法的正確性和有效性.至于傳統(tǒng)AST方法存在的問(wèn)題的嚴(yán)格數(shù)學(xué)證明，是我們下一步的工作目標(biāo).

需要指出的是，本文采用快速掃描法(Zhao, 2005)求解程函方程和伴隨方程，確定走時(shí)場(chǎng)和伴隨場(chǎng)，從而計(jì)算AST方法的迭代梯度.

1 目前傳統(tǒng)AST方法存在的問(wèn)題

1.1 傳統(tǒng)AST方法簡(jiǎn)介

觀測(cè)和理論走時(shí)達(dá)到最佳匹配是走時(shí)反演的基本原則.目前傳統(tǒng)伴隨狀態(tài)法初至波走時(shí)層析普遍采用面積分來(lái)定義目標(biāo)函數(shù)(Leung and Qian，2006；Taillandier et al.，2009；Noble et al.，2010；Huang and Bellefleur，2012; Li and Leung, 2013; Li et al., 2014；Bretaudeau et al.，2014；Waheed et al., 2014; 謝春等，2014，2015；Benaichouche et al.，2015；Waheed et al., 2016；李勇德等，2017)：

(1)

其中，r表示檢波點(diǎn)的位置，檢波點(diǎn)位于區(qū)域Ω的邊界?Ω上，Tobs(r)表示r處的實(shí)際觀測(cè)走時(shí)，c(x)表示空間位置x處的介質(zhì)速度.t(r)表示在當(dāng)前速度模型下空間位置r處計(jì)算的理論走時(shí)，在高頻假設(shè)下，走時(shí)場(chǎng)t(r)滿(mǎn)足程函方程：

(2)

利用伴隨狀態(tài)法，可以確定上述目標(biāo)函數(shù)(1)的梯度(Leung and Qian，2006；Taillandier et al.，2009；Noble et al.，2010；Li and Leung, 2013; Li et al., 2014；Waheed et al., 2014, 2016)：

(3)

其中，伴隨場(chǎng)λ(x)滿(mǎn)足：

(4)

(5)

(4)式中的n為邊界?Ω上位置r處的邊界外法向單位矢量.

利用(4)式，根據(jù)檢波點(diǎn)處的走時(shí)殘差設(shè)置檢波點(diǎn)處的伴隨場(chǎng)λ的初值后，再利用方程(5)將走時(shí)殘差反向傳播即可獲得伴隨場(chǎng)λ(x)，進(jìn)而通過(guò)公式(3)就可以得到目標(biāo)函數(shù)對(duì)速度的梯度，然后通過(guò)線性搜索算法可以很方便地獲得一個(gè)優(yōu)化的步長(zhǎng)αn(其中n為迭代次數(shù))，這樣就可以利用下面的局部?jī)?yōu)化算法進(jìn)行速度迭代，得到最終的速度反演結(jié)果.

(6)

可見(jiàn)，這種AST方法的核心是利用伴隨方程計(jì)算伴隨場(chǎng)λ(x)，關(guān)鍵是伴隨方程.自從Sei和Symes(1994)提出AST走時(shí)反演方法之后，有關(guān)AST方法的所有文獻(xiàn)中，所用的伴隨場(chǎng)方程的邊界條件均是(4)式，該式依賴(lài)于邊界的法向量.

1.2 傳統(tǒng)AST方法存在的問(wèn)題

下面具體分析上述傳統(tǒng)AST方法存在的問(wèn)題.

問(wèn)題一：傳統(tǒng)AST方法采用面積分來(lái)定義目標(biāo)函數(shù)，決定了其無(wú)法處理檢波器在模型內(nèi)部的情況.

方程(1)采用面積分來(lái)定義目標(biāo)函數(shù)，就預(yù)先假定了所有檢波器均位于模型表面，這就造成了傳統(tǒng)AST方法只適用于檢波點(diǎn)在邊界上、無(wú)法考慮檢波器在模型內(nèi)部的情況.當(dāng)進(jìn)行井中觀測(cè)以及井間觀測(cè)時(shí)，檢波器在模型內(nèi)部，這時(shí)方程(1)的面積分定義目標(biāo)函數(shù)的方式就不適用，伴隨方程(4)中的所謂的邊界法向當(dāng)然無(wú)從談起.可見(jiàn)，對(duì)于檢波點(diǎn)在地下的情況(如井中觀測(cè))，這種傳統(tǒng)的AST方法是不適用的.

問(wèn)題二：傳統(tǒng)AST方法中，伴隨場(chǎng)的計(jì)算依賴(lài)于模型表面的法向量，這是有悖射線理論實(shí)際的.

設(shè)計(jì)一個(gè)除地表形狀外中間對(duì)稱(chēng)的高速異常體模型(圖1)，初始速度也是兩邊對(duì)稱(chēng)的縱向梯度模型(圖2).在該模型的中心地表激發(fā)，在兩側(cè)具有相同橫向距離和同等深度、但地表形態(tài)不同的兩個(gè)位置(圖1中兩個(gè)R點(diǎn)處)進(jìn)行接收.由于炮點(diǎn)兩側(cè)真實(shí)和初始速度模型以及兩個(gè)觀測(cè)點(diǎn)都是完全對(duì)稱(chēng)的，只要地表起伏不是極端劇烈，從地震波傳播的客觀現(xiàn)實(shí)來(lái)講，第一輪迭代的梯度在激發(fā)點(diǎn)兩側(cè)也應(yīng)該是對(duì)稱(chēng)的.但是，在首次迭代中，盡管左右兩個(gè)檢波點(diǎn)具有相同的走時(shí)差t(r)-Tobs(r)，但由于左右兩個(gè)檢波點(diǎn)處的地表法向量n不同，因此，利用依賴(lài)于地表法向量的傳統(tǒng)AST方法得到的伴隨場(chǎng)兩側(cè)并不對(duì)稱(chēng)(圖3)，由此得到的梯度場(chǎng)當(dāng)然也不對(duì)稱(chēng)(圖4).從這個(gè)簡(jiǎn)單的例子可以看出，傳統(tǒng)AST方法的伴隨方程存在明顯的不合理之處.

圖1 不同地表法向量的對(duì)稱(chēng)接收點(diǎn)試驗(yàn)的真實(shí)模型Fig.1 The true velocity model with two symmetrical receivers but different surface normal vectors

圖2 初始速度模型Fig.2 The initial velocity model

圖3 傳統(tǒng)AST方法得到的歸一化后的伴隨場(chǎng)分布Fig.3 The normalized adjoint field distribution by the conventional AST method

圖4 傳統(tǒng)AST方法計(jì)算的梯度Fig.4 The gradient obtained by the conventional AST method

問(wèn)題三：當(dāng)?shù)乇砗途型瑫r(shí)觀測(cè)時(shí)，傳統(tǒng)AST方法無(wú)法正確計(jì)算梯度.

2 改進(jìn)的AST方法

若使用傳統(tǒng)AST方法中的方程(1)來(lái)定義走時(shí)層析目標(biāo)函數(shù)，檢波點(diǎn)只能位于區(qū)域Ω的邊界?Ω上.我們將?Ω′面定義為?Ω面向外的一個(gè)拓展面，拓展的空間部分充填空氣，速度為340 m·s-1.這樣，原來(lái)地表上以及井中的檢波點(diǎn)均變?yōu)閰^(qū)域Ω′的內(nèi)部點(diǎn)，如圖5.由于實(shí)際地表之上本來(lái)就是空氣層，并且在區(qū)域Ω的邊界?Ω為凸邊界的條件下，這樣的拓展并不會(huì)改變區(qū)域Ω內(nèi)的地震波走時(shí).

圖5 計(jì)算區(qū)域的拓展示意圖Fig.5 The extended sketch map of the computational domain

現(xiàn)在基于拓展的新空間Ω′及其邊界?Ω′，利用體積分定義伴隨狀態(tài)法初至波走時(shí)層析目標(biāo)函數(shù)：

(7)

將控制方程(2)與伴隨變量λ(即乘子)引入目標(biāo)函數(shù)(7)中，形成新的目標(biāo)函數(shù)：

(8)

在最優(yōu)化過(guò)程中，走時(shí)場(chǎng)變量t、伴隨變量λ和速度c相互獨(dú)立，所以當(dāng)目標(biāo)函數(shù)達(dá)到極小值時(shí)，目標(biāo)函數(shù)仍滿(mǎn)足以下三個(gè)關(guān)系：

(9)

(10)

(11)

由(9)式可以得到控制方程，即程函方程(2)式.而由(11)式可以得到目標(biāo)函數(shù)(7)的梯度：

(12)

可見(jiàn)，由于每一次迭代中的當(dāng)前速度模型c(x)已知，伴隨場(chǎng)λ(x)的求取是確定目標(biāo)函數(shù)梯度的關(guān)鍵.而由(10)式就可以導(dǎo)出控制伴隨場(chǎng)λ(x)的伴隨方程.

我們使用擾動(dòng)法求解(10)式，得到了改進(jìn)的AST方法的伴隨方程(推導(dǎo)過(guò)程見(jiàn)附錄A)：

x∈Ωand ?Ω.

(13)

由于每一次迭代中的理論走時(shí)場(chǎng)t(x)可以通過(guò)求解程函方程(2)式確定，這樣，求解(13)式的伴隨方程即可確定伴隨場(chǎng)λ(x)，再利用(12)式即可獲得目標(biāo)函數(shù)的梯度，進(jìn)而利用(6)式就可以進(jìn)行速度迭代反演.

需要指出的是，我們這種改進(jìn)的AST方法與傳統(tǒng)AST方法的唯一差別，就是確定伴隨場(chǎng)λ(x)的伴隨方程不同.傳統(tǒng)AST方法的伴隨方程是(4)和(5)式，利用(4)式設(shè)置邊界上檢波點(diǎn)處的伴隨場(chǎng)λ的初值，再利用(5)式確定伴隨場(chǎng)λ(x).而我們這種改進(jìn)的AST方法的伴隨方程是統(tǒng)一的(13)式，地表和地下檢波點(diǎn)具有相同的處理方式，并且伴隨場(chǎng)的分布與模型邊界的法向量無(wú)關(guān)，自然克服了前面指出的傳統(tǒng)AST方法存在的三個(gè)缺陷.

也就是說(shuō)，改進(jìn)的AST方法適用于各種地表和井中同時(shí)觀測(cè)方式，地下和井中觀測(cè)數(shù)據(jù)的處理方式得到了統(tǒng)一，從理論方法上避免了Waheed等(2016)處理井中觀測(cè)時(shí)的井中和地表觀測(cè)數(shù)據(jù)對(duì)梯度貢獻(xiàn)不一致的問(wèn)題.無(wú)論是內(nèi)部點(diǎn)還是邊界上，改進(jìn)方法的伴隨方程都遵從統(tǒng)一形式，伴隨場(chǎng)的計(jì)算與模型邊界的法向量無(wú)關(guān)，不需要區(qū)分檢波點(diǎn)位置處的地表法向量，這樣更符合實(shí)際物理意義，保證了梯度場(chǎng)的正確性.

仍然采用圖1所示的起伏地表高速異常體模型，初始模型仍然為圖2所示的縱向常梯度模型，使用上述改進(jìn)的AST方法計(jì)算得到的平滑后的首輪迭代梯度如圖6所示.可以發(fā)現(xiàn)，炮點(diǎn)兩側(cè)的梯度場(chǎng)是對(duì)稱(chēng)的，這是符合理論預(yù)期的，證明了我們得到的改進(jìn)AST方法中新的伴隨方程的正確性.

圖6 改進(jìn)的AST方法計(jì)算的歸一化后的梯度分布Fig.6 The normalized gradient obtained by the improved AST method

3 模型試驗(yàn)

下面，通過(guò)兩個(gè)理論模型試驗(yàn)來(lái)說(shuō)明改進(jìn)的AST方法的正確性及反演效果.

3.1 簡(jiǎn)單起伏地表模型

為此，設(shè)計(jì)一個(gè)地表起伏比較劇烈的理論模型(見(jiàn)圖7a)，在2000 m·s-1的均勻速度模型內(nèi)部存在一個(gè)速度為3000 m·s-1高速異常體.在地表等間距激發(fā)的201炮產(chǎn)生的地表處的初至波走時(shí)數(shù)據(jù)作為“觀測(cè)數(shù)據(jù)”，相鄰炮點(diǎn)的水平間距為20 m.每炮都在地表均勻設(shè)置401個(gè)檢波器，相鄰檢波點(diǎn)水平間距為10 m.

圖7 簡(jiǎn)單起伏地表速度模型及改進(jìn)前后AST反演結(jié)果(a) 簡(jiǎn)單起伏地表速度模型； (b) AST-IV反演結(jié)果； (c) AST-DV反演結(jié)果.Fig.7 The simple topographic velocity model and the AST inversion results(a) The simple topographic velocity model; (b) The AST-IV result; (c) The AST-DV result.

模型離散為401×101網(wǎng)格點(diǎn)，縱橫向網(wǎng)格間距均為10 m，采用速度值為2000 m·s-1的常速模型作為反演的初始模型.分別使用改進(jìn)的不依賴(lài)于地表法向量的AST方法(AST-IV)和依賴(lài)于地表法向量的傳統(tǒng)AST方法(AST-DV)進(jìn)行試驗(yàn)，反演結(jié)果分別見(jiàn)圖7b和圖7c.其中，AST-IV和AST-DV分別迭代了693代和575代，它們的目標(biāo)函數(shù)都不再下降.

可以看到，改進(jìn)方法比傳統(tǒng)方法能更精確地反演出近地表附近的局部高速異常，而傳統(tǒng)方法的反演結(jié)果非常差.由于二者的唯一差別，就是確定伴隨場(chǎng)λ(x)的伴隨方程不同，說(shuō)明傳統(tǒng)AST方法在伴隨方程的求解過(guò)程中，檢波點(diǎn)處的走時(shí)殘差反向傳播是不正確的，導(dǎo)致其計(jì)算的梯度和模型修改量不正確，說(shuō)明了傳統(tǒng)AST方法的伴隨場(chǎng)依賴(lài)于地表法向量是不合理的.

圖8展示了改進(jìn)AST方法與傳統(tǒng)AST方法不同迭代次數(shù)(第1、20、100、300次)的歸一化梯度值.可以發(fā)現(xiàn)，隨著迭代的進(jìn)行，改進(jìn)AST方法的梯度形態(tài)逐漸變得與高速異常體形態(tài)一致，而傳統(tǒng)AST方法的梯度沒(méi)有得到本質(zhì)的改善，反演過(guò)早就陷入了局部極值.這進(jìn)一步證實(shí)了傳統(tǒng)AST方法在伴隨方程的求解過(guò)程中，檢波點(diǎn)處的走時(shí)殘差反向傳播是不正確的，導(dǎo)致其梯度和模型修改量形態(tài)異常，造成該問(wèn)題的根源就是來(lái)自于傳統(tǒng)AST方法的伴隨場(chǎng)錯(cuò)誤地依賴(lài)于地表法向量，當(dāng)?shù)乇砥鸱鼊×視r(shí)，這個(gè)問(wèn)題就很突出.而改進(jìn)的AST方法的伴隨場(chǎng)與地表法向量無(wú)關(guān)，它可以有效修正梯度，從而最終能得到正確的反演結(jié)果.

圖8 傳統(tǒng)與改進(jìn)的AST方法不同迭代次數(shù)的歸一化梯度從上到下分別為第1、20、100、300代的梯度(左側(cè)為傳統(tǒng)的AST方法，右側(cè)為改進(jìn)的AST方法).Fig.8 The normalized gradients of different iterations for the conventional and the improved AST methodsFrom top to bottom show the gradients corresponding to 1, 20, 100 and 300 iterations. The AST-DV and the AST-IV results are shown on the left and right respectively.

從本試驗(yàn)的目標(biāo)函數(shù)收斂曲線(圖9)上可以發(fā)現(xiàn)，改進(jìn)AST方法不僅收斂速度明顯快于傳統(tǒng)AST方法，且反演結(jié)果的走時(shí)殘差遠(yuǎn)小于傳統(tǒng)AST方法的走時(shí)殘差.這主要是由于傳統(tǒng)AST方法的伴隨場(chǎng)錯(cuò)誤地依賴(lài)于地表法向量，對(duì)于這個(gè)地表起伏劇烈的模型，造成了梯度不正確，使反演過(guò)早地陷入了局部極值.

圖9 歸一化后的目標(biāo)函數(shù)收斂曲線虛線為傳統(tǒng)AST方法，實(shí)線為改進(jìn)的AST方法.Fig.9 The normalized objective function iteration curvesThe dashed and the solid lines represent the conventional and the improved AST methods respectively.

3.2 復(fù)雜起伏地表模型

下面使用一個(gè)相對(duì)復(fù)雜的模型來(lái)驗(yàn)證改進(jìn)方法的反演效果.圖10顯示了該模型5～35 km范圍內(nèi)的近地表附近的速度變化.該模型深部結(jié)構(gòu)比較簡(jiǎn)單，近地表速度橫向變化很大，還存在一些低速和高速異常體，30 km的橫向范圍內(nèi)地表最大起伏只有0.5 km，起伏還是比較平緩的.采用的縱向常梯度初始速度模型如圖11所示，模型離散為4001×151網(wǎng)格點(diǎn)，縱橫向網(wǎng)格間距均為10 m.

圖10 起伏地表真實(shí)速度模型Fig.10 The true complex topographic velocity model

在地表5～35 km范圍內(nèi)共激發(fā)751炮，相鄰炮點(diǎn)的水平間距為40 m.每炮都在地表設(shè)置500個(gè)檢波器，均采用雙邊接收方式，相鄰檢波點(diǎn)水平間距20 m，最大偏移距5 km.基于準(zhǔn)確速度模型、利用最短路徑法(Moser，1991)計(jì)算的走時(shí)作為“觀測(cè)”的初至波到達(dá)時(shí).分別使用依賴(lài)于地表法向量的傳統(tǒng)AST方法和改進(jìn)的不依賴(lài)于地表法向量的AST方法進(jìn)行試驗(yàn)，兩種方法在反演中均使用快速掃描法(Zhao, 2005)計(jì)算走時(shí)場(chǎng)以及伴隨場(chǎng).

兩種方法的反演結(jié)果分別見(jiàn)圖12和圖13.由于本模型的地表起伏相對(duì)比較平緩(30 km的橫向范圍內(nèi)地表最大起伏只有0.5 km)，傳統(tǒng)和改進(jìn)的AST方法都反演出了近地表速度的基本變化趨勢(shì)，但改進(jìn)方法比傳統(tǒng)方法能更精確地反演出近地表附近的高速和低速局部異常.另外，在傳統(tǒng)方法中有較嚴(yán)重的能量聚焦現(xiàn)象(圖12)，在同一速度層中，不同的速度參數(shù)呈條帶狀分布，橫向速度波動(dòng)較大，說(shuō)明在傳統(tǒng)AST方法伴隨方程的求解過(guò)程中，檢波點(diǎn)處的走時(shí)殘差反向傳播是不正確的，導(dǎo)致其梯度和模型修改量形態(tài)異常，而這些不合理現(xiàn)象在改進(jìn)的AST方法中得到了較好的改善(圖13).

圖11 初始速度模型Fig.11 The initial velocity model

圖12 傳統(tǒng)AST方法反演結(jié)果Fig.12 The inversion result by the conventional AST method

圖13 改進(jìn)AST方法反演結(jié)果Fig.13 The inversion result by the improved AST method

對(duì)兩種反演方法第80代的反演結(jié)果，抽取地表以下20、30、40、50 m深度的速度剖面(見(jiàn)圖14)，可以發(fā)現(xiàn)，無(wú)論是什么深度，改進(jìn)后AST方法(紅線)的反演精度均高于傳統(tǒng)方法(藍(lán)線)，改進(jìn)方法的反演結(jié)果更加接近真實(shí)模型.

圖14 第80代反演結(jié)果地表下不同深度的速度剖面對(duì)比真實(shí)模型(黑實(shí)線)、初始模型(黑虛線)、傳統(tǒng)AST方法反演結(jié)果(藍(lán)線)、改進(jìn)AST方法反演結(jié)果(紅線).Fig.14 The inverted lateral velocity slices with depth of (a) 20 m, (b) 30 m, (c) 40 m and (d) 50 m below the surface for the 80th iterationThe black solid line, the black dashed line, the blue line and the red line represent the true model, the initial model, the AST-DV inversion result and the AST-IV inversion result respectively.

從炮點(diǎn)在15、20、25 km處的三個(gè)單炮的走時(shí)殘差上也可以明顯看出(圖15)，改進(jìn)AST方法的反演結(jié)果的初至波走時(shí)殘差更小，說(shuō)明改進(jìn)后AST方法與傳統(tǒng)AST方法相比，反演精度得到了提高.

圖15 方法改進(jìn)前后反演結(jié)果的剩余走時(shí)對(duì)比(a) 炮點(diǎn)在15 km處的單炮; (b) 炮點(diǎn)在20 km處的單炮; (c) 炮點(diǎn)在25 km處的單炮. 真實(shí)模型(黑線)、傳統(tǒng)AST方法(藍(lán)線)、改進(jìn)AST方法(紅線).Fig.15 The residual first-arrival traveltime of the final inversion results for different shots with sources located at (a) 15 km, (b) 20 km and (c) 25 kmThe black line, the blue line and the red line represent the true model, the AST-DV inversion result and the AST-IV inversion result respectively.

4 結(jié)論

本文分析了目前傳統(tǒng)的伴隨狀態(tài)法初至波走時(shí)層析成像方法的不足：一是伴隨方程依賴(lài)于地表法向量不合理，二是無(wú)法合理處理井中觀測(cè)問(wèn)題.在此基礎(chǔ)上，論文進(jìn)一步提出了不依賴(lài)地表法向量的改進(jìn)的伴隨狀態(tài)法初至波走時(shí)層析成像方法，從理論上得到了不依賴(lài)于地表法向量的伴隨方程，克服了傳統(tǒng)方法中伴隨場(chǎng)依賴(lài)于地表法向量的缺陷，使得檢波點(diǎn)處的走時(shí)殘差可以正確地反傳播至地下，進(jìn)而可以得到更加合理的速度修正方向，提高了速度反演的精度.同時(shí)，提出的改進(jìn)方法還可以適應(yīng)地表和井中同時(shí)觀測(cè)的任意觀測(cè)系統(tǒng).兩個(gè)理論模型反演試驗(yàn)結(jié)果，說(shuō)明了改進(jìn)的AST方法的正確性和有效性.

附錄A 改進(jìn)的AST方法伴隨方程的推導(dǎo)

下面使用擾動(dòng)法(Nocedal and Wright, 2006)求解方程(10)，從而得到改進(jìn)的AST方法的伴隨方程.

即：

(A1)

(A2)

(A3)

其中，n為擴(kuò)展的邊界?Ω′的外法向單位矢量.

(A4)

和區(qū)域Ω′內(nèi)部的反傳播方程：

x∈Ω′.

(A5)

由于區(qū)域Ω及其邊界?Ω均位于擴(kuò)展的區(qū)域Ω′之內(nèi)，因此，無(wú)論在區(qū)域Ω內(nèi)還是在邊界?Ω上，伴隨場(chǎng)λ(x)均滿(mǎn)足(A5)式，即伴隨方程為：

x∈Ωand ?Ω.

該式即為改進(jìn)后的AST方法的伴隨方程(13)式.

可以看到，在實(shí)際需要計(jì)算的區(qū)域Ω+?Ω范圍內(nèi)，無(wú)論是在模型邊界上還是在模型內(nèi)部，改進(jìn)后的伴隨方程具有統(tǒng)一的形式，而且與地表的法向量無(wú)關(guān).