周武，羅威，藍(lán)星，簡(jiǎn)興祥

(1.甘肅省交通規(guī)劃勘察設(shè)計(jì)院股份有限公司，甘肅 蘭州 730030； 2.成都理工大學(xué) 地球物理學(xué)院，四川 成都 610000； 3.四川省冶勘設(shè)計(jì)集團(tuán)有限公司，四川 成都 610000)

0 引言

大地電磁測(cè)深(MT)是探測(cè)地球深部電性結(jié)構(gòu)的主要方法，通過在地表同步觀測(cè)電場(chǎng)、磁場(chǎng)分量，定性或定量地對(duì)深部地質(zhì)構(gòu)造進(jìn)行研究，在油氣勘探、固體礦產(chǎn)資源勘察、深部地質(zhì)結(jié)構(gòu)探測(cè)、地?zé)岷偷叵滤Y源調(diào)查、地震預(yù)測(cè)和地質(zhì)災(zāi)害防治等領(lǐng)域應(yīng)用廣泛。在MT二維正演模擬方面，主要計(jì)算方法有邊界元法、有限單元法和有限差分法等。其中，有限差分法以實(shí)現(xiàn)簡(jiǎn)單、計(jì)算量較小和實(shí)用性強(qiáng)的優(yōu)點(diǎn)使用最為廣泛[1-3]，但過去的二維有限差分模擬大都將電場(chǎng)和磁場(chǎng)采樣點(diǎn)放置于相同位置，而忽略了可能存在的散度問題。交錯(cuò)采樣網(wǎng)格的最大特點(diǎn)是能自動(dòng)保證電磁場(chǎng)分布遵守能量守恒定律[4]，由于三維電磁場(chǎng)散度一般不為零，交錯(cuò)采樣網(wǎng)格主要在大地電磁三維正演中被采用[5-6]，而二維情況下的大地電磁交錯(cuò)網(wǎng)格有限差分正演問題研究目前相對(duì)較少[7]。

在大地電磁反演問題上，目前基于一階梯度方向或二階牛頓方向收斂的各種最優(yōu)化反演理論已經(jīng)較為成熟，因此只要擁有良好的正演方法，反演需要考慮的關(guān)鍵問題就是如何計(jì)算或處理靈敏度矩陣或海賽矩陣，然后再套用相對(duì)成熟的反演方法。目前，大地電磁中主要存在3種主流線性反演方法：OCCAM及其改進(jìn)型、非線性共軛梯度法(NLCG)和有限內(nèi)存擬牛頓法(LBFGS)。OCCAM反演[8]及其改進(jìn)型是基于梯度收斂的最優(yōu)化法，與最小二乘不同，OCCAM在每次迭代中都要求取最佳的目標(biāo)函數(shù)正則化因子，因此其反演比較穩(wěn)定，但計(jì)算速度很慢，同時(shí)由于在反演過程中必須存儲(chǔ)全部靈敏度矩陣，內(nèi)存需求較大。NLCG[9]是介于最速下降法和牛頓法之間的方法，利用梯度的共軛方向提升收斂效率，避免求取二階海賽矩陣；為避免存儲(chǔ)靈敏度矩陣，NLCG利用正演的稀疏矩陣求取梯度向量，因此對(duì)內(nèi)存需求很小，但NLCG在每次迭代中都需進(jìn)行線搜索求取最佳步長(zhǎng)，因此一次反演迭代中有可能進(jìn)行多次正演。LBFGS[10]是由擬牛頓法發(fā)展起來的一種反演方法，擬牛頓法根據(jù)更新海賽矩陣的秩不同有非常多的變種，主要有DFP法和BFGS法。LBFGS又是BFGS法的變種，其特點(diǎn)是使用有限次迭代的梯度和修正量信息來更新海賽矩陣，對(duì)內(nèi)存需求小，同時(shí)基本和BFGS迭代一樣穩(wěn)定。LBFGS和NLCG相比，二者內(nèi)存需求差不多，求取梯度的過程一樣，唯一的不同是LBFGS在迭代中后期基本上每次迭代的修正方向都為牛頓方向，其修正步長(zhǎng)基本為1，這樣就可一定程度上節(jié)省線搜索時(shí)間，其用于線搜索的正演次數(shù)約為NLCG的一半。近年來，LBFGS法在大地電磁二三維反演中開展了諸多研究[7,11-13]，已成為目前大地電磁二三維反演的主要方法之一。

