劉新彤,劉四新, 孟 旭,傅 磊

1.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026 2.澳門科技大學(xué)太空科學(xué)研究所，澳門 999078

0 引言

跨孔探地雷達層析成像技術(shù)不斷發(fā)展，已經(jīng)在水文、環(huán)境和工程物探中得到了驗證與有效應(yīng)用。在數(shù)據(jù)采集過程中，跨孔雷達將兩個天線放置在相隔一定距離的兩個井孔中，一個作為發(fā)射源，另一個作為接收器。發(fā)射天線發(fā)出的雷達波穿過井孔之間的介質(zhì)被接收天線采集，該雷達波包含井孔之間介質(zhì)的物性參數(shù)(如介電常數(shù)和電導(dǎo)率)空間分布信息。傳統(tǒng)的走時層析成像和衰減層析成像技術(shù)雖然可以分別有效地反演出地下介質(zhì)的介電常數(shù)和電導(dǎo)率分布情況，但由于上述方法基于射線理論，僅僅使用了雷達波的部分信息，因而所得到的結(jié)果分辨率有限?？缈桌走_全波形反演技術(shù)在反演過程中使用雷達波的全部信息，所得到的地下介質(zhì)物性參數(shù)分布結(jié)果的分辨率可以達到亞波長級別。

近些年來，國際地球物理探測學(xué)界對全波形反演(full waveform inversion， FWI)方法持續(xù)關(guān)注，波形反演技術(shù)在很多領(lǐng)域都取得了顯著進展，尤其是在地震方面，已經(jīng)得到了相對成熟的理論架構(gòu)。早期，Tarantola[1]將廣義最小二乘反演理論應(yīng)用到了全波形反演方法中，通過合成數(shù)據(jù)與觀測數(shù)據(jù)計算出殘差并進行反傳獲得反傳波場，利用反傳波場與正傳波場的互相關(guān)來構(gòu)建梯度，形成了一套完整的時間域波形反演方法理論架構(gòu)。隨后，Pratt等[2-4]將全波形反演方法擴展到了頻率域，提出頻率域全波形反演方法；馮晅等[5]提出基于逐減隨機震源采樣法的頻率域波形反演方法。與地震全波形反演方法相比，探地雷達全波形反演方法的研究起步較晚，且目前仍然處于發(fā)展初期。Kuroda等[6]采用全波形反演方法來解釋跨孔雷達數(shù)據(jù)，獲得了介電常數(shù)成像結(jié)果。Ernst[7]等使用全波形反演方法處理跨孔雷達數(shù)據(jù)，采用級聯(lián)更新的方式同時反演了介電常數(shù)和電導(dǎo)率，并取得了較好的效果。Meles等[8]利用矢量全波形反演技術(shù)處理探地雷達數(shù)據(jù)，進行了同時迭代介電常數(shù)和電導(dǎo)率的全波形反演。在國內(nèi)，波形反演也取得了一定的進展。吳俊軍等[9]利用攝動法求解全波形反演的目標(biāo)函數(shù)并詳細推導(dǎo)了基礎(chǔ)理論，采用最速下降法在迭代過程中同時更新介電常數(shù)和電導(dǎo)率；孟旭等[10]進一步將跨孔雷達全波形反演方法發(fā)展到了頻率域，推導(dǎo)了基于對數(shù)目標(biāo)函數(shù)的頻率域全波形反演的基礎(chǔ)理論，使用共軛下降法同時反演介電常數(shù)和電導(dǎo)率。在處理實際數(shù)據(jù)方面，Ernst等[11]提出一種二維全波形反演策略，采用級聯(lián)方式反演兩種物性參數(shù)，反演跨孔雷達實際數(shù)據(jù)。但是，到目前為止，由于實際數(shù)據(jù)信噪比較低，國內(nèi)外使用全波形反演方法處理探地雷達實際數(shù)據(jù)并沒有得到良好的應(yīng)用。

