姜奮勇，葉益信，陳海文，楊爍健

(東華理工大學(xué) 地球物理與測控技術(shù)學(xué)院，江西 南昌 330013)

0 引言

在過去的幾十年中，大地電磁法被廣泛應(yīng)用于深部構(gòu)造研究以及油氣和礦產(chǎn)勘查中[1-3]。常規(guī)的MT數(shù)據(jù)處理和解釋往往忽略地形因素的影響，而在實際MT數(shù)據(jù)采集中，難以避免的會遇到地形起伏的問題，而起伏地表會引起MT響應(yīng)的畸變[4-5]。因此，將地形因素引入到大地電磁數(shù)值模擬中是非常必要的，許多學(xué)者利用有限元法對起伏地形下二維大地電磁響應(yīng)進行了數(shù)值模擬研究，并對畸變后的大地電磁響應(yīng)做了進一步的分析[3,6]。常規(guī)基于規(guī)則網(wǎng)格的算法對復(fù)雜起伏地形的模擬精度不高，而且使用起來不方便，因此越來越多的學(xué)者采用非結(jié)構(gòu)三角網(wǎng)格有限元方法計算電磁場[7-9]。非結(jié)構(gòu)三角網(wǎng)格能夠準(zhǔn)確地模擬起伏地形和復(fù)雜地質(zhì)構(gòu)造，并且可以通過自適應(yīng)有限元法自動調(diào)節(jié)網(wǎng)格，以保證正演響應(yīng)的精度[10-15]。已有研究將非結(jié)構(gòu)三角網(wǎng)格自適應(yīng)有限元法用于求解大地電磁的響應(yīng)，取得了很好的效果，并結(jié)合Occam反演算法成功用于帶地形大地電磁數(shù)據(jù)二維反演中，證明了其算法具有較強的實用性和較高的精度[16-17]。在一些研究區(qū)內(nèi)地形比較復(fù)雜，水平地形的反演往往會產(chǎn)生一些不準(zhǔn)確的結(jié)果，要實現(xiàn)高質(zhì)量的探測，就要對研究區(qū)域采集的MT數(shù)據(jù)進行帶地形的反演，基于非結(jié)構(gòu)網(wǎng)格自適應(yīng)有限元的帶地形反演結(jié)果將會提供更加準(zhǔn)確的電阻率信息，這對以后的工程勘查和礦產(chǎn)勘探是很有幫助的。

本文對基于自適應(yīng)非結(jié)構(gòu)網(wǎng)格有限元的MT二維帶地形反演做了應(yīng)用研究，對起伏地形模型進行了MT二維反演試算，得到了比較好的結(jié)果。最后將本算法應(yīng)用于克拉瑪依后山區(qū)域?qū)崪yMT數(shù)據(jù)的反演，同樣取得了較好的結(jié)果。

1 反演理論

1.1 Occam反演

MT數(shù)據(jù)反演方法采用的是Occam算法[18-19]，該算法是基于正則化約束的最小二乘法，其泛函U表達式為：

U=μ‖Rm‖2+‖W(d-F(m))‖2，

(1)

式中：m為N維的模型參數(shù)矢量，一般為電阻率值；R為粗糙度算子矩陣；μ為拉格朗日乘子，用以平衡模型粗糙度和數(shù)據(jù)擬合誤差，當(dāng)μ取較大值時，反演以搜索光滑模型為主，反之則以搜索最小擬合誤差為主；W為關(guān)聯(lián)的對角加權(quán)矩陣；d為觀測數(shù)據(jù)矢量，F(xiàn)(m)為模型m對應(yīng)的正演響應(yīng)。

給定初始模型mk的函數(shù)，通過以下迭代方法實現(xiàn)目標(biāo)函數(shù)的最小化：

mk+1=[μRTR+(WJk)TWJk]-1×(WJk)TWd,

(2)

其中，修改數(shù)據(jù)向量

(3)

雅可比矩陣J為M×N大小的矩陣，其中每個分量的表達式為：

(4)

1.2 靈敏度計算

靈敏度矩陣的計算對于電磁反演是非常重要的步驟，能夠精確求解靈敏度矩陣，對于最終的電磁數(shù)據(jù)反演分辨率至關(guān)重要。求解靈敏度矩陣的方法有很多，最經(jīng)典的方法是將問題線性化，然后迭代求解[20]。而在我們的反演問題中，往往模型參數(shù)的數(shù)量要遠大于數(shù)據(jù)的數(shù)量，所以這種方法比較低效，本文采用的是伴隨方程法推導(dǎo)靈敏度矩陣[21-22]。

