胥利文 陳 浩 張秀梅 王秀明

(1 中國科學(xué)院聲學(xué)研究所 北京 100190)

(2 中國科學(xué)院大學(xué) 北京 100049)

(3 北京市海洋深部鉆探測量工程技術(shù)中心 北京 100190)

0 引言

波動方程逆散射問題的研究是聲波在固體介質(zhì)中傳播和成像理論的重要研究內(nèi)容。一維波動方程逆散射問題的經(jīng)典解之一就是Marchenko 方程[1]，它將介質(zhì)邊界上單邊記錄的波場響應(yīng)和來自于介質(zhì)內(nèi)部聲源的散射響應(yīng)通過積分耦合，可以用來反演介質(zhì)內(nèi)部的構(gòu)造。最近，Broggini 等[2]展示了在不需要知道介質(zhì)信息的情況下，利用單邊記錄的反射響應(yīng)重構(gòu)介質(zhì)中虛擬點源的格林函數(shù)。Wapenaar 等[3]通過引入波場的因果性條件，將上述方法推廣到多維介質(zhì)中。van der Neut等[4]提出多維互相關(guān)理論，解釋了如何通過波場迭代求解多維Marchenko方程。Cui 等[5]將方程求解過程中定義的聚焦函數(shù)解釋為邊界源在負時刻的入射波場，分析了基于初始聚焦函數(shù)求取單程格林函數(shù)過程中幅度對結(jié)果的影響。除應(yīng)用在彈性介質(zhì)中，Marchenko自聚焦方法還可用于復(fù)雜的黏彈性介質(zhì)成像中[6?7]。

本文主要研究在一維介質(zhì)中，介質(zhì)表面到介質(zhì)內(nèi)任一點之間的格林場是如何通過耦合Marchenko方程組求解得到的，首先解釋了如何利用單程波場的互易原理推導(dǎo)耦合Marchenko 方程組；接著基于波場的因果性，設(shè)計時間窗算子迭代求解聚焦函數(shù)，并利用聚焦函數(shù)求解期望的格林函數(shù)；然后在求解過程中引入卷積理論，分別在數(shù)值解和解析解模型中計算期望的格林函數(shù)；最后基于解析解分析數(shù)值解中各個波形序列構(gòu)造的具體過程。

1 Marchenko方程

本文將聲學(xué)介質(zhì)限制在一維狀態(tài)，模型如圖1所示，藍實線表示聲速，紅虛線表示密度，位于介質(zhì)表面的黑色五角星和倒三角形分別表示聲源和接收器，在深度z=900 m 處的紅色五角星表示介質(zhì)中的虛擬點源，也是下文所述的聚焦函數(shù)的聚焦點。其中圖1(a)為狀態(tài)B，代表實際介質(zhì)；圖1(b)為狀態(tài)A，表示參考介質(zhì)。兩種介質(zhì)的不同之處在于參考介質(zhì)中用紅色五角星表示的聲源下方是均勻介質(zhì)，聲波入射之后無反射，相當于實際介質(zhì)的截斷版本。在兩種狀態(tài)下，對于兩個在時間軸上正向傳播的波場，可以通過卷積型互易原理聯(lián)系[8]：

對于時間軸上正向傳播的波場和反向傳播的波場，通過相關(guān)型互易原理聯(lián)系[8]：

式(1)～(2)中，pA和pB分別是兩種介質(zhì)中記錄的波場，“+”表示聲波朝深度增加的方向傳播，“?”表示聲波朝深度減小的方向傳播，本文分別用它們表示下行和上行的單程波場，()H為復(fù)共軛。等式(1)和等式(2)的左邊積分區(qū)域是沿著介質(zhì)表面?z0(z= 0 m)積分，右邊是沿著定義的聚焦平面?zi(z=900 m)積分。

在狀態(tài)A 中的介質(zhì)表面上方，如果布置一個向深度方向發(fā)射聲波的脈沖點聲源，則介質(zhì)表面記錄的下行波場為p+A=δ(x ?x0)，其中x0是聲源的水平坐標位置，在本文所示的一維介質(zhì)中x0= 0，因此p+A=δ(x)。同時介質(zhì)表面也會記錄來自于介質(zhì)內(nèi)部的反射響應(yīng)，即上行波場p?A=RA(z0,z0,t)。在定義的聚焦點z=900 m處記錄的下行波場是參考介質(zhì)的透射響應(yīng)p+A=TA(zs,z0,t)，同時因為在聚焦點下方是均勻介質(zhì)，無反射波，因此聚焦點處記錄的上行波場為0，即p?A=0。