FWI方法廣泛采用局部優(yōu)化的方式進行反演，但由于其本身對初始模型高度依賴，面對先驗信息缺失的地下介質(zhì)時，無法從準(zhǔn)確的初始模型出發(fā)；且在實際采集到的數(shù)據(jù)中包含噪聲以及存在低頻信息缺失的情況下，局部優(yōu)化方法很容易陷入局部極值[12]。在目前的波形反演研究領(lǐng)域中，普遍選擇為其構(gòu)建較為準(zhǔn)確的初始模型，來降低非線性程度以獲取更好的結(jié)果。而對于構(gòu)建大尺度的初始模型，低頻信息的使用尤為關(guān)鍵。對于實際數(shù)據(jù)，由于帶寬或數(shù)據(jù)采集濾波器設(shè)計等原因，低頻信息往往難以采集；雖然在實驗室條件下能夠采集到一定量的低頻數(shù)據(jù)，但是難度極高，成本巨大，所采集到的數(shù)據(jù)并不完整且應(yīng)用難度較大，無法推廣到探地雷達野外數(shù)據(jù)的處理中。以中心頻率為100 MHz的鉆孔雷達為例，難以獲得30 MHz以下的有效信息。因此，重新構(gòu)建探地雷達數(shù)據(jù)中的有效低頻信息，成為了利用全波形反演方法處理跨孔雷達數(shù)據(jù)的重點。

本文在時間域跨孔雷達波形反演的基礎(chǔ)上，使用包絡(luò)目標(biāo)函數(shù)，采用求導(dǎo)法，詳細推導(dǎo)了包絡(luò)波形反演(EWI)方法的基礎(chǔ)理論，將EWI成功推廣到了探地雷達(ground penetrating radar，GPR)領(lǐng)域，并采用共軛梯度法同時反演了兩種物性參數(shù)(介電常數(shù)和電導(dǎo)率)；通過模擬原始數(shù)據(jù)低頻信息缺失的情況驗證該方法在低頻缺失情況下的反演能力，反演結(jié)果在為地下介質(zhì)做出評價的同時，可以作為初始模型，降低非線性程度，進一步參與到傳統(tǒng)波形反演當(dāng)中， 來獲取分辨率更高的結(jié)果。該全波反演方法的正演部分采用高階時間域有限差分算法，使用了CPML(convolutional perfect matched layer)吸收邊界，并在每一次迭代過程中，對梯度進行了中值濾波等優(yōu)化方式。由于反演過程中數(shù)據(jù)量巨大，采用了Matlab的分布式計算引擎(matlab distributed computing engine， MDCE)并行運算的方式，以提高在普通微機上的運行效率。

1 正演問題表達形式

在討論探地雷達地球物理問題時，地下介質(zhì)的物性參數(shù)分布情況使用介電常數(shù)和電導(dǎo)率來描述，并假設(shè)磁導(dǎo)率是均勻不變的。在時間域，對于任意時間點t和空間x，Maxwell方程如下[7]：

(1)

式中：ε(x)和σ(x)分別表示介電常數(shù)和電導(dǎo)率的空間分布；s表示某一個特定的源；μ0表示磁導(dǎo)率，這里假定為常數(shù)；Es和Hs分別表示由電流密度源矢量Js產(chǎn)生的電場和磁場。式(1)可簡寫成如下形式：

(2)

式中，M(ε,σ)代表Maxwell方程的一個線性算子。忽略磁場項：

MEs=Js或Es=M-1Js。

(3)

在時間域，空間任意位置的電場可以用源與對應(yīng)位置格林函數(shù)之間的卷積來表示：

Es=g*Js。

(4)

式中，g為M算子對應(yīng)的格林函數(shù)。顯然，M-1包含離散近似的格林函數(shù)[2-4]。因此，

M-1=g。

(5)

2 目標(biāo)函數(shù)及梯度計算

2.1 傳統(tǒng)全波形反演