本文首先利用交錯(cuò)采樣網(wǎng)格實(shí)現(xiàn)大地電磁二維正演，然后基于LBFGS法，實(shí)現(xiàn)大地電磁二維交錯(cuò)采樣網(wǎng)格有限差分LBFGS反演，目標(biāo)是實(shí)現(xiàn)一種精度高、效率高的大地電磁二維正反演算法，并能有效應(yīng)用于實(shí)測(cè)大地電磁資料的反演。

1 大地電磁交錯(cuò)網(wǎng)格有限差分二維正演

在大地電磁研究的頻率范圍內(nèi)(約105～10-4Hz)，由于頻率相對(duì)較低，滿足似穩(wěn)電磁場(chǎng)，忽略位移電流后的Maxwell積分形式方程為：

∮H·dl=?σE·dS，

(1)

∮E·dl=?iωμ0H·dS。

(2)

式中：E和H分別為電場(chǎng)強(qiáng)度和磁場(chǎng)強(qiáng)度矢量，真空磁導(dǎo)率μ0=4π×10-7H/m，σ為介質(zhì)電導(dǎo)率，ω為圓頻率。

對(duì)于二維大地大地電磁，需要垂直測(cè)線方向無限延伸。對(duì)于第(i,k)個(gè)單元的電場(chǎng)分量，根據(jù)式(1)可由周邊的4個(gè)磁場(chǎng)求得，這4個(gè)磁場(chǎng)分量，根據(jù)式(2)又可由各自周邊的2個(gè)電場(chǎng)求得。綜合起來，即這單元的電場(chǎng)可由其周邊的4個(gè)電場(chǎng)分量求得，然后按這種關(guān)系以此類推，將所有單元的電場(chǎng)與各自周邊4個(gè)電場(chǎng)的線性方程組合起來，就構(gòu)建了剛度矩陣。圖1為直角坐標(biāo)系中的第(i,k)個(gè)單元差分示意。

1.1 TE模式正演

對(duì)于TE模式，以Ey(i,k)為中心，式(1)Maxwell第1個(gè)方程展開為：

Hz(i,k)Δz(k)+Hx(i,k+1)Δx(i)-

Hz(i+1,k)Δz(k)-Hx(i,k)Δx(i)=

σ(i,k)Ey(i,k)Δz(k)Δx(i) 。

(3)

對(duì)式(3)中的4個(gè)磁場(chǎng)，由式(2)又可分別用周圍的電場(chǎng)求得：

Ey(i,k)-Ey(i-1,k)=

(4)

Ey(i,k-1)-Ey(i,k)=

(5)

Ey(i+1,k)-Ey(i,k)=

(6)

圖1 MT二維交錯(cuò)網(wǎng)格有限差分離散示意Fig.1 MT 2-D staggered grid finite difference separation diagram

Ey(i,k)-Ey(i,k+1)=

。(7)

將這4個(gè)關(guān)系式代入式(3)中即可消去磁場(chǎng)分量，最終得到Ey(i,k)與其周邊4個(gè)電場(chǎng)的關(guān)系式。用C1～C5的系數(shù)表示為

C1Ey(i-1,k)+C2Ey(i,k)+C3Ey(i+1,k)+

C4Ey(i,k-1)+C5Ey(i,k+1)=0 ，

(8)

其中：

C2=-C1-C5-C4-C3-σ(i,k)iωμ0Δz(k)Δx(i) ，

1.2 TM模式正演

對(duì)于TM模式，以Hy(i,k)為中心，式(2)Maxwell第2個(gè)方程展開為

Ez(i,k)Δz(k)+Ex(i,k+1)Δx(i)-

Ez(i+1,k)Δz(k)-Ex(i,k)Δx(i)=

iωμ0Hy(i,k)Δz(k)Δx(i) 。

(9)

