樊亞楠，李貅*，戚志鵬，魯凱亮

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安 710054 2 長安大學(xué)地球物理場多參數(shù)綜合模擬實驗室(中國地球物理學(xué)會重點實驗室)，西安 710054

0 引言

瞬變電磁法采集的信號為異常體引起的純二次場，具有成本低、裝置輕便、工作效率高等優(yōu)勢，因而被廣泛應(yīng)用于地下水探測(Fitterman and Stewart，1986；Foley et al.，2016)、金屬礦勘探(Nabighian and Asten，2002；Vallée et al.，2011；底青云等，2019)、含水采空區(qū)探測(Li D S et al.，2019)以及隧道超前預(yù)報(Li Y et al.，2019；Liu et al.，2020)等領(lǐng)域中(孫懷鳳等，2021).然而，瞬變電磁場滿足時間域擴散方程，刻畫的是電磁場的感應(yīng)擴散特征，不易于對電性界面成像(李貅等，2010，2012)，而波動場滿足波動方程，反映的是波的傳播特性，易于對地質(zhì)界面成像(薛國強等，2006).因此，通過波場反變換將瞬變電磁場轉(zhuǎn)換為虛擬波場，即可借助地震勘探中成熟的成像方法進行瞬變電磁擬地震成像(薛國強等，2011；戚志鵬等，2013；Xue et al.，2013,2019；鐘華森等，2016).

在地震勘探中，Born近似算法可以解決積分方程解非線性的問題(Cohen and Bleistein，1977，1979；黃聯(lián)捷和楊文采，1991)，在無初始速度模型的條件下，可以直接得到深度剖面，克服了偏移算法需要事先給定初始速度模型的缺陷(丁科和宋守根，2004)，因此，該方法被廣泛應(yīng)用于地震勘探成像領(lǐng)域.Cohen和Bleistein(1979)、Bleistein和Cohen(1982)、Bleistein等(2001)假定反射波足夠弱，針對地震反演中的二維速度變化問題將非線性問題進行線性化處理，率先建立了Born近似逆散射理論.Beylkin(1984)利用Born近似小擾動理論，將反演問題線性化，引入傅里葉算子求解第一類Fredholm積分方程.Beylkin和Burridge(1990)提出了基于單散射的逆散射算法，進一步驗證了Born近似算法在對地震數(shù)據(jù)進行線性化時所起的重要作用；黃聯(lián)捷和楊文采(1991)推導(dǎo)出了2.5維聲波方程的逆散射表達式，利用傅里葉變換快速獲得速度擾動量，使Born近似算法的適用性得以提高.丁科和宋守根(2004)利用逆散射序列對含有地震多次波信息的數(shù)據(jù)進行奇性界面反演，使該方法不僅適用于小擾動量介質(zhì)模型的反演，而且對于大擾動量地質(zhì)模型同樣適用；Beylkin(1985)、Mao等(2013)、Ouyang等(2015)提出了二階Born近似的非線性反演算法，克服了一階Born近似僅適用于弱散射介質(zhì)的局限性.

在電法勘探領(lǐng)域，樊亞楠等(2019)在波場反變換基礎(chǔ)上，利用Born近似成像算法實現(xiàn)了對地質(zhì)界面快速成像的目的.由于Born近似算法在弱散射條件下，忽略散射場使方程線性化，導(dǎo)致對大擾動量的成像會產(chǎn)生較大的誤差.因此本文在波場反變換的基礎(chǔ)上，借鑒地震勘探上的二階Born近似算法進行擬地震成像研究，使得在強散射介質(zhì)中能夠準確定位地質(zhì)界面的位置和形態(tài).