圖1 1D 數(shù)值解模型Fig.1 1D numerical model

同理，在狀態(tài)B 中的介質(zhì)表面上方，也布置一個向深度方向發(fā)射聲波的點聲源，則介質(zhì)表面記錄的下行波場為p+B=δ(x)，上行波場是實際介質(zhì)的反射響應(yīng)，即p?B=R(z0,z0,t)。在聚焦點處，可以將記錄的全波場G(zs,z0,t)，分解為上行波場p?B=G?(zs,z0,t)和下行波場p+B=G+(zs,z0,t)，即G=G++G?。上述兩種狀態(tài)中的表達式R(z0,z0,t)、T(zs,z0,t)和G(zs,z0,t)，括號中的第一項表示接收器，第二項表示聲源，都表示聲源發(fā)射、接收器接收的波場響應(yīng)。

狀態(tài)A 中參考介質(zhì)的透射響應(yīng)TA(zs,z0,t)和反射響應(yīng)RA(z0,z0,t)是未知的，在狀態(tài)A 中定義一種新的函數(shù)：聚焦函數(shù)f1(z0,zs,t)，它的下行分量f1+(z0,zs,t)被定義為透射響應(yīng)的時反，表示從介質(zhì)表面z0發(fā)射聲波，聚焦于點zs處，所以可以將公式中的透射響應(yīng)TA(zs,z0,t)的時反用下行聚焦函數(shù)代替。同理定義聚焦函數(shù)的上行分量f?1(z0,zs,t)，表示聚焦點在zs處的向外擴散的波場，它是聚焦函數(shù)的下行分量通過參考介質(zhì)表面反射序列得到的響應(yīng)，所以將推導(dǎo)過程中的中間變量RA(z0,z0,t)f+1(z0,zs,t)用f?1(z0,zs,t)代替。聚焦函數(shù)表達式中的第二項zs表示此函數(shù)的聚焦點，這兩種聚焦函數(shù)的定義會方便后續(xù)的求解過程。

將上述不同狀態(tài)、不同位置處記錄的波場分量代入等式(1)和等式(2)，可以將“狀態(tài)B” 中深度為0 m 處記錄的反射響應(yīng)R(z0,z0,t)和深度為900 m 處的接收器接收的介質(zhì)表面聲源發(fā)射的單程格林函數(shù)G±(zs,z0,t)，通過定義在“狀態(tài)A”中的聚焦函數(shù)f1±(z0,zs,t)聯(lián)系起來。因為本文將方程求解過程限制在一維介質(zhì)，所以得到一維時域耦合Marchenko方程[3?5]：

其中，“?”表示卷積，f1±(z0,zs,t)表示聚焦于深度為zs處的聚焦函數(shù)，G±(zs,z0,t)為聲源位于z0、 在zs處接收的格林函數(shù)。為了書寫方便，下文將和R(z0,z0,t)分別簡寫為G+、G?、f1?、f1+和R。Marchenko 耦合方程組(3)和方程組(4)是欠定的，兩個方程，四個未知數(shù)：G+、G?、f1?、f1+，接下來將首先展示如何利用實際介質(zhì)中的反射響應(yīng)求取兩個定義在參考介質(zhì)中的聚焦函數(shù)。

2 聚焦函數(shù)重構(gòu)

位于zs= 900 m 處的脈沖點源激發(fā)的波場直達波將會在t=td(直達波傳輸時)到達介質(zhì)表面z0= 0 m 處的接收器，即格林函數(shù)中的所有序列將會在t≥td出現(xiàn)在接收記錄中。在t

依據(jù)上文和van der Neut 等[4]關(guān)于聚焦函數(shù)的分析，下行聚焦函數(shù)是參考介質(zhì)透射響應(yīng)的時反波場，它由?td時刻的首波和在?td與td之間的尾波序列構(gòu)成，將時間窗算子應(yīng)用于下行聚焦函數(shù)可以得到

而上行聚焦函數(shù)是在?td與td之間的波場，加窗后為

因此，對于耦合Marchenko方程組(3)和方程組(4)，將時間窗算子作用于等式兩邊最終可以得到

