魯凱亮, 李貅*, 戚志鵬*, 樊亞楠, 周建美,李文翰, 李賀, 張明晶, 王揚州

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院， 西安 710054 2 長安大學(xué)地球物理場多參數(shù)綜合模擬實驗室(中國地球物理學(xué)會重點實驗室)， 西安 710054 3 山東大學(xué)齊魯交通學(xué)院， 濟南 250061 4 山東省交通規(guī)劃設(shè)計院， 濟南 250031 5 北京探創(chuàng)資源科技有限公司， 北京 100071

0 引言

目前瞬變電磁反演結(jié)果為連續(xù)電阻率的空間分布，在識別電性分界面等構(gòu)造時具有一定的困難.而瞬變電磁的擬地震成像是實現(xiàn)三維成像的有效手段之一(Xue et al.，2013；戚志鵬等，2013；鐘華森等，2016).在實現(xiàn)擬地震成像之前，需要把瞬變電磁擴散場轉(zhuǎn)換到虛擬波場.只有在此基礎(chǔ)上，才能實現(xiàn)穩(wěn)定高效的瞬變電磁法擬地震解釋.

瞬變電磁擴散場主要描述的是低頻電磁場傳播過程中的擴散和感應(yīng)特征，具有較強的體積效應(yīng)，因而對電性界面的分辨能力較差(Khan et al.，2018；底青云等, 2019; Di et al.，2020；Xue et al.，2019；Li et al.，2018)；而波動場能夠比較好地描述波場的一系列特征，所以基于波動方程的擬地震成像可以很好地刻畫地層的界面信息(Qi et al.，2015).而直接將瞬變電磁的擴散場引入至波場的處理流程中是不合適的，需要通過數(shù)學(xué)手段進行變換.Lee等(1989)給出了虛擬波場到擴散場的積分公式，實現(xiàn)了穩(wěn)定的正變換；隨后，Lee和Xie(1993)根據(jù)波場正變換的積分公式，嘗試進行波場反變換的計算，他們指出波場反變換需要求解第一類的Fredholm積分方程，屬于嚴重的不適定性問題.之后，Gershenson(1997)分別使用了直接求解法、反褶積法和拉氏正弦變換法計算虛擬波場，并對這三種方法的優(yōu)劣進行了比較.陳本池等(1999)使用正則化方法進行虛擬波場反變換的求解，并對層狀模型進行試算，結(jié)果表明從擴散場轉(zhuǎn)換至虛擬波場后，在電性分界面處會出現(xiàn)“反射”和“折射”等現(xiàn)象，這一發(fā)現(xiàn)為電磁場的擬地震處理解釋提供了重要的理論依據(jù).近些年，李貅等(2005)，朱宏偉等(2010)，張軍等(2011)將正則化算法應(yīng)用于波場反變換中，取得較好的結(jié)果；戚志鵬等(2013)通過使用預(yù)條件正則化共軛梯度法，實現(xiàn)了全時域的虛擬波場反變換算法，通過引入地震勘探中的脈沖反褶積方法，提高了虛擬波場的縱向分辨率.樊亞楠等(2019)使用掃時波場變換算法，進一步提高了波場反變換的穩(wěn)定性.雖然以上方法都取得了長足的發(fā)展，但是虛擬波場反變換的精度和效率還有待進一步提高.

在計算虛擬波場反變換時，將其離散化后的線性方程組具有高度的病態(tài)性.近年來，精細積分法作為求解高度病態(tài)線性方程組的方法，在各個領(lǐng)域得到了日益廣泛的應(yīng)用.該方法通過把求逆矩陣的過程轉(zhuǎn)化為一個求解無窮積分的過程，大大降低了求解病態(tài)線性方程組的難度.精細積分法(鐘萬勰，1994)在計算矩陣指數(shù)函數(shù)時，具有效率高、精度高等優(yōu)點，已經(jīng)在各類工程問題中取得了較大進展(Gu et al.，2001)，而目前將精細積分法引入到虛擬波場反變換的工作鮮有提及，本文將對這一方法的引入初做嘗試.