本文基于波場反變換，將波動方程中的總波場和總波速分為兩部分，并運用Green定理推導(dǎo)出了一維常背景速度地質(zhì)模型的Born近似表達式，然后考慮散射序列二階項對成像結(jié)果的影響，根據(jù)遞推公式，推導(dǎo)出了適用于大擾動量模型的二階Born近似成像算法.通過對三層、四層模型的解析解以及二維、三維模型的計算可知，Born近似算法對地下電性界面成像結(jié)果較差，而二階Born近似算法能夠較為準確地定位地質(zhì)界面的位置和形態(tài).因此，建立在波場反變換基礎(chǔ)上的二階Born近似成像算法可以實現(xiàn)對大擾動地質(zhì)模型界面的準確定位和快速成像.

1 基本原理

1.1 Born近似算法

由文獻(Lee et al.，1989)可知，瞬變電磁擴散場與虛擬波動場之間的關(guān)系式為

(1)

其中Hz(r,t)表示磁場強度的z分量，單位為A·m-1，其滿足的擴散方程為

(2)

U(r,τ)表示虛擬波場，其滿足的時間域波動方程為

(3)

虛擬波場滿足的波動方程在頻率域的表達式可以寫為

(4)

為確保波動方程在無邊界介質(zhì)中解的唯一性，防止波從無窮遠處向源點傳播，設(shè)方程(4)滿足如下的索莫非輻射邊界條件：

(5)

根據(jù)擾動理論，可將地下變化的波速分為兩部分：全局性大尺度緩慢連續(xù)變化的背景波速v0(r)和局部快速變化的擾動速度α(r)，為了保持與波動方程形式相一致，可以將兩者關(guān)系表示為如下的形式：

(6)

同時，緩慢連續(xù)變化的背景波速v0(r)不產(chǎn)生新的震相，只與入射波場有關(guān)，將這部分波場記為背景波場ui(r,rs,ω)；局部快速變化的波速擾動產(chǎn)生新的震相，與入射波場和擾動波速有關(guān)，把入射波、繞射波和隨機散射波統(tǒng)稱為散射波場，記為us(r,rs,ω)，因此

u(r,rs,ω)=ui(r,rs,ω)+us(r,rs,ω).

(7)

將式(6)、(7)代入式(4)中，可以得到：

(8a)

+us(r,rs,ω)],(8b)

為求解方程(8b)，引入Green函數(shù)G(r,rg,ω)，其滿足方程：

(9)

運用Green定理，經(jīng)過一系列推導(dǎo)可得：

+us(r,rs,ω)]d3r.

(10)

由式(10)可知，界面速度擾動量α(r)和散射波場us(rg,rs,ω)之間存在非線性的關(guān)系，要準確求解速度擾動量α(r)，必須對其進行線性化近似，考慮當(dāng)擾動量較小時，即us(r,rs,ω)?ui(r,rs,ω)，忽略右端項中的散射場us(r,rs,ω)，這時可以得到：

(11)

其中ui(r,rs,ω)=G(r,rs,ω)，因此：

(12)

式(12)即為虛擬波動場關(guān)于電性界面速度擾動量的Born近似表達式.

在一維背景速度為常數(shù)、零偏移距條件下，即v0(x)=v0=const，xs=xg，式(12)可以寫為：

(13)

根據(jù)文獻(Bleistein et al.，2001)可知，在一維常速度背景下Green函數(shù)的解析表達式為：

(14)

假定接收點位于地表，即x=0，并將式(14)代入式(13)中，可以得到：

(15)

根據(jù)Fourier變換對的表達式：

(16)

可以求得速度擾動量的表達式：

(17)

由(17)式可知，通過Fourier變換可以快速由虛擬波場的數(shù)據(jù)求得界面的速度擾動量.

根據(jù)圖1階躍函數(shù)以及圖2 Delta函數(shù)的頻譜圖可知，階躍函數(shù)的頻譜主要集中在低頻區(qū)域，而且隨著頻率的升高振幅逐漸減小，而Delta函數(shù)的頻譜振幅始終等于1.圖3、圖4分別展示了全頻帶、缺少零頻信息、10～100 Hz、10～40 Hz的階躍函數(shù)以及Delta函數(shù)，黑色實線表示函數(shù)的真實值，黑色圓圈表示有限帶寬的階躍函數(shù)和Delta函數(shù).隨著帶寬越來越窄，階躍函數(shù)變形越來越嚴重，而Delta函數(shù)幾乎不受頻帶的影響，因此可知有限帶寬的Delta函數(shù)比有限帶寬的階躍函數(shù)更容易識別電性界面.