傳統(tǒng)全波形反演的主要目的是尋找目標(biāo)函數(shù)最小時對應(yīng)的介電常數(shù)ε與電導(dǎo)率σ的空間分布。傳統(tǒng)目標(biāo)函數(shù)表達式為

(6)

式中：d表示接收器；τ表示觀測時間；E(ε,σ)和Eobs分別為電場正演數(shù)據(jù)和實際數(shù)據(jù)。式(6)表示源s在接收點d處、觀測時間為τ時全部誤差的總和。目標(biāo)函數(shù)對應(yīng)的梯度由正傳合成波場與反傳伴隨波場的互相關(guān)獲得：

(7)

式中：p代表模型參數(shù)ε和σ；v為虛擬源向量，對于不同的模型參數(shù)具有不同的表達式；r(ξ)=E(ξ)-Eobs(ξ)，代表反傳波場的后向殘場源。

2.2 包絡(luò)波形反演

本文通過將包絡(luò)目標(biāo)函數(shù)引入到Meles等[7]的矢量波形反演技術(shù)當(dāng)中，重新構(gòu)建包絡(luò)全波形反演基礎(chǔ)理論。不同于Meles等使用擾動目標(biāo)函數(shù)一階近似來推導(dǎo)梯度的方法，本文使用對包絡(luò)目標(biāo)函數(shù)求導(dǎo)的方式來推導(dǎo)梯度。式(6)為反演中最廣泛使用的傳統(tǒng)目標(biāo)函數(shù)，本文給出基于該目標(biāo)函數(shù)基礎(chǔ)上的包絡(luò)目標(biāo)函數(shù)，即合成包絡(luò)波場與實際包絡(luò)波場差的2-范數(shù)：

(8)

(9)

進一步推導(dǎo)獲得

(10)

式中，對合成波場的偏微分可以通過使用模型參數(shù)p對公式(3)兩端求偏導(dǎo)獲得：

(11)

其中，

(12)

在時間域，空間任意位置的電場可以認(rèn)為是源與對應(yīng)位置格林函數(shù)之間的卷積：

(13)

其中，

(14)

將式(5)代入式(13)并寫成積分的形式：

(15)

(16)

(17)

為了實現(xiàn)最大程度的收斂，本文使用共軛梯度法來更新模型：

[ε(x)k+1]=[ε(x)k]-ζε,k·[Cε(x)k]，

(18)

[σ(x)k+1]=[σ(x)k]-ζσ,k·[Cσ(x)k]。

(19)

式中：ζε,k和ζσ,k分別為第k次迭代中介電常數(shù)和電導(dǎo)率的迭代步長；Cε(x)k和Cσ(x)k分別為介電常數(shù)和電導(dǎo)率對應(yīng)的共軛梯度方向，由當(dāng)前和上一次的梯度計算得到。

Cε(x)k=Sε(x)k+

(20)

Cσ(x)k=Sσ(x)k+

(21)

當(dāng)k=1時，滿足Cε(x)1=Sε(x)1和Cσ(x)1=Sσ(x)1。通過沿各自Cε、Cσ方向?qū)ふ覙O值點，我們能夠在同一次迭代中實現(xiàn)介電常數(shù)和電導(dǎo)率的同步反演。迭代步長由以下公式獲得：

[Eenv(ε+κεCε,k,σ)-Eenv(ε,σ)]，

(22)

[Eenv(ε,σ+κσCσ,k)-Eenv(ε,σ)]。

(23)

式中，κε、κσ為不同的小穩(wěn)定因子。在反演過程中必須小心地為其選擇適當(dāng)?shù)闹挡⑶以撝惦S著迭代不斷更新。

3 數(shù)值模擬