1 波場變換基本原理

1.1 基本原理

在導(dǎo)電介質(zhì)中，忽略位移電流，瞬變電磁場滿足擴散方程.為了不失一般性，取f(x,y,z,t)為瞬變電磁場的電磁分量.

根據(jù)文獻(Lee et al.，1989)可知，瞬變電磁擴散場轉(zhuǎn)換到虛擬波場的表達式為

(1)

式(1)為第一類Fredholm型積分方程，從擴散場求解虛擬波場是一個典型的不適定問題.對式(1)進行離散，可以得到：

(2)

其中：

(3)

離散后得到的線性方程組是高度病態(tài)的，且隨著階數(shù)的增加，矩陣的條件數(shù)急劇增大(戚志鵬等，2013)，本文采用精細積分法求解(1)式.

圖1 核函數(shù)曲線圖Fig.1 Kernel function graph

1.2 精細積分法

考慮病態(tài)線性方程組

Ax=b,

(4)

式中A為n階正定矩陣，b為n維實向量，記H=-A、r=b，則式(4)可以轉(zhuǎn)化為

Hx+r=0,

(5)

考慮如下一階微分方程

(6)

式(6)是瞬態(tài)熱傳導(dǎo)方程解的形式；穩(wěn)態(tài)的熱傳導(dǎo)方程的解具有式(5)的形式，而式(6)的解可以表示為

上式中，H是負定矩陣，r為常向量.從上式易知，當(dāng)時間t′→∞時，方程(6)趨向于方程(5)，且容易看出exp(Ht′)→0，所以式(7)中的積分項就逼近式(5)的解，即

(9)

所以，式(4)的解就可以寫為(富明慧，2018)

(10)

從上面的分析我們可以得出一個結(jié)論，病態(tài)線性方程組的求解可以轉(zhuǎn)化為一個求積分的穩(wěn)定過程.

(11)

1.3 迭代終止準(zhǔn)則

前文主要介紹了精細積分法，這一節(jié)給出迭代終止條件.理論上講，迭代次數(shù)越多，結(jié)果越精確，但是當(dāng)積分區(qū)間達到一定數(shù)值范圍后，隨著計算誤差的累積以及矩陣的高病態(tài)性，反而會導(dǎo)致精度隨積分區(qū)間的增加迅速下降，所以有必要在虛擬波場計算的過程中添加一個迭代終止條件.

