石振明 盧崔燦 劉 鎏 王 堯 彭 銘

(①同濟大學地下建筑與工程系，上海 200092，中國)(②同濟大學巖土及地下工程教育部重點實驗室，上海 200092，中國)(③佛羅里達大學土木與海岸工程系，佛羅里達 32611，美國)

0 引 言

巖溶是水對可溶性巖石(碳酸鹽巖、石膏、巖鹽等)進行以化學溶蝕作用為主，物理侵蝕為輔的地質(zhì)作用現(xiàn)象(張咸恭等，2000)。全球巖溶塌陷廣泛分布，而我國巖溶地質(zhì)面積約占國土面積的1/3(陳學軍等，2019)。復雜巖溶地質(zhì)環(huán)境中存在無規(guī)則巖溶分布與填充情況、復雜的巖溶水賦存條件、起伏的土巖界面(姚邦杰等，2019)，給大型工程的選址、設計、施工和運營造成極大的風險，探明地下巖溶分布具有重要意義。

目前，應用于工程中的巖溶探測地球物理方法可分為地表物探方法和孔中物探方法，前者包括高密度視電阻率法(Ei-Qady et al.，2005; 底青云等，2018)、地震反射法(Li et al.，2017)和高頻面波法(夏江海等，2015)；后者包括跨孔彈性波CT法(朱文仲等，2008)、鉆孔雷達法(Tan et al.，2012)、管波測試法(李學文等，2005)和樁底溶洞聲吶探測法(石振明等，2016)。跨孔彈性CT法由于結(jié)合鉆孔，其探測深度和精度較符合巖土勘察的需求，是適合巖溶強發(fā)育區(qū)的勘察方法。

基于初至旅行時的層析成像方法，由于其高效性，是目前近地表建模應用最廣泛的方法。經(jīng)典的初至層析反演方法基于射線追蹤方程，在反演中需要顯式計算炮點到檢波點的旅行時及射線路徑(Stewart，1991)，反演通常采用最小奇異值分解(LSQR)(Zelt et al.，1998)，迭代重構(SIRT)(Yang et al.，2009)等方法。而基于程函方程的層析反演方法僅需要計算初至時間場而不需要計算射線路徑，借助伴隨狀態(tài)方法求取模型更新量，在復雜區(qū)域理論上具有更好的反演效果，更適合巖溶強發(fā)育區(qū)的探測(Plessix，2006)。

本文針對跨孔地震波巖溶探測，展開了基于程函方程的層析反演的數(shù)值模擬研究。利用高階時域有限差分解彈性波方程，考慮多種巖溶發(fā)育工況，模擬現(xiàn)場實測數(shù)據(jù)，并對模擬數(shù)據(jù)進行基于程函方程的層析反演，驗證了該方法在跨孔巖溶探測中的有效性，并討論了其局限性和未來研究方向。

1 數(shù)值研究方法

基于程函方程的初至層析反演是一種以程函方程作為反演中的正演方法的成像技術。本文的數(shù)值研究思路如圖1，先設置真實模型，用時域有限差分計算得地震記錄，模擬現(xiàn)場觀測數(shù)據(jù)，從地震記錄中提取初至旅行時；然后設置初始模型，基于程函方程計算時間場，計算旅行時殘差，并更新模型，將更新模型作為新的初始模型代入，最后重復迭代算法直至滿足迭代終止條件，得到最終的反演地層速度模型，用于巖溶地質(zhì)的識別。

圖1 基于程函方程的跨孔地震巖溶探測反演研究思路

1.1 正演方法

根據(jù)以下二維彈性介質(zhì)中的彈性波方程模擬跨孔地震波探測過程：

(1)

式中：x，y為兩個方向的直角坐標；t為時間(s)；ρ為密度(kg·m-3)；vx，vy分別為質(zhì)點振動速度在x，y兩個方向上的分量(m·s-1)；σxx，σyy，σxy為應力分量(N·m-2)；fx，fy為體力在x，y兩個方向上的分量(N·m-3)；σxx0，σyy0，σxy0為表面力分量(N·m-2)；λ，μ為拉梅常數(shù)(N·m-2)(Landau et al.，1986)。