圖1 階躍函數(shù)及其頻譜圖(a)階躍函數(shù)；(b)階躍函數(shù)的頻譜圖.Fig.1 Step function and its spectrum(a)Step function;(b)Spectrum diagram of step function.

圖2 Delta函數(shù)及其頻譜圖(a)Delta函數(shù)；(b)Delta函數(shù)的頻譜圖.Fig.2 Delta function and its spectrum(a)Delta function;(b)Spectrum diagram of Delta function.

圖3 不同頻率的階躍函數(shù)(a)全頻帶；(b)缺少0頻信息；(c)10～100 Hz；(d)10～40 Hz.Fig.3 Step functions of different frequencies(a)Full frequency band;(b)Lack of 0 frequency information;(c)10 to 100 Hz;(d)10 to 40 Hz.

圖4 不同頻率的Delta函數(shù)(a)全頻帶；(b)缺少0頻信息；(c)10～100 Hz；(d)10～40 Hz.Fig.4 Delta functions of different frequencies(a)Full frequency band;(b)Lack of 0 frequency information;(c)10 to 100Hz;(d)10 to 40 Hz.

由于瞬變電磁擬地震成像僅需要定位電性界面的位置，因此，本文使用Delta函數(shù)形式的反射率函數(shù)識別電性界面，記為β(x)，通過對速度擾動量函數(shù)α(x)求一階導(dǎo)數(shù)，即可得到界面反射率函數(shù)的表達式：

(18)

1.2 二階Born近似算法

已知Born近似算法忽略散射場，雖然使積分方程線性化，但是由于其只利用了散射序列中的第一項，丟失了部分有效信息，只適用于小擾動量介質(zhì)的成像(丁科和宋守根，2004).因此，本文將對散射序列中的各項進行分析研究.

已知波動方程(4)的解具有如下的形式：

×u(r,rs,k)d3r,(19)

=u0+u1+u2+…

(20)

式(20)稱為散射序列或Born序列.已知u(r,rs,ω)=ui(r,rs,ω)+us(r,rs,ω)，因此散射場的表達式為：

=u1+u2+…

(21)

根據(jù)(21)式可知，去掉第二項以后，即為Born近似的表達式，同時在(21)式中存在如下的遞推關(guān)系：

(22)

對式(21)取前兩項，可得：

×G(r′,r″,k)u(r″,rs,k)d3r″d3r′.

(23)

當(dāng)偏移距為零時，即xs=xg=0，則Green函數(shù)的表達式為：

(24)

將式(24)代入式(23)中，可以得到：

×u(r″,0,k)dx″dx′.

(25)

已知Green函數(shù)的另一種表達形式為：

(26)

將式(26)代入式(25)，整理可得下面的方程組：

(27)

由式(27)可知，各變量之間存在循環(huán)迭代的關(guān)系.令us(0,0,k)=εus(0,0,k)，則式(27)可寫為：

(28)

當(dāng)n=1時，可得：

(29)

當(dāng)n=2時，可得：

(30)

當(dāng)n=3時，可得：

(31)

通過循環(huán)迭代依次可以求得W(k)、A(kx,k)和T(kx,0,k)，從而得到二階Born近似反射率函數(shù)β(x)的表達式：

(32)

2 算法驗證

2.1 電阻率變化對結(jié)果的影響

