改進的伴隨狀態法初至波走時層析成像方法

2021-03-08 09:45:36董良國張建明韓佩恩

地球物理學報 2021年3期

董良國，張建明，韓佩恩

同濟大學海洋地質國家重點實驗室，上海 200092

0 引言

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

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

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

基于程函方程的AST方法最早由Sei和Symes(1994)所提出，自此之后，幾乎所有的關于AST方法的理論方法及應用的文獻，都是基于面積分來定義目標函數，由此得到的伴隨方程都依賴于地表的法向量，本文稱這種方法為傳統AST方法.

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

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

我們已經發現，這種傳統AST方法無論在理論方法上還是在應用過程中還存在一定的問題(Han et al., 2019)，本文的目的是對這些問題進行進一步的討論.首先對這種傳統AST方法簡要介紹之后，從理論模型試驗以及地震波傳播物理規律的角度指出了這些問題的存在，并提出了改進的方法，即不依賴于地表法向量的改進的AST走時反演層析成像方法，并通過兩個模型反演試驗證明了改進方法的正確性和有效性.至于傳統AST方法存在的問題的嚴格數學證明，是我們下一步的工作目標.

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

1 目前傳統AST方法存在的問題

1.1 傳統AST方法簡介

觀測和理論走時達到最佳匹配是走時反演的基本原則.目前傳統伴隨狀態法初至波走時層析普遍采用面積分來定義目標函數(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表示檢波點的位置，檢波點位于區域Ω的邊界?Ω上，Tobs(r)表示r處的實際觀測走時，c(x)表示空間位置x處的介質速度.t(r)表示在當前速度模型下空間位置r處計算的理論走時，在高頻假設下，走時場t(r)滿足程函方程：

(2)

利用伴隨狀態法，可以確定上述目標函數(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)

其中，伴隨場λ(x)滿足：

(4)

(5)

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

利用(4)式，根據檢波點處的走時殘差設置檢波點處的伴隨場λ的初值后，再利用方程(5)將走時殘差反向傳播即可獲得伴隨場λ(x)，進而通過公式(3)就可以得到目標函數對速度的梯度，然后通過線性搜索算法可以很方便地獲得一個優化的步長αn(其中n為迭代次數)，這樣就可以利用下面的局部優化算法進行速度迭代，得到最終的速度反演結果.

(6)

可見，這種AST方法的核心是利用伴隨方程計算伴隨場λ(x)，關鍵是伴隨方程.自從Sei和Symes(1994)提出AST走時反演方法之后，有關AST方法的所有文獻中，所用的伴隨場方程的邊界條件均是(4)式，該式依賴于邊界的法向量.

1.2 傳統AST方法存在的問題

下面具體分析上述傳統AST方法存在的問題.

問題一：傳統AST方法采用面積分來定義目標函數，決定了其無法處理檢波器在模型內部的情況.

方程(1)采用面積分來定義目標函數，就預先假定了所有檢波器均位于模型表面，這就造成了傳統AST方法只適用于檢波點在邊界上、無法考慮檢波器在模型內部的情況.當進行井中觀測以及井間觀測時，檢波器在模型內部，這時方程(1)的面積分定義目標函數的方式就不適用，伴隨方程(4)中的所謂的邊界法向當然無從談起.可見，對于檢波點在地下的情況(如井中觀測)，這種傳統的AST方法是不適用的.

問題二：傳統AST方法中，伴隨場的計算依賴于模型表面的法向量，這是有悖射線理論實際的.

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

圖1 不同地表法向量的對稱接收點試驗的真實模型Fig.1 The true velocity model with two symmetrical receivers but different surface normal vectors

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

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

圖4 傳統AST方法計算的梯度Fig.4 The gradient obtained by the conventional AST method

問題三：當地表和井中同時觀測時，傳統AST方法無法正確計算梯度.

2 改進的AST方法

若使用傳統AST方法中的方程(1)來定義走時層析目標函數，檢波點只能位于區域Ω的邊界?Ω上.我們將?Ω′面定義為?Ω面向外的一個拓展面，拓展的空間部分充填空氣，速度為340 m·s-1.這樣，原來地表上以及井中的檢波點均變為區域Ω′的內部點，如圖5.由于實際地表之上本來就是空氣層，并且在區域Ω的邊界?Ω為凸邊界的條件下，這樣的拓展并不會改變區域Ω內的地震波走時.

圖5 計算區域的拓展示意圖Fig.5 The extended sketch map of the computational domain

現在基于拓展的新空間Ω′及其邊界?Ω′，利用體積分定義伴隨狀態法初至波走時層析目標函數：

(7)

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

(8)

在最優化過程中，走時場變量t、伴隨變量λ和速度c相互獨立，所以當目標函數達到極小值時，目標函數仍滿足以下三個關系：

(9)

(10)

(11)

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

(12)

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

我們使用擾動法求解(10)式，得到了改進的AST方法的伴隨方程(推導過程見附錄A)：

x∈Ωand ?Ω.

(13)

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

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

也就是說，改進的AST方法適用于各種地表和井中同時觀測方式，地下和井中觀測數據的處理方式得到了統一，從理論方法上避免了Waheed等(2016)處理井中觀測時的井中和地表觀測數據對梯度貢獻不一致的問題.無論是內部點還是邊界上，改進方法的伴隨方程都遵從統一形式，伴隨場的計算與模型邊界的法向量無關，不需要區分檢波點位置處的地表法向量，這樣更符合實際物理意義，保證了梯度場的正確性.

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

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

3 模型試驗

下面，通過兩個理論模型試驗來說明改進的AST方法的正確性及反演效果.

3.1 簡單起伏地表模型

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

圖7 簡單起伏地表速度模型及改進前后AST反演結果(a) 簡單起伏地表速度模型； (b) AST-IV反演結果； (c) AST-DV反演結果.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網格點，縱橫向網格間距均為10 m，采用速度值為2000 m·s-1的常速模型作為反演的初始模型.分別使用改進的不依賴于地表法向量的AST方法(AST-IV)和依賴于地表法向量的傳統AST方法(AST-DV)進行試驗，反演結果分別見圖7b和圖7c.其中，AST-IV和AST-DV分別迭代了693代和575代，它們的目標函數都不再下降.

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

圖8展示了改進AST方法與傳統AST方法不同迭代次數(第1、20、100、300次)的歸一化梯度值.可以發現，隨著迭代的進行，改進AST方法的梯度形態逐漸變得與高速異常體形態一致，而傳統AST方法的梯度沒有得到本質的改善，反演過早就陷入了局部極值.這進一步證實了傳統AST方法在伴隨方程的求解過程中，檢波點處的走時殘差反向傳播是不正確的，導致其梯度和模型修改量形態異常，造成該問題的根源就是來自于傳統AST方法的伴隨場錯誤地依賴于地表法向量，當地表起伏劇烈時，這個問題就很突出.而改進的AST方法的伴隨場與地表法向量無關，它可以有效修正梯度，從而最終能得到正確的反演結果.

圖8 傳統與改進的AST方法不同迭代次數的歸一化梯度從上到下分別為第1、20、100、300代的梯度(左側為傳統的AST方法，右側為改進的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.

從本試驗的目標函數收斂曲線(圖9)上可以發現，改進AST方法不僅收斂速度明顯快于傳統AST方法，且反演結果的走時殘差遠小于傳統AST方法的走時殘差.這主要是由于傳統AST方法的伴隨場錯誤地依賴于地表法向量，對于這個地表起伏劇烈的模型，造成了梯度不正確，使反演過早地陷入了局部極值.

圖9 歸一化后的目標函數收斂曲線虛線為傳統AST方法，實線為改進的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 復雜起伏地表模型

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

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

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

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

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

圖12 傳統AST方法反演結果Fig.12 The inversion result by the conventional AST method

圖13 改進AST方法反演結果Fig.13 The inversion result by the improved AST method

對兩種反演方法第80代的反演結果，抽取地表以下20、30、40、50 m深度的速度剖面(見圖14)，可以發現，無論是什么深度，改進后AST方法(紅線)的反演精度均高于傳統方法(藍線)，改進方法的反演結果更加接近真實模型.

圖14 第80代反演結果地表下不同深度的速度剖面對比真實模型(黑實線)、初始模型(黑虛線)、傳統AST方法反演結果(藍線)、改進AST方法反演結果(紅線).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.

從炮點在15、20、25 km處的三個單炮的走時殘差上也可以明顯看出(圖15)，改進AST方法的反演結果的初至波走時殘差更小，說明改進后AST方法與傳統AST方法相比，反演精度得到了提高.

圖15 方法改進前后反演結果的剩余走時對比(a) 炮點在15 km處的單炮; (b) 炮點在20 km處的單炮; (c) 炮點在25 km處的單炮. 真實模型(黑線)、傳統AST方法(藍線)、改進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 結論

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

附錄A 改進的AST方法伴隨方程的推導

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