對(duì)式(9)中的4個(gè)電場(chǎng)，由式(1)又可分別用周圍的磁場(chǎng)求得：

Hy(i,k)-Hy(i-1,k)=

(10)

Hy(i,k-1)-Hy(i,k)=

(11)

Hy(i+1,k)-Hy(i,k)=

(12)

Hy(i,k)-Hy(i,k+1)=

(13)

將這4個(gè)關(guān)系式代入式(9)中即可消去電場(chǎng)分量，最終得到Hy(i,k)與其周邊4個(gè)磁場(chǎng)的關(guān)系式。用C1～C5的系數(shù)表示為：

C1Hy(i-1,k)+C2Hy(i,k)+C3Hy(i+1,k)+

C4Hy(i,k-1)+C5Hy(i,k+1)=0 ，

(14)

其中：

C2=-C1-C5-C3-C4-iωμ0Δz(k)Δx(i) ，

1.3 方程組

根據(jù)前述TE和TM離散方案，最終可以獲得關(guān)于Ey和Hy的離散方程：

(15)

式中：Ey和Hy為求解向量，be、bh對(duì)應(yīng)的右端項(xiàng)向量，Ae、Ah分別為式(8)、式(14)中不同單元相關(guān)的C系數(shù)組成。

要想求解研究區(qū)域內(nèi)所有采樣點(diǎn)上的場(chǎng)值，還需給出研究區(qū)邊界上場(chǎng)值，交錯(cuò)采樣網(wǎng)格二維邊界條件可簡(jiǎn)要概括為表1所示形式[2]，其中頂邊界為第一類邊界條件(其中TE模式頂面位于空氣頂層，TM模式頂面位于地表)，側(cè)邊界為第二類邊界條件，底邊界對(duì)應(yīng)第三類邊界條件(其中σ為底層介質(zhì)電導(dǎo)率，Zbottom為底層頂面阻抗)。對(duì)構(gòu)建的方程組，由于直接解法可獲得精確解，因此這里采用基于LU分解的直接法求解。

表1 邊界條件形式

2 LBFGS反演

根據(jù)Tikhonov的正則化反演理論[14]，反演問題可表述為

ψ(m)=ψd(m)+λψm(m)=

(WdΔd)T(WdΔd)+λ(Wmm)T(Wmm) 。

(16)

式中：m為模型;ψd(m)為數(shù)據(jù)目標(biāo)函數(shù);ψm(m)為模型約束目標(biāo)函數(shù);λ為正則化因子;Δd為觀測(cè)數(shù)據(jù)與模型響應(yīng)的差向量;Wd為數(shù)據(jù)加權(quán)矩陣;Wm為模型約束矩陣，目的是使反演得到的模型具備所設(shè)定的相應(yīng)結(jié)構(gòu)特征，一般采用二階拉普拉斯算子。

在確立反演目標(biāo)函數(shù)后，接下來的工作就是選擇合適的方法進(jìn)行最優(yōu)化反演。大地電磁反演方法較多，維性、反演效率以及對(duì)計(jì)算資源需求等因素影響著反演方法的選擇。經(jīng)典的最優(yōu)化方法基本都是基于梯度下降法和牛頓法的改進(jìn)或重構(gòu)，在大地電磁中幾乎所有的最優(yōu)化反演方法都是這兩種方法的實(shí)用化或變種，其中有限內(nèi)存擬牛頓法(LBFGS)因其特殊的優(yōu)勢(shì)，近年來逐漸應(yīng)用廣泛并受到認(rèn)可[12-15]。LBFGS算法流程如下:

1)設(shè)定k=0，初始模型m0和精度閾值ε；

2)設(shè)置初始的正定矩陣H0=I；

5)若‖gk+1‖<ε，則算法結(jié)束；

6)利用m個(gè)已保存的s、y向量對(duì)計(jì)算Hk+1：

7)迭代次數(shù)k=k+1，跳轉(zhuǎn)到步驟3繼續(xù)。

3 算例分析

3.1 層狀模型正演