目前有如下四種迭代終止條件：(1)對迭代次數(shù)設(shè)置一個迭代上限m；(2)終止條件設(shè)置為errk<ε，ε為給定的誤差容限，errk=‖xk+1-xk‖；(3)當(dāng)?shù)鷼埐钣袉握{(diào)下降和單調(diào)上升的兩個計算過程時，理想的迭代終點在殘差的拐點處，即當(dāng)errk/errk-1≥1時迭代終止；(4)在第k次計算過程中，往后逐次取n個迭代步，如果這n個迭代步的相對殘差都在增長，默認第k次的解為最優(yōu)解.

對于第一種終止條件，更多的是以經(jīng)驗為主，要么沒有達到精度要求，要么增加不必要的計算量；第二種終止條件在應(yīng)用于病態(tài)方程組時經(jīng)常失效(富明慧和李勇息，2018)；第三種終止條件太過理想，在實際計算中的殘差曲線往往并非單調(diào)，而且經(jīng)常有波動.綜上，本文選取第四種終止條件.

設(shè)

(12)

式中n為選擇的檢測窗口，為整數(shù)，一般情況下取2，最大不超過10(富明慧和李勇息，2018).

2 高斯脈沖的波場反變換

2.1 算法精度

首先研究一個簡單脈沖，令u(τ)=1

從圖2的結(jié)果可以看出，對于簡單脈沖u(τ)=1，本文的方法與正則化共軛梯度法都能得到相對較好的結(jié)果，但是相較而言，本文結(jié)果的精度比較高.

圖2 波場反變換結(jié)果和相對誤差 (a) 本文方法結(jié)果； (b) 本文方法相對誤差； (c) PRCG結(jié)果； (d) PRCG相對誤差.Fig.2 Inverse wave field transformation results and relative error (a) The present method; (b) The relative error in paper； (c) The PRCG method; (d) The relative error by PRCG.

接下來本文使用單個高斯脈沖和組合高斯脈沖的方式來得到更加復(fù)雜的虛擬波場.令u(τ)為高斯脈沖，它的表達式如下：

(13)

式中，a表示高斯脈沖的峰值大小，c表示高斯脈沖的寬度，τ0表示高斯脈沖的峰值位置，h0表示高斯脈沖在垂直方向的偏移量.通過給定不同的a、c、τ0和h0的值就可以得到不同樣式的高斯脈沖，而且還可以對上式進行組合得到更加復(fù)雜的高斯脈沖序列.

使用精細積分法分別計算了D型、G型、Q型、A型、H型和K型模型的虛擬波場(圖3)，并與前人研究結(jié)果進行了精度對比.從相對誤差圖上可以看出，本文方法的相對誤差較小，均在4%以下，而正則化共軛梯度法的相對誤差隨著模型的復(fù)雜而逐漸增大，最大相對誤差達到了將近50%，從這里可以看出本文方法具有較高的精度.

圖3 不同型模型波場反變換結(jié)果和相對誤差 (A) D型模型； (B) G型模型； (C) Q型模型; (D) A型模型; (E) H型模型; (F) K型模型.(a1—f1) 本文方法結(jié)果； (a2—f2) 本文方法相對誤差； (a3—f3) PRCG結(jié)果；(a4—f4) PRCG相對誤差.Fig.3 The inverse transformation results and relative errors of thedifferent model (A) D model; (B) Gmodel; (C) Q model; (D) A model; (E) H model; (F) K model; (a1—f1) The present method; (a2—f2) The relative error in paper; (a3—f3) The PRCG method; (a4—f4) The relative error by PRCG.

2.2 模型分辨率分析

令第一個脈沖的脈寬為0.02，振幅為2；第二個脈沖的脈寬為0.04，振幅為1.6；使第一個脈沖的波峰位置為τ01=0.03，改變第二個波峰τ02的位置，使其分別為0.04和0.08.結(jié)果如圖4所示：

本節(jié)研究了本文方法對Q型、A型、H型和K型模型的分辨能力(圖4).固定第一個脈沖的波峰位置為τ01=0.03，改變第二個波峰τ02的位置，使其分別為0.08和0.04，使得兩個波峰的位置越來越近.從上述結(jié)果可以看出，本文方法具有較好的分辨能力，而隨著兩個波峰的位置越來越近，正則化共軛梯度法的分辨能力越來越低.

圖4 波場反變換結(jié)果 (A) Q型模型； (B) A型模型; (C) H型模型; (D) K型模型; (a1—d1) 本文方法結(jié)果τ02=0.08；(a2—d2)PRCG結(jié)果τ02=0.08； (a3—d3) 本文方法結(jié)果τ02=0.04； (a4—d4) PRCG結(jié)果τ02=0.04.Fig.9 The inverse transformation results and relative errors of the different model (A) Q model； (B) A moel; (C) H model; (D) K model; (a1—d1) The present method τ02=0.08; (a2—d2) The PRCG method τ02=0.08; (a3—d3) The present method τ02=0.04; (a4—d4) The PRCG method τ02=0.04.

2.3 加噪試驗

分別對D型、G型、Q型、A型、H型和K型模型的擴散場中分別加入2%和5%的噪聲，并計算對應(yīng)的虛擬波場.計算結(jié)果如圖5.

圖5 不同模型加噪波場反變換結(jié)果 (A) D型模型; (B) G型模型; (C) Q型模型; (D) A型模型; (E) H型模型; (F) K型模型; (a1—f1) 本文方法結(jié)果SNR=50∶1； (a2—f2) PRCG結(jié)果SNR=50∶1； (a3—f3) 本文方法結(jié)果SNR=20∶1； (a4—f4) PRCG結(jié)果SNR=20∶1.Fig.5 The inverse transformation results and relative errors of the different model (A) D model; (B) G model; (C) Q model; (D) A model; (E) H model; (F) K model; (a1—f1) The present method SNR=50∶1; (a2—f2) The PRCG method SNR=50∶1; (a3—f3) The present method SNR=20∶1; (a4—f4) The PRCG method SNR=20∶1.

圖5各模型的擴散場信號中加入2%和5%的噪聲，測試了本文方法和正則化共軛梯度法的抗噪性.從圖中可以看出，本文的方法抗噪能力比較強，對于含有2%和5%噪聲的擴散場，依舊能求取較為準(zhǔn)確的虛擬波場；對于含有2%噪聲的擴散場，PRCG可以較為準(zhǔn)確地求取對應(yīng)的虛擬波場，對于含有5%噪聲的擴散場，求得的虛擬波場結(jié)果較差.當(dāng)信噪比超過5%時，兩種方法的虛擬波場結(jié)果波動劇烈，與真值差距較大.因此對于野外信號的采集必須控制噪聲在5%以下，并且對于局部不光滑數(shù)據(jù)不能直接求解，需要對數(shù)據(jù)進行平滑降噪后才能計算虛擬波場.

3 三維模型計算

為了進一步驗證本文使用的波場反變換算法的有效性，本節(jié)設(shè)計如下三維模型計算虛擬波場，并與正則化共軛梯度法進行比較.圖6為均勻半空間中含有兩個低阻薄板狀異常體的模型示意圖.均勻半空間的電阻率為200 Ωm，上層異常體上表面距離地表150 m，幾何尺寸為100 m×100 m×20 m，電阻率為5 Ωm；下層異常體上表面距地表250 m，幾何尺寸為100 m×100 m×30 m，電阻率為1 Ωm.接地導(dǎo)線源的長度為200 m，最短偏移距為100 m，發(fā)射電流為10A；布設(shè)13條測線，測線的線距為10 m，起始測線位置為y=-80 m，終止測線位置為y=50 m，每條測線長度為500 m，點距為25 m，每條測線21個測點，接收的電磁場信號為Ex.本文三維正演采用擬態(tài)有限體積法(Zhou et al.，2018).

圖6 三維模型示意圖 (a) XOZ平面圖； (b) XOY平面圖.Fig.6 Schematic diagram of three-dimensional model (a) XOZ plan; (b) XOY plan.

圖7為波場反變換結(jié)果，本文在這里只展現(xiàn)三條測線的波場反變換結(jié)果，分別是測線y=-70 m，測線y=0 m，測線y=30 m.選取這三條測線的原因是測線y=-70 m能反映出深部異常信息；測線y=0 m可以反映淺部和深部的異常信息；而測線y=50 m可以反映出淺部異常信息，其他測線的結(jié)果與這三條測線的結(jié)果類似，所以在此不進行展示.

圖7 本文方法波場反變換結(jié)果 (a) 測線y=-70 m； (b) 測線y=0 m； (c) 測線y=30 m.Fig.7 Results of wave field transformation in this paper (a) Line y=-70 m; (b) Line y=0 m; (c) Line y=30 m.

圖7a是測線y=-70 m的波場反變換結(jié)果，從圖中可以看出，除了直達波以外，還有深部異常體的上界面信息，由于瞬變電磁的高頻成分被地層迅速吸收，所以導(dǎo)致對深部異常體的分辨率較差；圖7b是測線y=0 m的波場反變換結(jié)果，從圖中可以看出，虛擬波場可以很好地反映出淺部異常體的上下界面，波形的同相軸比較明顯；對于深部異常體，圖7b也能反映異常體的上界面信息；圖7c是測線y=50 m的波場反變換結(jié)果，能夠反映出淺部異常體的界面信息.從三幅圖中可以看出，使用本文方法計算的結(jié)果，虛擬波場比較平滑，在整個時間段，無意義的波動起伏也比較少.

圖8為使用PRCG方法得到的波場反變換結(jié)果.可以看出，只有測線y=-70 m的結(jié)果比較好；測線y=0 m的結(jié)果主要反映的是深部異常體上界面的部分信息；測線y=50 m對淺層異常體界面反映較好，但是由于PRCG方法的精度較低，出現(xiàn)了多余的虛擬波場.

圖8 PRCG波場反變換結(jié)果 (a) 測線y=-70 m； (b) 測線y=0 m； (c) 測線y=30 m.Fig.8 Results of wave field transformation by PRCG (a) Line y=-70 m; (b) Line y=0 m; (c) Line y=30 m.

圖9是本文方法的虛擬波場三維切片圖，共有四張切片，對應(yīng)的位置分別為y=-50 m，y=40 m，x=120 m和x=240 m.從圖中可以看出，位置(1)和位置(2)處可以很清楚地反映異常體在xoy平面的位置，位置(3)和位置(4)處也能較好地反映出深部異常體在xoy平面的位置.通過計算，三維模型的虛擬波場與設(shè)計模型相符，可以進一步說明本文方法的有效性和準(zhǔn)確性.

圖9 波場反變換結(jié)果三維切片圖Fig.9 Three-dimensional slice map of wavefield transformation results

4 實測數(shù)據(jù)分析

上文利用各種算例與模型對本文方法進行了驗證，本節(jié)對甘肅省某煤田采空區(qū)實測數(shù)據(jù)進行處理.采區(qū)地層大部分被第四系黃土覆蓋，基本隱蔽,部分基巖裸露.地層從上到下由上白堊統(tǒng)、下白堊統(tǒng)、上侏羅統(tǒng)、中下侏羅統(tǒng)及上三疊統(tǒng)組成，中下侏羅統(tǒng)為主要煤系地層，上三疊統(tǒng)為煤系基底.采區(qū)主要含煤地層由石英細粗砂巖、砂礫巖、炭質(zhì)泥巖、炭質(zhì)粉砂巖組成.采區(qū)地表南部為丘陵，部分丘陵的白堊紀巖層出露.地表北部開闊平緩，為農(nóng)田及低矮丘陵，小沖溝發(fā)育.地表中部、北部有兩條季節(jié)性沙河.

野外工作使用V8多功能電法儀器，發(fā)射系統(tǒng)采用地面長導(dǎo)線源，發(fā)射電流15 A，基頻8.3 Hz，共15條測線；每條測線400 m，點距10 m，共41個測點，沿南北向布設(shè)，接收Ex電磁信號.本文選取380號測線數(shù)據(jù)進行視電阻率成像(張瑩瑩等，2015)與波場反變換，其余測線不在此文贅述.

由圖10a視電阻率剖面圖可以看出，隨著深度的增加，視電阻率值總體呈現(xiàn)先減小后增大的趨勢.在測線的280～320 m、340～400 m、440～500 m處，高程在800～1000 m之間，有三個低阻異常區(qū)，推斷可能為煤礦采空區(qū)塌陷，導(dǎo)致裂隙發(fā)育形成富水區(qū)域，從而形成低阻異常區(qū)，與已知的鉆孔資料較為吻合(已知鉆孔的煤線高程為1042.93 m).由圖10b可以看出，虛擬波場的同相軸與視電阻率的趨勢基本一致，整體呈南低北高，與煤層底板等值線圖較吻合(圖11).在測線的280～320 m、340～380 m、400～430 m，430～500 m處，由于煤層被開采，地層錯斷較為嚴重，在虛擬波場圖上表現(xiàn)出橫向不連續(xù)，間斷跳躍的特征.本文通過與已知資料對比，本文的計算結(jié)果與已知資料相吻合，說明本文方法在識別電性界面方面有較好的效果.

圖10 380號測線視電阻率圖(a) 與波場反變換結(jié)果(b)Fig.10 Apparent resistivity profile and pseudo wave field of line 380

圖11 工區(qū)煤層底板等值線圖Fig.11 Contour map of coal seam floor in working area

5 結(jié)論

本文通過使用精細積分法，成功實現(xiàn)了從瞬變電磁擴散場到虛擬波場全時域穩(wěn)定的反變換算法.在保證計算精度的前提下，使用指數(shù)級的積分步長可以顯著提高計算效率.通過設(shè)置合理的終止迭代條件，使得本文的方法適用性更廣.

實現(xiàn)波場反變換后，本文對算法進行了驗證.通過對典型地電模型的計算，證明了本文的方法具有較好的精度(得到虛擬波場值與理論值的誤差小于4%)、分辨能力以及抗噪性，可以說明本文方法適用于復(fù)雜介質(zhì)模型.通過計算三維模型的虛擬波場值，可以看出虛擬波場對淺層異常體的分辨率較好，隨著深度增加，高頻成分很快被吸收，低頻成分占主導(dǎo)，波場幅值減小，子波寬度增加，分辨率逐漸降低.最后通過對實測數(shù)據(jù)的處理，進一步說明了本文方法的可行性，為后續(xù)的瞬變電磁構(gòu)造成像方法打下了良好的基礎(chǔ).

從文中可以看出，當(dāng)擴散場包含的噪聲越來越大時，波場反變換的結(jié)果與實際結(jié)果相差也就越大.因此，在后續(xù)工作中可以引入窗口掃描、相關(guān)疊加的方法進一步提高該方法的抗噪性.

致謝感謝中山大學(xué)富明慧教授在作者研究過程中給予指導(dǎo)，同時感謝論文評審專家和編輯提出的寶貴修改意見

附錄A

(A1)

(A2)

類似地，我們可以得到：

(A3)

此時對兩邊同時乘以b，得：

(A4)

當(dāng)k→∞時，我們將得到精確解xk.但是精細積分的步長ζ一般都比較小(富明慧和李勇息，2018)，如果對積分區(qū)間[0,∞)采用線性等間隔離散，計算量將大到無法接受.所以本文參考前人工作(鐘萬勰，1994)，采用指數(shù)增長的積分步長.

對于矩陣指數(shù)函數(shù)，鐘萬勰提出了一種精確且高效的計算方法，在各類偏微分方程中獲得了廣泛應(yīng)用.與上文類似，設(shè)

(A5)

對上式兩邊同時乘以b，得：

x1=(I+exp(-Aζ))x0.

(A6)

由于ζ是極小量，在計算矩陣指數(shù)函數(shù)時，可使用Taylor展開的前幾項：

(A7)

對式(A7)在區(qū)間[0,ζ]積分，得

(A8)

綜上，x的初始值x0=F(ζ)b.

令積分步長以2k-1ζ變化，則有

k=1,2,…

(A9)

且當(dāng)2kζ→Fk時，有

(A10)

(A11)

上式中Tk即為矩陣指數(shù)函數(shù)的值，這樣計算矩陣指數(shù)函數(shù)的好處在于：在計算過程中當(dāng)k很小時，能將大量和小量分離，避免小量與大量相加時出現(xiàn)數(shù)據(jù)“湮沒”，保證計算精度(富明慧和李勇息，2018).