頻率域麥克斯韋方程組表示為：

(5)

式中：ε為介電常數(shù)；μ為磁化率；σ為可變電導(dǎo)率；E為電場強度；H為磁場強度；Je和Jm分別為電源和磁源。

設(shè)置合適的邊界條件，將φj(r)作為一個基函數(shù)，則電導(dǎo)率被表示為一組基函數(shù)的線性組合：

(6)

式中：σj為系數(shù)。將式(6)帶入方程組(5)，然后對σj求偏導(dǎo)，得到靈敏度方程：

(7)

(8)

考慮邊界問題并帶入其中，得到下式：

(9)

(10)

1.3 粗糙度計算

模型粗糙度算子R主要用于反演結(jié)果的穩(wěn)定，避免產(chǎn)生虛假的電性結(jié)構(gòu)，目前比較常用的是將模型梯度的L2范數(shù)作為模型的粗糙度算子[23-24]，其表達式為：

(11)

基于結(jié)構(gòu)網(wǎng)格的對積分的離散近似較為容易，而在非結(jié)構(gòu)網(wǎng)格中情況卻較為復(fù)雜，因此本文采用加權(quán)平方和的方法[17]處理梯度點積的積分：

(12)

其中：

Δmij=mi-mj,

(13)

(14)

(15)

式中：Vi表示第i個三角單元的面積；Ni表示與第i個非結(jié)構(gòu)三角單元共享一個頂點所有三角單元的集合；Δrij表示三角形之間質(zhì)心的距離。

2 模型實例

2.1 陸地起伏地形模型

為了驗證該方法反演的效果，本文首先設(shè)計了一個起伏地形的復(fù)雜模型。模型結(jié)構(gòu)如圖1所示，在不同的地層中分別設(shè)置了兩個低阻異常體，1號低阻異常體電阻率為3 Ω·m，2號低阻異常體電阻率為0.1 Ω·m。MT接收點從左往右依次設(shè)置了21個測點，如圖1中白色小三角形所示，每個測點間距800 m。MT的頻率范圍為0.001～1 000 Hz，共選取了48個頻點，取對數(shù)等間隔分布。

圖1 陸地起伏地形模型Fig.1 Land complex topography model

對上述起伏地形模型進行自適應(yīng)有限元正演計算，得到模型正演響應(yīng)與自適應(yīng)細化網(wǎng)格，如圖2所示，圖2a為第一次自適應(yīng)細化網(wǎng)格；圖2b為第五次自適應(yīng)細化網(wǎng)格，表1為細化次數(shù)及每次細化生成的節(jié)點和單元數(shù)量。可以看出，由于測點附近的網(wǎng)格剖分對測點接收到的電磁場分量影響較大，單元的誤差估計值較大，所以在自適應(yīng)網(wǎng)格優(yōu)化過程中，對淺地表和測點附近的網(wǎng)格進行了加密，而邊界和深部的網(wǎng)格單元誤差很小，所以并沒有進行加密。

表1 陸地起伏地形模型網(wǎng)格自適應(yīng)細化次數(shù)Table 1 Mesh adaptive refinement of land complex topography model

圖2 陸地起伏地形模型不通細化次數(shù)網(wǎng)格剖面Fig.2 Sections of different mesh refinement times of land complex topography model

模型進行正演計算后，得到了TE、TM模式的視電阻率和相位值，然后添加4%的隨機噪聲，則生成了反演計算所需的觀測數(shù)據(jù)。反演的初始模型包含空氣層和均勻地層，其中空氣不參與反演，給定均勻地層的初始電阻率為1.0 Ω·m。計算區(qū)域所劃分的單元越多，所需要的計算量越大，因此只對目標(biāo)反演區(qū)域采用較精細的網(wǎng)格剖分，其余部分則采用粗網(wǎng)格以減少不必要的計算量。如圖3所示，在y范圍為-10～10 km及地表到10 km深的目標(biāo)區(qū)域，調(diào)用Triangle程序?qū)υ搮^(qū)域進行精細網(wǎng)格剖分，共生成19 514個三角單元，而在其余區(qū)域共剖分產(chǎn)生了1 130個三角形單元。