利用高階交錯網(wǎng)格的時域有限差分解彈性波方程，將應力分量和速度分量分別定義在兩套網(wǎng)格上，并引入半網(wǎng)格點，在半網(wǎng)格點上進行空間導數(shù)的計算，即在網(wǎng)格上進行線性應力速度關系的離散化和動量方程的離散化。正演的初始時間所有網(wǎng)格點速度和應力均為0。為提高計算效率，必須引入完全匹配層(PML)作為吸收邊界，在吸收層內(nèi)部設置衰減因子對波場進行衰減。

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

旅行時層析反演方法是通過迭代計算尋找正演旅行時與觀測旅行時之間差異最小的地層參數(shù)模型(郭振波等，2019)，其目標函數(shù)可定義為正演和觀測旅行時之差的最小二乘如式(2)：

(2)

t(x)·t(x)=s2(x)

(3)

式中：t(x)為波場從源點到空間任意點x的旅行時；s(x)為x處的慢度。利用迎風有限差分對式(3)進行離散之后的形式如式(4)：

(Dxts)(Dxts)+(Dyts)(Dyts)=diag(m)m

(4)

生產(chǎn)力的發(fā)展是社會文明形式的最終決定力量。換句話說，每一種社會文明都不是憑空而來無緣而去，說到底還是有社會生產(chǎn)力的發(fā)展水平所決定的。生產(chǎn)力大發(fā)展過程實際上是生產(chǎn)力不斷積累過程，他的發(fā)展是川流不息的，他不會因為某個階段的消滅而消滅，也不會因為某種社會形態(tài)的消亡而消亡。當生產(chǎn)力發(fā)展到一個更高的階段，原有的舊的社會文明形態(tài)已經(jīng)不適應社會生產(chǎn)力的發(fā)展，于是就將逐漸被淘汰，社會就將誕生新的文明形式以適應生產(chǎn)力的發(fā)展，同時這種新的文明形式將隨著生產(chǎn)力的繼續(xù)發(fā)展而不斷成長壯大。生態(tài)文明的誕生和發(fā)展，就是生產(chǎn)力發(fā)展到新的階段的必然要求。

(5)

將式(5)帶入目標函數(shù)式(2)后對介質(zhì)參數(shù)m求一階導數(shù)如式(6)：

(6)

若先對式(4)中m在等式兩側(cè)求導，可先得到：

(7)

將式(7)代入式(6)，可得基于程函方程的層析反演的目標函數(shù)的一階導數(shù)為：

(8)

如想得到目標函數(shù)O(m)的最小值，其一階導數(shù)應等于0，即需要找到適合的介質(zhì)參數(shù)m使得g=0; 可以將式(2)在m0處進行二階泰勒展開后對m求導，并令其等于0。可得到：

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

(9)

式中：H為目標函數(shù)對介質(zhì)參數(shù)的二階導數(shù)，即Hessian矩陣。

式(8)結(jié)合式(9)經(jīng)整理(Taillandier et al.，2009)可得其等價形式非線性迭代求解式(10)。

(10)

利用高斯-牛頓法求解式(10)，在求解過程中Hessian矩陣近似為單元矩陣(即最速梯度法)。從而得到介質(zhì)參數(shù)m更新量，在反演迭代過程中，速度模型不斷更新，程函正演的初至旅行時與從觀測記錄中提取的初至旅行時的殘差逐漸減小，當殘差小于設定值時即滿足迭代終止條件時，停止迭代，此時的速度模型即最終反演結(jié)果。

2 數(shù)值試驗