本節(jié)主要探究電阻率變化對Born近似算法以及二階Born近似算法結(jié)果的影響，設(shè)置了H、K型模型，電性界面距地面分別為100 m、200 m，保持第一層與第三層的電阻率值不變，使模型第二層的電阻率依次變化.H型模型第一層與第三層的電阻率均設(shè)置為100 Ωm，設(shè)置第二層電阻率與背景電阻率的差異越來越大，依次為90 Ωm、80 Ωm、70 Ωm、60 Ωm；K型模型第一層與第三層的電阻率均設(shè)置為10 Ωm，第二層電阻率依次設(shè)置為20 Ωm、30 Ωm、40 Ωm、50 Ωm.模型示意圖如圖5所示，均選擇第一層作為背景速度v0，分別計算了Born近似算法和二階Born近似算法的結(jié)果，如圖6，7所示.

圖5 三層模型示意圖(a)H型模型；(b)K型模型.Fig.5 Three-layer model diagram(a)H-type model;(b)K-type model.

圖6 H型模型Born近似、二階Born近似以及真實值結(jié)果圖(a)ρ1=100 Ωm，ρ2=90 Ωm，ρ3=100 Ωm；(b)ρ1=100 Ωm，ρ2=80 Ωm，ρ3=100 Ωm；(c)ρ1=100 Ωm，ρ2=70 Ωm，ρ3=100 Ωm；(d)ρ1=100 Ωm，ρ2=60 Ωm，ρ3=100 Ωm.Fig.6 Born approximation,second-order Born approximation and real value of H-type model

圖中，黑色實線、黑色加號、黑色圓圈分別表示三層模型電性界面反射率函數(shù)的真實值、Born近似值、二階Born近似值，其幅值表示電性界面反射率函數(shù)的大小.從H型模型的結(jié)果圖中可知，黑色加號與黑色實線僅在第一個電性界面處重合，隨著中間層與背景電阻率的差異變大，兩者所指示的第二個電性界面的位置差也越來越大，K型模型結(jié)果類似.說明了Born近似算法僅能正確定位第一個電性界面的位置，隨著電阻率差異變大，對第二個電性界面的定位偏差也會變大，這是由于在第一個電性界面處，第一層真實的速度與背景速度是一致的，因此對于第一個電性界面能夠正確反映；對于第二個電性界面，由于第二層的速度與背景速度是不一致的，因此會造成一定的誤差.反觀黑色圓圈與黑色實線在兩個電性界面處幾乎重合，說明了二階Born近似算法可以較好地反映地下電性界面的位置，其結(jié)果與真實值的偏差幾乎不會隨著電阻率差異變大而變大.以上分析說明了Born近似算法僅適用于小擾動量的成像問題，對于大擾動量模型，幾乎不能正確定位界面的位置；而二階Born近似算法幾乎不受電阻率變化的影響，在大擾動量模型中具有一定的優(yōu)越性.

圖7 K型模型Born近似、二階Born近似以及真實值結(jié)果圖(a)ρ1=10 Ωm，ρ2=20 Ωm，ρ3=10 Ωm；(b)ρ1=10 Ωm，ρ2=30 Ωm，ρ3=10 Ωm；(c)ρ1=10 Ωm，ρ2=40 Ωm，ρ3=10 Ωm；(d)ρ1=10 Ωm，ρ2=50 Ωm，ρ3=10 Ωm.Fig.7 Born approximation,second-order Born approximation and real value of K-type model

2.2 層厚變化對結(jié)果的影響

本節(jié)主要探究模型層厚變化對結(jié)果的影響，設(shè)置了HK型模型，模型的四層電阻率分別設(shè)置為：100 Ωm、10 Ωm、100 Ωm、1 Ωm，模型示意圖如圖8所示，分別設(shè)置兩組模型：第一組模型第一層、第三層厚度均設(shè)置為100 m，第二層厚度依次設(shè)置為40 m、30 m、20 m、10 m；第二組模型第一層厚度為100 m，第二層厚度為50 m，第三層層厚依次設(shè)置為40 m、30 m、20 m、10 m.在計算過程中，均將第一層作為背景速度，相關(guān)計算結(jié)果如圖9，10所示.