為了求解聚焦函數(shù)和，將介質(zhì)中虛擬點源到介質(zhì)表面的直達波的時反波場作為初始的下行聚焦函數(shù)，同時假定的初始值為0。將接收的反射響應(yīng)R(t)和初始聚焦函數(shù)的估計代入等式(7)，可以得到的首次更新，然后將代入等式(8)，得到的第一次更新值。接著，將第一次迭代求得的與求和，可以得到的全波場，再代入等式(7)就可以得到的第二次估計，然后再代入等式(8)得到的第二次更新，就這樣依次進行迭代求解，直到聚焦函數(shù)收斂且達到穩(wěn)定狀態(tài)。

假定圖1(a)所展示的數(shù)值解模型是無損的，且不考慮自由表面多次波的影響，利用時域有限差分算法[9]正演模擬反射響應(yīng)R(t)和初始聚焦函數(shù)。數(shù)值求解Marchenko 方程時，反射響應(yīng)R(z0,z0,t)需要子波反卷積，這里采用寬帶零相位子波作為聲源[10]，避免子波反卷積對反射響應(yīng)造成的誤差影響，而且較寬的頻帶將不會在時域上產(chǎn)生較長的時寬，子波時域圖如圖2(a)所示。初始聚焦函數(shù)可由宏觀速度模型計算[3]，這里直接運用圖1(a)的參考介質(zhì)正演。模擬初始聚焦函數(shù)的源子波設(shè)置為中心頻率為25 Hz 的雷克子波，子波時域圖如圖2(b)所示。

圖2 模擬反射響應(yīng)R 和初始聚焦函數(shù)的源子波Fig.2 Source wavelet for simulating the reflection response R and the initial focus function

圖3展示了得到的介質(zhì)表面記錄反射響應(yīng)和直達波，它們將作為迭代求解Marchenko 方程的初始輸入數(shù)據(jù)。通過分析圖3(a)中反射響應(yīng)的記錄和一維介質(zhì)的分層信息，得到圖4中的多次波序列解析解模型，圖4中從左到右的彩色箭頭分別表示圖3(a)中的前七個波序列。對于文中的一維模型，zi表示深度，下標i= 0,1,2,···對應(yīng)1D 模型的深度，其中z0是接收面，t1、t2、t3、t4分別代表聲波在各層的單程傳輸時，r01、r12、r23、r34、r45和τ01、τ12、τ23、τ34、τ45分別是按深度增加方向的各層的反射系數(shù)和透射系數(shù)，r54、r43、r32、r21、r10和τ54、τ43、τ32、τ21、τ10分別是按深度減小方向的各層的反射系數(shù)和透射系數(shù)。

圖3 迭代求解的初始輸入Fig.3 Initial input of iterative solution

圖4 1D 解析解模型Fig.4 1D analytical model

圖1(a)中的數(shù)值解模型中表面記錄的反射響應(yīng)可以展開為圖4解析解模型中各個序列的和，

為了更清晰地展示反射響應(yīng)中各個序列幅度對合成波場的影響，將等式(9)中用反射系數(shù)和透射系數(shù)表示的各個序列幅度用RAi(i=1,2,3···)表示，即

聚焦點位于深度z2和z3的中心處，聚焦點到表面接收點的直達波可以表示為

依據(jù)上文的解釋，將直達波的時反序列作為初始聚焦函數(shù)，即：

幅度τ21τ32用表示。

第一次迭代，將反射響應(yīng)和初始聚焦函數(shù)代入等式(7)可以求得f1?(t)，這個步驟相當于將反射響應(yīng)沿著時間軸逆向移動直達波的到時，即移動相位，加時間窗算子θ{·}后可得到這兩個序列構(gòu)造的f1?(t)：

由等式(11)得到的f1?(t)有兩個序列，幅度分別為因為聲波透射系數(shù)恒為正，依據(jù)圖1(a)模型中各層的聲阻抗的值，可得，可知f1?(t)中兩個序列的幅值一正一負。接著，將f1?(t)時反后代入等式(8)，相當于反射響應(yīng)分別與f1?(?t)中的兩個序列分別卷積(相移)再求和，加窗后可得

本模型中，通過計算可以發(fā)現(xiàn)(?t)的相位是由反射響應(yīng)中第一個序列與f1?(t)中的第二個序列構(gòu)造的，即將反射響應(yīng)中的第一個序列的相位移動，依據(jù)反射系數(shù)和透射系數(shù)的值可以判斷的幅值

圖5 第一次迭代后的結(jié)果Fig.5 The result after the first iteration