巖溶最基本的兩種發(fā)育形態(tài)為溶洞和斷層，因此數(shù)值試驗設置了單一溶洞、單一斷層、溶洞斷層組合、多溶洞組合共4組工況。4種不同的地質(zhì)體組合的工況以均質(zhì)地層為背景，在地層中有兩個深度為10im，其相鄰水平間隔為10im的鉆孔，地層的材料參數(shù)如表1。在模型的左側(cè)鉆孔中布置了一排豎向間距為0.5im的發(fā)射器，共20個，震源是中心頻率為3ikHz的Spike波。模型的右側(cè)鉆孔中布置了一排豎向間距為0.2im的接收器，共50個。

表1 地層參數(shù)設置

2.1 工況1：單一溶洞

圖2為單一溶洞正反演結(jié)果。其中：圖2a表示正演模擬的真實模型；由圖2a用時域有限差分計算可得地震記錄如圖2c、圖2d，分別代表從上到下第1炮(深度為0.5im)、第20炮(深度為10im)x分量的地震記錄，以此模擬現(xiàn)場觀測記錄；從地震記錄中提取初至旅行時，得到圖2e中的觀測旅行時；圖2b表示反演的初始模型；將程函方程作為反演中的正演方法，由圖2b可計算得到圖2e中的正演旅行時，比較圖2a 與圖2b，發(fā)現(xiàn)在反演前圖2e中的觀測旅行時和正演旅行時差別較大；在反演過程中，通過多次迭代，更新速度模型，并更新正演旅行時，可以看到圖2h中標準化殘差隨著迭代次數(shù)增加而不斷減小，最終滿足標準化殘差小于0.005的迭代終止條件，得到圖2f迭代結(jié)束后旅行時，圖2g最終反演結(jié)果。

圖2 工況1：單一溶洞試驗結(jié)果圖

對比圖2a 和圖2g可以看出，反演結(jié)果揭示了鉆孔間深度6～9im處有一溶洞發(fā)育，其填充物速度明顯低于背景地層?；诔毯匠痰姆囱莼窘沂玖巳芏吹奈恢谩⒋笮『托螤?，反演迭代次數(shù)為10次，收斂迅速。就溶洞的形狀而言，反演的縱向分辨率略高于橫向分辨率。同時，洞中填充物的波速較真實模型偏高。

2.2 工況2：單一斷層

單一斷層的反演數(shù)值研究過程如圖3a～圖3h所示。在真實模型中設置一傾斜的線性斷層，深度發(fā)育在1～4.5im，斷層內(nèi)填充含水的低速體。反演經(jīng)過了15次迭代達到了收斂。對比圖3a 和圖3g可以看出，反演結(jié)果成功揭示了斷層的方向和位置，反演結(jié)果中斷層填充物的波速較真實模型偏高。與工況1的反演結(jié)果相比，圖3g中有較多的偽像，相比溶洞，初至旅行時層析反演對斷層較不敏感。

圖3 工況2：單一斷層試驗結(jié)果圖

從理論上解釋，這是因為經(jīng)過斷層的最短射線路徑相對溶洞較小，盡管斷層的反射波對整個地震記錄上產(chǎn)生了很大變化，但是由斷層產(chǎn)生的初至旅行時變化較小，導致初至層析反演過程中容易陷入局部最優(yōu)解的問題。

2.3 工況3：溶洞斷層組合

為了研究基于程函方程的跨孔層析反演方法對多個異常體同時成像的效果，工況3組合了工況1和工況2的兩種巖溶地質(zhì)體，在工況3模型中設置了相同的溶洞和斷層。通過對比圖2、圖3和圖4可以看出工況3的地震記錄、反演結(jié)果與工況1和工況2疊加的結(jié)果相似，溶洞與斷層的存在對各自反演結(jié)果的相互干擾較小。同時工況3的迭代次數(shù)并沒有隨著模型變得復雜而明顯增加，收斂速度與前兩個較簡單工況相當，表明基于程函方程的跨孔層析反演方法對于復雜地質(zhì)條件的反演具有較好的穩(wěn)定性。