圖8 四層HK型模型示意圖(a)改變第二層厚度；(b)改變第三層厚度.Fig.8 The schematic diagram of the four-layer HK model(a)Changing the thickness of the second layer;(b)Changing the thickness of the third layer.

圖9 改變第二層厚度的HK模型結(jié)果圖(a)h1=100 m，h2=40 m，h3=100 m；(b)h1=100 m，h2=30 m，h3=100 m；(c)h1=100 m，h2=20 m，h3=100 m；(d)h1=100 m，h2=10 m，h3=100 m.Fig.9 Results of HK model with changing the second layer′s thickness

圖中，黑色實線、黑色加號、黑色圓圈分別表示地下電性界面反射率函數(shù)的真實值、Born近似值、二階Born近似值.根據(jù)兩組模型的計算結(jié)果可知，各層電阻率不變、只改變第二層或者第三層厚度的情況下，黑色加號與黑色實線僅在第一個電性界面處重合，對于第二個、第三個電性界面，兩者的偏差較大，說明Born近似算法在厚度不斷變化的模型中除了第一個電性界面外，其余界面均不能正確定位；然而黑色圓圈與黑色實線幾乎完全重合，不受模型層厚變化的影響，均能正確定位電性界面的具體位置，這進一步說明二階Born近似算法相對比Born近似算法具有一定的優(yōu)越性.

3 理論模型計算

3.1 二維模型

圖10 改變第三層厚度的HK模型結(jié)果圖(a)h1=100 m，h2=50 m，h3=40 m；(b)h1=100 m，h2=50 m，h3=30 m；(c)h1=100 m，h2=50 m，h3=20 m；(d)h1=100 m，h2=50 m，h3=10 m.Fig.10 Results of HK model with changing the third layer′s thickness

圖11 二維模型示意圖(a)XOZ平面；(b)XOY平面.Fig.11 Diagram of two-dimensional model(a)XOZ plane;(b)XOY plane.

圖12 二維模型的多測道圖Fig.12 Multi-channel diagram of two-dimensional model

圖13 二維模型的虛擬波場圖Fig.13 Pseudo wave-field diagram of two-dimensional model

在成像的計算過程中將第一層作為背景速度，圖14為二維模型的Born近似和二階Born近似算法的成像結(jié)果圖，虛線表示二維模型真實的地質(zhì)界面位置，從圖中可知兩種算法均能反映出模型有三個電性界面.Born近似的結(jié)果僅能準確地定位第一個地質(zhì)界面，上盤距地面大致100 m，下盤距地面大致150 m，第二個、第三個電性界面的成像結(jié)果與實際地質(zhì)界面的位置分別大致相差30 m和40 m，而且斷層的形態(tài)與真實模型有一定差距；二階Born近似算法成像的結(jié)果對地下三個電性界面的位置定位準確，成像的結(jié)果與虛線幾乎重合，橫向上在X=-50～50 m處斷層的形態(tài)較為明顯，各界面的埋深與模型設(shè)置幾乎是一致的.可見，相比于Born近似算法的成像結(jié)果，二階Born近似算法的成像效果較好.

圖14 二維模型的成像結(jié)果(a)Born近似結(jié)果；(b)二階Born近似結(jié)果.Fig.14 Imaging results of the two-dimensional model(a)Born approximation results;(b)Second-order Born approximation results.

3.2 兩個低阻薄板模型

圖15 兩個低阻薄板模型示意圖(a)XOZ平面；(b)XOY平面.Fig.15 Schematic diagram of two low-resistance thin plates(a)XOZ plane;(b)XOY plane.

圖16 三維模型的瞬變電磁數(shù)據(jù)(a)Y=40 m；(b)Y=0 m；(c)Y=-60 m.Fig.16 Transient electromagnetic data of three-dimensional model