圖5是基于數(shù)值解模型求解的第一次迭代后的結(jié)果，圖中繪制了兩條紅虛線表示時間窗算子θ，兩條線在時間軸上鏡像對稱，窗算子θ去除正負窗口外的所有數(shù)據(jù)，而僅僅保存窗口之間的數(shù)據(jù)。利用時間窗算子截取窗內(nèi)部的序列只得到兩個序列(圖5(a))，這就是首次迭代更新的上行聚焦函數(shù)f1?(t)。與圖3(b)中的初始聚焦函數(shù)相比較，發(fā)現(xiàn)第一個序列的幅值是正數(shù)相乘的結(jié)果，而第二個序列是負數(shù)相乘的結(jié)果，與對等式(11)兩個序列的幅度分析結(jié)果一致。接著將這兩個序列時反，得到f1?(?t)，將它和反射響應(yīng)卷積可以得到圖5(b)的結(jié)果。兩個序列和反射響應(yīng)序列的卷積可以分為單個序列卷積的和，也相當于將反射響應(yīng)序列在時間軸上移動相位后再求和。窗算子作用后只留一個序列，即(?t)，時反后就可得到(t)，幅值為負數(shù)相乘的結(jié)果。數(shù)值求解的結(jié)果與等式(11)和等式(12)的結(jié)果相同。本質(zhì)上來說，來自反射響應(yīng)的序列隨著首波到達時間而移動，以構(gòu)建序列。由于首波是時反非因果的，因此它會將反射系列中的所有序列向負時間軸移動。其中的一些序列變得太多，以致變?yōu)榉且蚬蛄?。接下來又計算了第二、第三次迭代，得到具有相同相位而幅度不同的f?1(t)和(t)。

第二次迭代， 首先將第一次迭代求得的(?t)時反后與(t)求和可得

由等式(7)可得

接著將等式(14)得到的f?1 (t)代入等式(8)可得

第三次迭代，首先將等式(15)得到的(?t)與(t)求和：

由等式(7)可得

接著將等式(17)得到的f1?(t)代入等式(8)可得

將經(jīng)歷多次迭代后的f?1(t)與(?t)作為最終的聚焦函數(shù)，求取期望的格林函數(shù)。對比三次迭代結(jié)果，等式(11)、等式(14)和等式(17)得到的f?1(t)，發(fā)現(xiàn)每次迭代得到的下行聚焦函數(shù)的序列都只有兩列，且相位完全相同，而且第一個序列的幅度也是完全相同的，等式(14)中第二個序列的幅度是等式(11)得到的序列幅度基礎(chǔ)上再加一個第三次迭代的結(jié)果是第一次結(jié)果的基礎(chǔ)上對上一次迭代的幅度進行微調(diào)，而不影響相位。同理，對于等式(12)、等式(15)和等式(18)得到的(?t)，只有一個序列，幅度在每次迭代后進行微調(diào)。

數(shù)值模擬實驗中，總計迭代計算了六次，抽取其中第一、第五和第六次結(jié)果，在圖6中對比，不同的顏色分別代表不同的迭代結(jié)果。對于圖6(a)中的f?1(t)只有兩個序列，第一個序列幅度和相位完全重合，第二個序列的幅度有輕微的差別，將圖中第二個序列放大得到6(a)圖中左邊插圖展示，可以發(fā)現(xiàn)幅度差別很小，第六次迭代結(jié)果已經(jīng)收斂于第五次迭代。在圖6(b)中，只有一個序列，幅度差異很小，不同的迭代結(jié)果在到時和波形上匹配很好，數(shù)值解的圖形很好匹配了解析解的分析結(jié)果。

3 單程格林函數(shù)重構(gòu)

將在時間軸上與窗算子θ相逆的算子記為Θ，Θ去除了直達波到達前的所有序列，而不包括直達波本身?？捎脕砬蠼馕挥谒阕应韧鈬母窳趾瘮?shù)。將Θ置于Marchenko 方程組(3)和方程組(4)兩端，可以得到

圖6 不同迭代結(jié)果對比Fig.6 Comparison of different iteration results

和

通過等式(7)和等式(8)的三次迭代求解得(t)和f1?(t)，結(jié)合已知的R(t)和(t)，代入等式(19)和等式(20)，就可以得到單程格林函數(shù)。為了更加清晰地描述迭代過程的特點，只取迭代結(jié)果的前三個序列，具體解析解求解過程如下。

首先(t)的初始值為0，第一步可用R(t)和卷積后，加窗算子Θ求取上行格林函數(shù)，得到