與傳統(tǒng)的FWI方法對比完整波形不同，我們的方法比較的是合成數(shù)據(jù)的包絡(luò)與觀測數(shù)據(jù)的包絡(luò)。這里我們以Ricker子波為例來分析包絡(luò)轉(zhuǎn)換的作用。圖1a顯示了一道Ricker子波的波形與它對應(yīng)的包絡(luò),通過對比可以看出，Ricker子波信號本身變化劇烈，而包絡(luò)波形變化相對緩慢。圖1b為子波頻譜與包絡(luò)頻譜的對比圖，可以很明顯看出存在低頻信息通過非線性轉(zhuǎn)換被加入到了包絡(luò)頻譜中。

傳統(tǒng)的跨孔雷達FWI在雷達數(shù)據(jù)低頻信息充足且初始模型足夠準(zhǔn)確的情況下可以獲得良好的反演結(jié)果；但是當(dāng)?shù)皖l信息缺少時，F(xiàn)WI方法往往無法恢復(fù)大尺度的背景信息，甚至錯誤地更新模型并導(dǎo)致反演最終無法收斂，非常容易陷入局部最小。雷達數(shù)據(jù)缺少低頻信息這種狀況在很多雷達實例中存在，面對類似情況的實例，F(xiàn)WI方法難以發(fā)揮作用；但EWI方法有能力在數(shù)據(jù)缺少大量低頻信息的情況下更新模型的長波長信息，其原因是它的伴隨波場含有低頻信息。我們?nèi)コ齊icker子波30 MHz以下的低頻成分，圖2a中為去除低頻成分后的Ricker子波及其包絡(luò)，可以看出去除低頻成分后，Ricker子波信號變動更劇烈，而包絡(luò)曲線依然平穩(wěn)。圖2b為其頻譜對比圖，與圖1b相比較，可以非常明顯地看出，在除去低頻成分后，Ricker子波頻譜低頻部分雖然被移除，但其包絡(luò)低頻部分幾乎不受影響，仍然有明顯的低頻信息被添加。

為了驗證EWI方法的成像能力，我們建立了如圖3所示的復(fù)雜小模型，模型尺寸為6 m×6 m，包含3個地層：第一層與第三層物性參數(shù)相同，相對介電常數(shù)均為5，電導(dǎo)率為0.001 S·m-1；中間層的相對介電常數(shù)為5.5，電導(dǎo)率為0.002 8 S·m-1，并埋藏兩個柱狀異常體，位于深度3 m處，間隔1 m，相對介電常數(shù)為7，電導(dǎo)率為0.008 S·m-1。在本次模擬中共13個發(fā)射源，每組發(fā)射對應(yīng)13個接收器。發(fā)射器與接收器位置均為深度0～6 m，并以0.5 m等間隔排列；天線中心頻率為100 MHz。

首先，我們對比了兩種方法基于均勻初始模型反演的結(jié)果。均勻初始模型的相對介電常數(shù)為5.5，電導(dǎo)率為0.002 8 S·m-1。初始模型物性參數(shù)與原始模型中間層背景物性參數(shù)一致。由于原始數(shù)據(jù)包含了完整且正確的長波長信息，所以FWI和EWI兩種反演方法均能夠獲得良好的反演結(jié)果。其中傳統(tǒng)的全波形反演方法所得到的結(jié)果分辨率略高于基于包絡(luò)的波形反演結(jié)果，這是因為FWI方法使用的信息更加完整，但也有更高的非線性程度。從圖4反演結(jié)果來看，F(xiàn)WI方法在波形信息完整時，反演結(jié)果要優(yōu)于EWI方法，尤其對于中間兩個管線異常體，傳統(tǒng)波形反演電導(dǎo)率結(jié)果優(yōu)勢明顯。

圖1 Ricker子波和包絡(luò)波形(a)及其頻譜(b)Fig.1 Ricker wavelet and its envelope (a) and spectrum (b)

圖2 去除30 MHz以下低頻成分Ricker子波及其包絡(luò)(a)和頻譜(b)Fig.2 Ricker wavelet without low frequency data under 30 MHz and its envelope (a) and spectrum (b)

a.相對介電常數(shù)；b.電導(dǎo)率。白色圓圈表示發(fā)射源，叉號表示接收器，下同。圖3 原始模型Fig.3 Original model