由于層狀介質(zhì)模型可通過一維正演求得其解析解，因此二維正演算法的數(shù)值解可采用一維解析解進(jìn)行對(duì)比驗(yàn)證。層狀模型第一層電阻率100 Ω·m，厚1 km，第二層電阻率1 000 Ω·m，厚9 km，第三層電阻率10 Ω·m。

計(jì)算過程中非擴(kuò)展區(qū)域網(wǎng)格數(shù)為橫向50列，每列寬度10 km，縱向80層，層深度對(duì)數(shù)遞增，網(wǎng)格剖分單元數(shù)4 000個(gè)，測(cè)點(diǎn)均勻分布在地表每列單元中心共100個(gè)，計(jì)算了104～10-4Hz共61個(gè)頻率的視電阻率和阻抗相位響應(yīng)。圖2給出了解析解與剖面中心位置測(cè)點(diǎn)的視電阻率和阻抗相位對(duì)比曲線，可以看出數(shù)值結(jié)果和解析解曲線基本重合，僅在約1 000 Hz以上高頻部分存在細(xì)小差異，但這是有限差分?jǐn)?shù)值計(jì)算普遍均存在的現(xiàn)象。二維數(shù)值解與解析解視電阻率平均相對(duì)誤差為1.19%，相位平均相對(duì)誤差為0.72%，表明相對(duì)誤差較小，二維正演結(jié)果精度較高。

圖2 層狀模型二維正演中心測(cè)點(diǎn)視電阻率和阻抗相位Fig.2 Apparent resistivity and impedance phase of 2-D forward modeling center of layered model

3.2 異常體模型正反演

為進(jìn)一步測(cè)試正演算法的正確性以及反演算法的穩(wěn)定性和效率，設(shè)計(jì)了較為復(fù)雜的異常體模型。模型的設(shè)計(jì)參考借鑒Egbert對(duì)Modem算法的測(cè)試模型[12]，設(shè)計(jì)由淺至深為2層，橫向上3列，共6個(gè)交錯(cuò)分布的高低阻異常體，模型共80層×50列=4 000個(gè)單元，50個(gè)地表測(cè)點(diǎn)按100 等距分布，圖3為模型的具體形態(tài)。模型背景電阻率100 Ω·m，縱向上在150～300 m、800～1 500 m深度存在異常體，位置交錯(cuò)分布，電阻率值分別為1 000 Ω·m、10 Ω·m。計(jì)算了104～1 Hz共30個(gè)頻點(diǎn)的正演響應(yīng)，圖4為正演響應(yīng)擬斷面，可以看到，無論視電阻率還是阻抗相位，都對(duì)高阻和低阻異常有較明顯的反映。

圖3 異常體模型Fig.3 Anomalous body model

圖4 異常體模型正演擬斷面Fig.4 Forward section of abnormal body model

對(duì)異常體模型進(jìn)行反演計(jì)算。反演前在正演數(shù)據(jù)中加入3%的高斯隨機(jī)噪音，并將其作為數(shù)據(jù)協(xié)方差，反演正則化因子設(shè)為 0.01，初始模型為100 Ω·m的均勻半空間。反演終止條件為擬合差RMS不大于1，這里的擬合差計(jì)算公式為：

(17)

分別對(duì)合成數(shù)據(jù)開展TE、TM和TE+TM反演，最終反演的步長(zhǎng)曲線和擬合差曲線如圖5所示。其中，TE模式反演步長(zhǎng)經(jīng)過8次后穩(wěn)定為1，TM模式9次后穩(wěn)定，TE+TM模式29次迭代后穩(wěn)定為1。3種模式反演的擬合差RMS從初始的約11開始，經(jīng)過幾次迭代后快速衰減到3左右，隨后緩慢降低，其中TE、TM和TE+TM分別在第17、23、40次迭代后滿足收斂條件。從反演步長(zhǎng)曲線和擬合差曲線看出反演較為穩(wěn)定，符合LBFGS收斂特征。