圖3 反演初始模型及其網(wǎng)格剖分Fig.3 The initial model and its grid subdivision used in the inversion

利用并行機進行反演試算，如表2所示，反演經(jīng)過69次迭代，耗時7 h，占用內(nèi)存約為12 GB,最終模型均方根擬合差(RMS)值為4.91，反演結(jié)果如圖4所示。從圖中可以看出，該算法能夠準(zhǔn)確地反演出兩個低阻異常體在復(fù)雜模型中的埋深、大小和形態(tài)，同時帶地形的MT反演能夠較好地區(qū)分不同電性的地層，這也說明MT無論是對局部低阻異常體還是對較大尺度的電性結(jié)構(gòu)都能有效地探測。

表2 反演參數(shù)Table 2 Inversion parameters

圖4 陸地起伏地形模型反演結(jié)果Fig.4 Inversion result of land complex topography model

圖5為真實模型響應(yīng)數(shù)據(jù)與反演最終模型響應(yīng)數(shù)據(jù)的對比，由圖5可以看出，真實模型的響應(yīng)數(shù)據(jù)與反演模型得到的視電阻率和相位信息高度吻合，驗證了反演結(jié)果的可靠性和準(zhǔn)確度。

2.2 海底起伏地形模型

設(shè)計了一個起伏海底模型，模型結(jié)構(gòu)如圖6所示，海水深度在2 242～3 012 m之間變化，設(shè)海水的電阻率為0.3 Ω·m,在不同的地層中分別設(shè)置了兩個低阻異常體，1號異常體電阻率為0.6 Ω·m，2號異常體電阻率為3 Ω·m。MT接收點從左往右依次設(shè)置了21個測點每個測點間距800 m。MT的頻率范圍為0.001～1 Hz，共選取了21個頻點，取對數(shù)等間隔分布。

對上述起伏海底模型進行自適應(yīng)有限元正演計算，得到模型正演響應(yīng)與自適應(yīng)細化網(wǎng)格，表3為起伏海底模型細化次數(shù)，從表3中可以看出，相比于陸地反演模型，海底模型的剖分次數(shù)更多，產(chǎn)生的節(jié)點數(shù)以及單元數(shù)也比陸地模型的數(shù)量多。

在并行機上進行反演試算，如表4所示，反演經(jīng)過21次迭代，耗時91.5 min,最終模型的均方根擬合差(RMS)值為1.151 8，占用內(nèi)存約為3.35 GB。反演結(jié)果如圖7所示，可以看出反演結(jié)果能夠準(zhǔn)確呈現(xiàn)出兩個異常體所處位置、體積形態(tài)等特征，并且具有較高的分辨率，驗證了算法對于起伏海底模型的適用性。圖8為真實模型響應(yīng)數(shù)據(jù)與反演最終模型響應(yīng)數(shù)據(jù)的對比，由圖8可以看出，真實模型的響應(yīng)數(shù)據(jù)與反演模型得到的視電阻率和相位信息高度吻合，體現(xiàn)了反演結(jié)果的真實性和可靠性。

表4 反演參數(shù)Table 4 Inversion parameters

圖7 起伏海底模型反演結(jié)果Fig.7 Inversion result of undulating seabed model

圖8 真實模型響應(yīng)(a)與反演模型響應(yīng)(b)擬斷面Fig.8 Quasi-sectional view of MT response with real model (a) and inverse model (b)

3 實測數(shù)據(jù)反演

3.1 研究區(qū)概況

研究區(qū)域地處西伯利亞板塊和哈薩克斯坦板塊交界地帶，其具體位置處于克拉瑪依后山一帶，鄰近克拉瑪依市區(qū)。綜合該區(qū)域地質(zhì)構(gòu)造及前人研究資料可知，研究區(qū)在地層劃分上屬于北疆—興安嶺地層大區(qū)北疆地層區(qū)，位于哈圖斷裂以南的克拉瑪依地層小區(qū)，以石炭系為主，未見泥盆系。區(qū)內(nèi)的地層分區(qū)在二疊紀(jì)之前表現(xiàn)明顯，早二疊世晚期的卡拉崗組陸相火山巖及火山碎屑巖的出現(xiàn)代表了地層分區(qū)性消失，進入統(tǒng)一的陸內(nèi)演化階段。測區(qū)內(nèi)中新生界地層除了古近系外，其余均有出露。