圖4 工況3：溶洞斷層組合試驗結(jié)果圖

同時也注意到圖4a 中，真實模型中斷層中填充物的波速低于溶洞中填充物的波速，而在圖4g中，反演的斷層中填充物的波速卻高于溶洞中填充物的波速，更接近地層背景波速，這也驗證了該初至旅行時的反演對斷層較不敏感，反演成像精度沒有溶洞高。

2.4 工況4：多溶洞組合

串珠狀溶洞是典型的復雜巖溶地質(zhì)現(xiàn)象，往往在單一鉆孔中會揭示多個溶洞的存在。工況4設置了多個溶洞，并考慮了不同填充物和形狀，以討論該方法對多個溶洞不同方面的成像效果。

其數(shù)值研究過程如圖5所示，對比圖5a與圖5g發(fā)現(xiàn)，反演結(jié)果有效區(qū)分了兩個不同溶洞的深度和水平位置、形狀、大小和填充物速度的高低。反演的迭代次數(shù)有所增加(16次)，但依然可以在短時間內(nèi)完成。同時對比圖2g、圖5g 可以看出兩個溶洞的存在對反演結(jié)果的相互干擾較小。反演依然存在低速異常體反演速度比真實速度高的問題。且橫向分辨率明顯沒有縱向分辨率高，這與射線密度有一定關系。

圖5 工況4：多溶洞組合試驗結(jié)果圖

根據(jù)4個工況的數(shù)值試驗結(jié)果，可知鉆孔間的溶洞和斷層等速度異常體對地震記錄影響很大，其中溶洞對初至旅行時的影響更為明顯。初至旅行時類的反演方法因為只提取初至信息，忽略了波形信息，反演速度較快，在不到20次迭代中均可以將標準化殘差降低至0.005以下，得到孔間巖溶異常體的分布。但由于僅僅提取初至時間，沒有充分利用反射波和透射波信息，對于突變型的地質(zhì)異常體，往往反演為一個速度漸變的異常體，這也導致：反演速度比真實速度偏高，更接近背景波速；其異常體的實際橫向尺寸偏大。但考慮到基于程函方程的跨孔層析反演方法在多個巖溶地質(zhì)體共存的復雜地質(zhì)條件下依然可以保持高效和穩(wěn)定，該方法可以為實際工程提供重要勘察信息。

3 結(jié) 論

本文針對跨孔地震波巖溶探測，通過數(shù)值模擬試驗，設置了4組巖溶發(fā)育工況，實現(xiàn)了基于程函方程的初至旅行時層析反演，得到以下結(jié)論。

(1)本文針對基于跨孔地震波巖溶探測，提出了基于程函方程的層析反演方法。并通過時域有限差分模擬真實探測數(shù)據(jù)，利用基于程函方程的初至層析反演對模擬數(shù)據(jù)進行成像，有效地揭示了鉆孔間的巖溶地質(zhì)體。

(2)基于程函方程的初至旅行時層析反演方法對溶洞位置、大小，斷層的方位反演較為準確，對于復雜地質(zhì)條件，多溶洞和裂隙組合反演較為穩(wěn)定，相互之間干擾很小。能定性區(qū)分溶洞填充物的速度與地層背景波速的高低，為工程處理巖溶發(fā)育提供依據(jù)。

(3)基于程函方程的初至層析反演方法有一定局限性，巖溶異常體中的反演速度較真實速度偏高，縱向分辨率比橫向分辨率高。相比溶洞，該方法對斷層較不敏感，因為斷層對初至旅行時的影響較小。

(4)基于程函方程的反演雖然未顯性地計算射線路徑，但可以通過等時線求垂線的方式求解，后續(xù)研究可探究射線密度對反演效果的影響，并考慮初始模型對探測精度的影響。同時，初至層析成像方法可以與全波形反演(FWI)結(jié)合，以充分利用現(xiàn)場探測數(shù)據(jù)信息，得到更好的反演精度。