TE、TM及TE+TM3種模式LBFGS最終反演結(jié)果如圖6所示。與真實(shí)模型對(duì)比看出，3種模式的反演結(jié)果總體上與真實(shí)模型都較為吻合：TE模式對(duì)異常體厚度反映較好，但橫向邊界相對(duì)一般；TM模式對(duì)異常橫向范圍反映較好，但縱向分辨相對(duì)一般；TE+TM模式綜合了TE和TM的優(yōu)點(diǎn)，對(duì)異常體反映較好。上述特征符合MT二維反演規(guī)律，表明算法可靠。

另外，為分析LBFGS反演的效率，將該模型同時(shí)開展了NLCG二維反演，表2給出了2種算法的時(shí)間對(duì)比。對(duì)于NLCG反演，TE、TM、TE+TM的10次迭代耗時(shí)分別為75、57、133 s，而LBFGS的10次迭代耗時(shí)分別為59、45、88 s，相對(duì)有約27%的效率提升。NLCG擬合差達(dá)到1的計(jì)算時(shí)間分別為175、92、506 s，而LBFGS只用了125、73、371 s，相對(duì)有約25%的效率提升。這表明LBFGS反演效率總體上比NLCG高。

圖5 異常體模型二維反演步長(zhǎng)和擬合差曲線Fig.5 2-D inversion step length and fitting difference curve of abnormal body model

圖6 異常體模型二維反演Fig.6 2-D inversion of anomalous body model

表2 LBFGS和NLCG反演效率對(duì)比

3.3 實(shí)測(cè)資料反演

圖8為采用本文算法對(duì)實(shí)測(cè)MT資料的TM模式反演結(jié)果，分析看出，剖面范圍內(nèi)發(fā)育多條斷層，斷層位置和形態(tài)與地表地質(zhì)調(diào)查以及氡氣測(cè)量結(jié)果較為一致。此外，地表的諸多斷裂在約1 000 m深度匯集，形成F1和F2兩個(gè)主要斷層，發(fā)育深度均大于2 000 m，其中F1傾向南西、 F2傾向北東，2個(gè)斷裂在深部交匯，斷層深部影響帶范圍較大。F1和F2斷層帶中間低阻現(xiàn)象較為明顯，表明斷層帶上盤局部巖體破碎帶較發(fā)育，為地?zé)崃黧w的賦存提供了一定的儲(chǔ)集空間，具有一定的導(dǎo)水、導(dǎo)熱作用。地?zé)崃黧w沿?cái)鄬訋仙髢?chǔ)存在熱儲(chǔ)含水層內(nèi)，在上覆隔水保溫蓋層的作用下，形成地?zé)豳Y源。通過實(shí)測(cè)MT資料的反演和解釋，可以看出本文算法的反演結(jié)果揭露的特征，與地質(zhì)和相關(guān)資料較為一致。

圖7 MT測(cè)線分布Fig.7 MT line distribution

圖8 MT反演結(jié)果Fig.8 MT inversion results

4 結(jié)論及討論

本文借鑒三維大地電磁正演模擬中使用的交錯(cuò)采樣網(wǎng)格方法，采用交錯(cuò)采樣網(wǎng)格開展了大地電磁二維有限差分正演研究，模型計(jì)算表明該正演算法具有較高的計(jì)算精度?；诮诲e(cuò)采樣網(wǎng)格大地電磁二維有限差分正演，利用有限內(nèi)存擬牛頓算法，實(shí)現(xiàn)了大地電磁二維快速反演。理論模型計(jì)算表明有限內(nèi)存擬牛頓反演算法具有較好的穩(wěn)定性，與非線性共軛梯度算法對(duì)比發(fā)現(xiàn)，LBFGS計(jì)算效率至少大于NLCG效率的1.2倍。將反演算法在宕昌縣官鵝溝地?zé)峥辈轫?xiàng)目中開展了大地電磁資料二維反演解釋，結(jié)果表明本文算法與現(xiàn)有資料和認(rèn)識(shí)較為一致，具有一定的實(shí)用價(jià)值。

由于算法暫未考慮地形以及各向異性介質(zhì)的影響，建議下一步開展基于交錯(cuò)采樣網(wǎng)格的起伏地形、各向異性正反演研究。此外，盡管LBFGS具有較高的計(jì)算效率，但隨著數(shù)據(jù)量和模型單元的增大，并行計(jì)算研究仍然值得開展。