正如上文所描述的，這個步驟相當于將反射響應(yīng)沿著時間軸逆向移動直達波的到時，即移動相位，加窗算子Θ作用后求取算子θ外圍的的波場序列。依據(jù)等式(12)求得的(t)，可以求得Θ{R(t)(t)}第一次迭代計算的結(jié)果：

對比等式(21)和等式(22)，可以發(fā)現(xiàn)等式(22)中前三個序列中沒有等式(21)中的第二個序列，同時對于等式(21)和等式(22)中的第一個序列，幅值，幅值，對于等式(21)中第三個序列和等式(22)中第二個序列的幅值，接著由等式(15)的結(jié)果，得到第二次迭代后的結(jié)果：

利用等式(18)，得到第三次迭代后的結(jié)果：

接著可以將等式(10)和等式(13)兩次求得的f1?(t)，結(jié)合已知的R({t)代入等式(20)，}首先可以分別得到三次迭代的ΘR(t)?f1?(?t)，具體結(jié)果如下。

第一次迭代：

第二次迭代：

第三次迭代：

對比等式(25)、等式(26)和等式(27)，三次結(jié)果的前三個序列的相位完全一樣，而幅度差異很小。接著可以基于等式(20)，利用已知的(?t)，減去上文求解的Θ{R(t)?f1?(?t)}，得到下行格林函數(shù)。(?t)的相位為與第一個序列的相位完全相同，所以相減操作是為了校正第一個序列的幅度，同時將后續(xù)的幅度乘以?1，使子波的幅值歸位。

圖7 格林函數(shù)序列不同迭代結(jié)果對比圖Fig.7 Comparison of different iteration results of Green’s function sequence

接下來，將上文基于Marchenko 方法迭代求解的上下行格林函數(shù)和正演模擬波場分離[8]的結(jié)果進行對比，其中紅虛線表示正演模擬的結(jié)果，藍實線表示由Marchenko方程求解結(jié)果，如圖8所示。

從圖8中，發(fā)現(xiàn)藍實線表示的基于Marchenko方程求解的上行和下行格林函數(shù)迭代求解的結(jié)果，和正演模擬的結(jié)果在波形和相位上吻合很好，說明了一維Marchenko 方程準確地分解了虛擬源位置接收的波場。如果將上行格林函數(shù)和下行格林函數(shù)求和，可以得到點源激發(fā)的全波場。

圖8 Marchenko 方程求解結(jié)果與正演結(jié)果對比Fig.8 Comparison of the results of the resolved Marchenko equation and the directly model

4 結(jié)論

本文通過一維聲學(xué)介質(zhì)的數(shù)值解和解析解模型對耦合Marchenko 求解單程格林函數(shù)過程進行了詳細的討論和說明?？梢缘玫揭韵挛鍌€結(jié)論：

(1)在欠定方程組不可求解的情況下，可以基于聚焦函數(shù)和格林函數(shù)的波場因果性質(zhì)，通過引入時間窗算子迭代求解耦合方程，時間窗算子內(nèi)部計算聚焦函數(shù)，外部計算格林函數(shù)；

(2)通過對不同迭代次數(shù)的聚焦函數(shù)分析，可以發(fā)現(xiàn)在第一次迭代后，期望的序列在相位上已經(jīng)很準確，而后續(xù)的迭代過程是對幅度的微小調(diào)整，而且聚焦函數(shù)內(nèi)部波形序列的相位差，在后續(xù)的步驟中可以去除由層間多次波產(chǎn)生的偽波序列；

(3)利用初始聚焦函數(shù)計算的上行格林函數(shù)估計包含了偽波，這些偽波需要利用R(t)(t)的迭代穩(wěn)定值去除；對于下行格林函數(shù)，可以發(fā)現(xiàn)在第一次迭代后，期望的序列在相位上已經(jīng)很準確，而后續(xù)的迭代過程是對幅度的微小調(diào)整；

(4)從解析解等式中可以發(fā)現(xiàn)，基于Marchenko方法迭代求解單程格林波場，反射響應(yīng)和直達波的幅度構(gòu)造了聚焦函數(shù)和上下行格林函數(shù)的幅度，因為Marchenko 方法中反射響應(yīng)是子波反卷積的結(jié)果，所以準確的子波反卷積操作對結(jié)果很重要；

(5)迭代求解Marchenko 方程本質(zhì)上是通過初始聚焦函數(shù)對記錄的單邊反射響應(yīng)進行相位調(diào)制，接著利用聚焦函數(shù)的作用，使反射響應(yīng)中的序列逐步移位到介質(zhì)中虛擬接收器接收的波場響應(yīng)序列。