a、b. FWI；c、d. EWI。圖4 完整頻率信息波形反演結(jié)果Fig.4 Waveform inversion results of complete frequency information

接著，我們提取了原始數(shù)據(jù)30 MHz以下的低頻信息，以驗證FWI和EWI兩種方法對低頻信息的反演能力。如圖5所示，F(xiàn)WI方法在單獨反演低頻信息時，反演能力明顯下降，圖5a、b是迭代20次的結(jié)果，雖然能夠反演出一定的地層和異常體信息，但已經(jīng)遠遠不及EWI方法(圖5c、d)；并沒有選擇更多的迭代次數(shù)是因為由于低頻缺失，F(xiàn)WI方法容易陷入局部最小，隨著迭代的進行，結(jié)果將會朝著錯誤的方向發(fā)展，最終無法收斂，完全無法準(zhǔn)確反映地下信息。而包絡(luò)波形反演能夠獲得良好的結(jié)果，目標(biāo)函數(shù)收斂，且最終結(jié)果(圖5c、d)與圖4中使用全部波形信息反演的結(jié)果相差不大。該組模擬證明，EWI方法對低頻信息的利用能力高于FWI方法。對比提取低頻信息的反演結(jié)果與完整波形信息的反演結(jié)果可以看出，無論使用FWI還是EWI方法，由于高頻信息的缺失，反演對于中間細節(jié)的刻畫明顯下降，極易陷入局部最小。

最后，我們模擬現(xiàn)實中低頻缺失的情況，以了解當(dāng)?shù)皖l缺失、FWI方法無法獲得理想反演結(jié)果的情況下，EWI方法是否仍然可以獲得良好結(jié)果，使跨孔雷達波形反演體系更加完善且更符合實際應(yīng)用。圖6中，由于相同的原因，F(xiàn)WI結(jié)果與圖5一致，均為20次迭代的結(jié)果。由圖6可見：在低頻信息缺失的情況下，無論是介電常數(shù)還是電導(dǎo)率，F(xiàn)WI均無法獲得與利用完整波形信息進行FWI時一致的結(jié)果，目標(biāo)函數(shù)不收斂；與FWI相對比，EWI的結(jié)果更好，且目標(biāo)函數(shù)收斂，能夠為地下地層與異常體從位置和物性參數(shù)方面提供定量解釋。

a、b. FWI；c、d. EWI。圖5 30 MHz以下低頻信息反演結(jié)果Fig.5 Inversion result of low frequency data under 30 MHz

a、b. FWI；c、d. EWI。圖6 去除30 MHz以下低頻信息反演結(jié)果Fig.6 Inversion result without low frequency under 30 MHz

4 結(jié)論

1)與全波形反演方法相比，包絡(luò)波形反演方法是一種對頻率成分更加穩(wěn)定且非線性更低的波形反演方法?？缈桌走_數(shù)值模擬證明該方法在跨孔雷達領(lǐng)域有著很高的應(yīng)用價值，尤其對于低頻信息缺失的數(shù)據(jù)，包絡(luò)波形反演可以在傳統(tǒng)波形反演不適用的情況下仍然提供良好的反演結(jié)果。

2)雖然對于完整的波形信息，全波形反演方法結(jié)果要優(yōu)于包絡(luò)波形反演的結(jié)果，尤其在電導(dǎo)率的結(jié)果上全波形反演方法優(yōu)勢明顯；但在低頻信息缺失時，包絡(luò)波形反演無論是介電常數(shù)還是電導(dǎo)率，結(jié)果明顯優(yōu)于全波形反演方法，可以為地下物性參數(shù)提供解釋。

3)包絡(luò)波形反演方法是對跨孔雷達波形反演的一種補充方法，使得跨孔雷達波形反演符合更多應(yīng)用實例。

