大地電磁交錯采樣有限差分二維正反演研究

周武，羅威，藍星，簡興祥

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

0 引言

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

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

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

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

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

∮H·dl=?σE·dS，

(1)

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

(2)

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

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

1.1 TE模式正演

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

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)

對式(3)中的4個磁場，由式(2)又可分別用周圍的電場求得：

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二維交錯網(wǎng)格有限差分離散示意Fig.1 MT 2-D staggered grid finite difference separation diagram

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

。(7)

將這4個關系式代入式(3)中即可消去磁場分量，最終得到Ey(i,k)與其周邊4個電場的關系式。用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模式正演

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

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)

對式(9)中的4個電場，由式(1)又可分別用周圍的磁場求得：

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個關系式代入式(9)中即可消去電場分量，最終得到Hy(i,k)與其周邊4個磁場的關系式。用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離散方案，最終可以獲得關于Ey和Hy的離散方程：

(15)

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

要想求解研究區(qū)域內(nèi)所有采樣點上的場值，還需給出研究區(qū)邊界上場值，交錯采樣網(wǎng)格二維邊界條件可簡要概括為表1所示形式[2]，其中頂邊界為第一類邊界條件(其中TE模式頂面位于空氣頂層，TM模式頂面位于地表)，側(cè)邊界為第二類邊界條件，底邊界對應第三類邊界條件(其中σ為底層介質(zhì)電導率，Zbottom為底層頂面阻抗)。對構建的方程組，由于直接解法可獲得精確解，因此這里采用基于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ù)目標函數(shù);ψm(m)為模型約束目標函數(shù);λ為正則化因子;Δd為觀測數(shù)據(jù)與模型響應的差向量;Wd為數(shù)據(jù)加權矩陣;Wm為模型約束矩陣，目的是使反演得到的模型具備所設定的相應結構特征，一般采用二階拉普拉斯算子。

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

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

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

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

6)利用m個已保存的s、y向量對計算Hk+1：

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

3 算例分析

3.1 層狀模型正演

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

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

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

3.2 異常體模型正反演

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

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

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

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

(17)

分別對合成數(shù)據(jù)開展TE、TM和TE+TM反演，最終反演的步長曲線和擬合差曲線如圖5所示。其中，TE模式反演步長經(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次迭代后滿足收斂條件。從反演步長曲線和擬合差曲線看出反演較為穩(wěn)定，符合LBFGS收斂特征。

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

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

圖5 異常體模型二維反演步長和擬合差曲線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反演效率對比

3.3 實測資料反演

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

圖7 MT測線分布Fig.7 MT line distribution

圖8 MT反演結果Fig.8 MT inversion results

4 結論及討論

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

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