圖9為研究區(qū)地質(zhì)及測線布置圖，MT測點采用約2 km的間距布置12個測點，測線方向為135°，測線長度約為22 km。從圖9中可以看出達爾布特斷裂橫穿測線而過，在測線上方則有哈圖斷裂橫穿而過，而在測線下方含有大量二疊紀(jì)侵入體。

圖9 研究區(qū)地質(zhì)與測線布置Fig.9 Research area geological and line arrangement map

3.2 數(shù)據(jù)反演與解釋

對實測數(shù)據(jù)進行處理，選取其中42個頻點的觀測數(shù)據(jù)進行反演，頻點選取范圍是1×10-2～1×102Hz。為了減少數(shù)據(jù)的噪聲對反演結(jié)果帶來的偏差，在反演進行之前，對實測數(shù)據(jù)進行了圓滑處理。數(shù)據(jù)在經(jīng)過36次迭代，耗時118.2 min，最終得到反演結(jié)果。

圖10為測線已知地質(zhì)剖面情況和MT數(shù)據(jù)非線性共軛梯度反演結(jié)果的電阻率剖面，圖11為測線淺部剖面Occam二維反演結(jié)果斷面。從圖11中可以看出，淺部剖面的反演結(jié)果具有清晰的電性特征，在剖面中段(10～14 km)存在明顯的近似垂直低阻特征，為達爾布特斷裂，斷裂向西北方向傾，最大深度可達2 km左右，并且有向深部延伸的趨勢，這與地質(zhì)資料相匹配，驗證了將本文算法程序應(yīng)用于MT數(shù)據(jù)得到的反演結(jié)果具有較高的準(zhǔn)確性。

圖10 測線地質(zhì)剖面(a)和MT數(shù)據(jù)非線性共軛梯度反演結(jié)果(b)Fig.10 Geological section (a)and nonlinear conjugate gradient inversion results (b) of MT data

圖11 MT數(shù)據(jù)Occam反演電阻率淺部斷面Fig.11 Shallow section of Occam inversion resistivity of MT data

圖12為測線反演結(jié)果的深部斷面，其中黑色倒三角為測點位，測線方向為135°，結(jié)合淺部數(shù)據(jù)，可以觀察到深度2 km以下的部分在斷裂的北側(cè)有向左下方延伸的低阻帶，并且在0 km附近的深部區(qū)域出現(xiàn)了明顯的低阻異常,與資料中哈圖斷層相吻合。而作為一個大斷裂，達爾布特斷裂兩側(cè)均分布有太勒古拉組中段，在斷裂帶北側(cè)電阻率偏低，南側(cè)偏高；在達爾布特斷裂南側(cè)剖面圖上(14～20 km)存在很高的電阻率特征，為克拉瑪依花崗巖巖體的侵入；在達爾布特斷裂北側(cè)高阻特征不明顯，這和地質(zhì)圖對應(yīng)得很好，進一步驗證了算法的準(zhǔn)確性。

圖12 MT數(shù)據(jù)Occam反演電阻率深部斷面Fig.12 Deep section of Occam inversion resistivity of MT data

4 結(jié)論

對基于非結(jié)構(gòu)自適應(yīng)有限元的帶地形二維MT反演算法進行了應(yīng)用研究，通過二維帶地形陸地及海底模型反演計算，表明該算法可以有效壓制靜態(tài)效應(yīng)，反演對TE數(shù)據(jù)的擬合要比對TM數(shù)據(jù)的擬合更好，同時帶地形模型反演對異常體反演效果較好，對異常體的大小和形態(tài)還原度很高，驗證了該算法對于陸地及海底MT數(shù)據(jù)的共同實用性，證明帶地形MT反演能夠準(zhǔn)確地揭示地下電性結(jié)構(gòu)和探明低阻異常體。

同時，應(yīng)用該方法對克拉瑪依后山區(qū)域?qū)崪y數(shù)據(jù)進行反演，最終得到電阻率結(jié)構(gòu)清晰的反演結(jié)果，并將反演結(jié)果與已知地質(zhì)資料以及非線性共軛梯度法的反演結(jié)果進行了對比驗證，能夠清楚顯示出地質(zhì)資料中表明的達爾布特斷裂等地質(zhì)構(gòu)造。