[1]Tarantola A. Inversion of Seismic-Reflection Data in the Acoustic Approximation [J]. Geophysics, 1984, 49: 1259- 1266.

[2] Pratt R G, Shin C, Hicks G J. Gauss-Newton and Full Newton Methods in Frequency Domain Seismic Waveform Inversion [J]. Geophysical Journal International, 1998, 133: 341-362.

[3] Pratt R G. Seismic Waveform Inversion in the Fre-quency Domai: Part 1: Theory and Verification in a Physical Scale Model [J]. Geophysics, 1999, 64: 888-901.

[4] Pratt R G, Shipp R M. Seismic Waveform Inversion in the Frequency Domain: Part 2: Fault Delineation in Sediments Using Crosshole Data [J]. Geophysics, 1999, 64: 902-914.

[5] 馮晅, 魯曉滿, 劉財,等. 基于逐減隨機震源采樣法的頻率域二維黏滯聲波方程全波形反演[J]. 吉林大學(xué)學(xué)報(地球科學(xué)版), 2016, 46(6):1865-1873.

Feng Xuan, Lu Xiaoman, Liu Cai, et al. Frequency-Domain Full Waveform Inversion of 2D Viscous Acoustic Wave Equation Using Decreasing Random Shot Subsampling Method[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(6): 1865- 1873.

[6] Kuroda S, Takeuchi M, Kim H J. Full Waveform In-version Algorithm for Interpreting Cross-Borehole GPR Data [J]. Seg Technical Program Expanded Abstracts, 1949, 24(1): 2668.

[7] Ernst J R, Maurer H, Green A G, et al. Full-Wave-form Inversion of Crosshole Radar Data Based on 2-D Finite-Difference Time-Domain Solutions of Maxwell’s Equations [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(9): 2807-2828.

[8] Meles G A, Kruk J V D, Greenhalgh S A, et al. A New Vector Waveform Inversion Algorithm for Simultaneous Updating of Conductivity and Permittivity Parameters from Combination Crosshole/Borehole-to-Surface GPR Data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(9): 3391-3407.

[9] 吳俊軍,劉四新,李彥鵬,等．跨孔雷達全波形反演成像方法的研究[J]. 地球物理學(xué)報, 2014, 57(5):1623-1635.

Wu Junjun, Liu Sixin, Li Yanpeng, et al. Study of Cross-Hole Radar Tomography Using Full-Waveform Inversion [J]. Chinese Journal of Geophysics, 2014, 57(5): 1623-1635.

[10] 孟旭, 劉四新, 傅磊,等. 基于對數(shù)目標(biāo)函數(shù)的跨孔雷達頻域波形反演[J]. 地球物理學(xué)報,2016,59(5):1875-1887.

Meng Xu, Liu Sixin, Fu Lei, et al. Frequency Domain Waveform Inversion of Cross-Hole GPR Data Based on a Logarithmic Objective Function [J]. Chinese Journal of Geophysics, 2016, 59(5):1875-1887.

[11] Ernst J R, Green A G, Maurer H, et al. Application of a New 2D Time-Domain Full-Waveform Inversion Scheme to Crosshole Radar Data [J]. Geophysics, 2007, 72(5): J53-J64.

[12] Virieux J, Operto S. An Overview of Full-Waveform Inversion in Exploration Geophysics [J]. Geophysics, 2009, 74 (6), WCC1-WCC26.

[14]Chi B X, Dong L G, Liu Y Z. Full Waveform Inversion Based on Envelope Objective Function [C]// 75th EAGE Conference & Exhibition Incorporating Spe Europec 2013. London: EAGE, 2013.

[15]Chi B X, Dong L G, Liu Y Z. Full Waveform Inversion Method Using Envelope Objective Function Without Low Frequency Data [J]. Journal of Applied Geophysics, 2014, 109: 36-46.

[16]Wu R S, Luo J R, Wu B Y. Seismic Envelope Inversion and Modulation Signal Model [J]. Geophysics, 2014, 79(3): WA13-WA21.