圖17 三維模型的虛擬波場圖(a)Y=40 m；(b)Y=0 m；(c)Y=-60 m.Fig.17 Pseudo wave-field diagram of three-dimensional model

圖19 測線Y=0 m的成像結(jié)果(a)Born近似結(jié)果；(b)二階Born近似結(jié)果.Fig.19 Imaging results of line Y=0 m(a)Born approximation results;(b)Second-order Born approximation results.

圖20 測線Y=-60 m的成像結(jié)果(a)Born近似結(jié)果；(b)二階Born近似結(jié)果.Fig.20 Imaging results of line Y=-60 m(a)Born approximation results;(b)Second-order Born approximation results.

在成像的計算過程中將背景電阻率作為背景速度，圖18—20分別展示了測線Y=40 m、Y=0 m、Y=-60 m的Born近似算法以及二階Born近似算法的擬地震成像結(jié)果圖，其中黑色虛線表示該測線下方異常體的真實位置.測線Y=40 m處，Born近似算法與二階Born近似算法均對上部異常體的上下界面以及下部異常體的上界面有所反映，但是Born近似算法的成像結(jié)果與真實界面的位置大約相差30～40 m，而且界面傾斜與真實模型形態(tài)相差較大；而二階Born近似算法的結(jié)果能夠準確定位界面的位置以及異常體的形態(tài).測線Y=0 m的Born近似成像的結(jié)果不能準確定位界面的位置，尤其是下部異常體的界面與真實界面大約相差50 m；二階Born近似擬地震的結(jié)果可以較好地定位兩個異常體的界面位置，而且其形態(tài)與真實模型相一致.測線Y=-60 m的擬地震B(yǎng)orn近似成像結(jié)果僅反映了上部異常體的部分界面信息以及下部異常體的界面位置，但結(jié)果與真實模型的界面位置存在一定的偏差，而二階Born近似算法的結(jié)果對界面的位置以及異常體的形態(tài)均能夠準確定位.測線Y=40 m雖然只穿過上部異常體、未穿過下部異常體，但是計算瞬變電磁擴散場時會受到下部異常體的影響，因此成像的結(jié)果中含有下部異常體的部分信息，同理可知測線Y=-60 m處雖只穿過下部異常體，但對上部異常體也有所反映.通過對三維低阻薄板模型不同位置的3條測線進行Born近似算法以及二階Born近似算法成像的結(jié)果可知，二階Born近似算法的結(jié)果既可以準確反映兩個低阻薄板上下界面的位置，又可以精確反映界面的形態(tài)，因此該方法相對比Born近似算法具有較高的精度和界面識別能力.

4 結(jié)論

本文基于瞬變電磁波場反變換，在波動方程中將總波場、總速度分為關(guān)于背景和擾動量的項，然后運用Green定理求解波動方程，并引入Born近似算法將積分方程線性化，建立了虛擬波場與界面反射率函數(shù)之間的關(guān)系.通過分析散射序列以及遞推關(guān)系，本文保留了散射序列的前兩項，推導(dǎo)出了二階Born近似算法的表達式.通過探究電阻率變化、層厚變化對Born近似算法以及二階Born近似算法的影響，可知二階Born近似算法幾乎不受影響，均能準確定位界面位置；根據(jù)理論模型的成像結(jié)果可知，相對比Born近似算法，二階Born算法既可以準確反映二維模型地質(zhì)界面的位置和形態(tài)，又能定位三維模型中的兩個低阻薄板的埋深以及相對位置.綜上可知，在波場反變換的基礎(chǔ)上，二階Born近似算法可以實現(xiàn)瞬變電磁擬地震快速、準確成像的目的.

本文主要研究散射序列的前兩項即二階Born近似與Born近似算法的結(jié)果對比，在后面的研究中將對散射序列的更高階項進行探究.

致謝我們特別向?qū)徃迦撕途庉媽Ρ疚奶岢龅慕ㄔO(shè)性意見表示感謝．