頻率域2.5維廣域電磁法正反演研究

周峰，張志勇，湯井田，李勇

(1.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西南昌，330013；2.中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙，410083；3.中國(guó)地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所自然資源部地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北廊坊，065000)

廣域電磁法(WFEM)是一種人工源頻率域電磁勘探方法，通過(guò)單分量測(cè)量并用其來(lái)定義具有全區(qū)特性的視電阻率，其特點(diǎn)是采用較小的收發(fā)距來(lái)獲取較大的勘探深度[1?2]。該方法與可控源音頻大地電磁法(CSAMT)[3]相比，摒棄了CSAMT 法在遠(yuǎn)區(qū)測(cè)量的限制、擴(kuò)大了頻率域有源電磁法觀(guān)測(cè)區(qū)范圍。目前，該方法已經(jīng)被廣泛應(yīng)用于頁(yè)巖氣勘察、礦產(chǎn)勘探、地?zé)嵋约八摹⒐こ炭辈榈阮I(lǐng)域[4?5]。

當(dāng)前，廣域電磁法數(shù)據(jù)處理與解釋仍然以一維模型為主，而實(shí)際的野外地電條件和觀(guān)測(cè)形式通常不能視為一維，從而影響該方法資料解譯的準(zhǔn)確性。為此，人們針對(duì)2.5維問(wèn)題有源電磁法開(kāi)展了大量的研究工作。正演方面，為了處理場(chǎng)源奇異性，LEE 等[6?8]開(kāi)展了基于二次場(chǎng)算法有源電磁法正演模擬研究，提高了2.5維有源電磁問(wèn)題正演求解精度，但存在不能處理帶地形復(fù)雜地電模型的不足。偽狄拉克函數(shù)的引入，成功解決了帶地形的2.5 維可控源正演計(jì)算問(wèn)題[9]；另外，海洋可控源電磁法2.5維正演計(jì)算中出現(xiàn)了一種差分算法，有助于場(chǎng)源奇異性問(wèn)題的解決[10]?？傮w上，場(chǎng)源奇異性是制約2.5維可控源電磁法正演求解精度的關(guān)鍵因素之一，該因素同樣會(huì)影響2.5維廣域電磁法正演計(jì)算。反演方面，UNSWORTH 等[11]提出一種子空間反演算法并成功應(yīng)用2.5維可控源電磁數(shù)據(jù)反演，該算法能夠有效提高線(xiàn)性方程組的求解效率。另外，靈敏度矩陣的快速求解也是反演問(wèn)題的關(guān)鍵，基于快速松弛算法的靈敏度矩陣快速近似求解策略被成功應(yīng)用于2.5維可控源電磁法反演，有效提高了反演的效率[12]，該反演算法成功應(yīng)用于實(shí)際數(shù)據(jù)的處理，取得了明顯的效果；此后，底青云[13]采用類(lèi)似的靈敏度矩陣計(jì)算開(kāi)展了復(fù)雜介質(zhì)2.5 維的CSAMT 反演研究。為了進(jìn)一步提高該問(wèn)題反演準(zhǔn)確性和穩(wěn)定性，雷達(dá)[14]結(jié)合OCCAM 反演算法實(shí)現(xiàn)了起伏地形條件下的2.5 維問(wèn)題CSAMT 反演研究，分析了起伏地形的影響；另外，戴世坤等[15]提出了基于場(chǎng)分量的可控源電磁法2.5維問(wèn)題共軛梯度反演算法。綜上所述，雖已存在可控源電磁法2.5維正演與反演方面的研究成果，但基于場(chǎng)分量來(lái)定義廣域視電阻率并以此為視參數(shù)來(lái)進(jìn)行反演研究較少，因此，開(kāi)展廣域電磁法2.5維反演研究對(duì)提高其資料解釋精度及推動(dòng)該方法發(fā)展具有重要的研究意義。

本文作者首先從總場(chǎng)算法出發(fā)，推導(dǎo)2.5維廣域電磁法滿(mǎn)足的微分方程；然后，采用偽Delta 函數(shù)來(lái)降低源帶來(lái)的奇異性，開(kāi)展起伏地形頻率域廣域電磁法2.5 維正演模擬研究。在正演的基礎(chǔ)上，構(gòu)建廣域視電阻率的正則化反演目標(biāo)函數(shù)，采用非線(xiàn)性共軛梯度算法實(shí)現(xiàn)反演目標(biāo)函數(shù)最優(yōu)化求解，以有效避免顯示、存儲(chǔ)靈敏度矩陣，提高反演的求解效率。最后，設(shè)計(jì)層狀介質(zhì)模型驗(yàn)證本文開(kāi)發(fā)的2.5維算法的正確性；并以正演為基礎(chǔ)，分析廣域視電阻率和Cagniard視電阻率的靜態(tài)效應(yīng)影響問(wèn)題；然后，開(kāi)展帶地形與不帶地形的低阻異常體以及組合異常體的地電模型的反演。

1 廣域電磁法2.5維問(wèn)題

1.1 基本理論

時(shí)間因子為eiωt，在準(zhǔn)靜態(tài)條件下，頻率域有源的Maxwell電磁滿(mǎn)足方程[16]為：

式中：E為電場(chǎng)強(qiáng)度；Η為磁場(chǎng)強(qiáng)度；Ms為外部磁流源；Js為外部電流源；=iωu，為阻抗率；ω為角頻率；u為磁導(dǎo)率；=σ-iωε，為導(dǎo)納率；σ為介質(zhì)的電導(dǎo)率；ε為介質(zhì)的介電常數(shù)。

假設(shè)地質(zhì)體走向?yàn)閥方向，對(duì)2.5維問(wèn)題來(lái)說(shuō)，對(duì)構(gòu)造方向y方向進(jìn)行一維傅里葉變換[17]，變換公式為

式中：ky為波數(shù)；F(x,y,z,ω) 為空間域場(chǎng)；(x,ky,z,ω)為通過(guò)構(gòu)造方向y方向進(jìn)行一維傅里葉變換得到的波數(shù)域場(chǎng)值。F(x,y,z,ω)可以表示為電場(chǎng)或磁場(chǎng)，通過(guò)式(3)，式(1)和式(2)可簡(jiǎn)化為：

化簡(jiǎn)式(7)和式(9)，可得：

將式(10)和式(11)代入式(12)和式(13)，可得：

將式(12)和式(13)代入式(10)和式(11)，可化簡(jiǎn)為：

然后，將式(16)和式(17)代入式(5)，可得

同理，可得

式中：k2e=k2y+。式(18)和式(19)為頻率域2.5維問(wèn)題滿(mǎn)足的耦合方程。

1.2 有限元矩陣合成

本文采用前期開(kāi)發(fā)的收縮二叉樹(shù)網(wǎng)格剖分技術(shù)[18?19]對(duì)求解區(qū)域進(jìn)行離散，采用Galerkin 加權(quán)余量法將式(18)和式(19)滿(mǎn)足的微分方程推導(dǎo)為有限元方程，具體公式為

運(yùn)用Galerkin加權(quán)余量法對(duì)式(20)處理，得

式中：Nei為第i單元De上的插值函數(shù)；N為單元總數(shù)。將式(20)代入式(21)，可得同理，式(19)可簡(jiǎn)化為

式(22)和式(23)分別為沿走向電場(chǎng)分量和磁場(chǎng)分量滿(mǎn)足的有限元方程，上述方程通過(guò)單元分析和剛度矩陣的合成得到一組稀疏復(fù)線(xiàn)性方程組：

式中：SE和SH分別為電流和磁流密度的場(chǎng)源積分值，分別為式(24)和式(25)的右端項(xiàng)；KE和KH為電場(chǎng)得到的剛度矩陣；KECH和KHCE均為電場(chǎng)Hy得到的剛度矩陣。通過(guò)分析可知，KECH=-KHCE，為了利用矩陣的對(duì)稱(chēng)性，則將式(24)和式(25)組合成為

式(26)右端源項(xiàng)利用偽Delta 函數(shù)進(jìn)行場(chǎng)源加載，最后采用PARDISO 直接求解器[20]進(jìn)行方程組求解，最后獲取需要的波數(shù)域的電磁場(chǎng)各場(chǎng)分量值。

1.3 源項(xiàng)處理

對(duì)于有源場(chǎng)的計(jì)算，當(dāng)采用總場(chǎng)算法時(shí)場(chǎng)源點(diǎn)是一個(gè)畸變點(diǎn)，在場(chǎng)源點(diǎn)附近電磁場(chǎng)變化很大，在有限的結(jié)點(diǎn)數(shù)量條件下通過(guò)線(xiàn)性、二次插值很難正確表示場(chǎng)源區(qū)附近的場(chǎng)值的變化。MITSUHATA[9]將pseudo-deta 函數(shù)引入電磁計(jì)算問(wèn)題：

式中：τ為偽狄拉克函數(shù)的步長(zhǎng)。

在二維條件下，三維場(chǎng)源在走向方向采用狄拉克函數(shù)，其他2個(gè)方向采用偽狄拉克函數(shù)，則場(chǎng)源可以表示為

除此之外，波束的選擇在2.5維正演計(jì)算中是關(guān)鍵問(wèn)題之一，需要效定波束的取值區(qū)間與波束的數(shù)量。在最有效的區(qū)間內(nèi)，取合適的波束數(shù)量，可以在保證計(jì)算精度的條件下最大程地度提高計(jì)算效率。波束的選擇與網(wǎng)格剖分相關(guān)，對(duì)式(18)，假設(shè)在均勻半空間條件下，只考慮走向y方向存在電流源且則有

MITSUHATA[9]在其采用四邊形網(wǎng)格二次插值的CSEM 2.5 維正演計(jì)算中認(rèn)為ky<1/Δ的有限元計(jì)算結(jié)果是滿(mǎn)足精度要求。在經(jīng)過(guò)大量試算后，在區(qū)間(1.0/xmax,1/Δ)中對(duì)數(shù)等距選擇ky，其中xmax為計(jì)算模型橫向最大尺度。

在此基礎(chǔ)之上，為了得到獲取廣域視電阻率，采用迭代公式進(jìn)行求解。

1)選取廣域視電阻率初值ρ(0)，一般可取對(duì)應(yīng)分量的波區(qū)視電阻率。

2)將ρ(0)代入視電阻率迭代計(jì)算公式(31)，求得第一次迭代后的視電阻率ρ(1)。

3)終止條件判斷。若(ρ(1)-ρ(0))/ρ(0)≤η或者達(dá)到最大迭代次數(shù)(式中η為給定計(jì)算精度)成立，則停止迭代，否則ρ(0)=ρ(1)，令返回步驟2)，直到判斷條件成立為止。

均勻半空間x方向水平電偶極子Ex分量表達(dá)式為

式中：FE-Ex(ikr)=1-3 sin2φ+e-ikr(1+ikr)。測(cè)量收發(fā)距觀(guān)測(cè)點(diǎn)與偶極距正向的夾角φ，同時(shí)，記錄供電電流I和供電偶極的長(zhǎng)度。

2 正則化反演基本理論

反演是電磁法資料解釋的關(guān)鍵環(huán)節(jié)，而正則化因子的選擇與目標(biāo)函數(shù)的優(yōu)化又是正則化反演的2個(gè)重要過(guò)程[21]，本文研究采用“L-curve”曲線(xiàn)方法[22]，實(shí)現(xiàn)了正則化因子求解；另外，采用非線(xiàn)性共軛梯度算法實(shí)現(xiàn)平滑約束穩(wěn)定因子的反演目標(biāo)函數(shù)求解，提高反演過(guò)程的穩(wěn)定性與計(jì)算效率。

采用正則化最小化目標(biāo)函數(shù)：

式中：m為模型參數(shù)；dobs為觀(guān)測(cè)數(shù)據(jù)；μ為正則化因子；φd(m,dobs)為數(shù)據(jù)擬合函數(shù)；為數(shù)據(jù)權(quán)重；G為正演算子；φm(m)為模型約束函數(shù)；Wm為模型權(quán)重矩陣；mapr為先驗(yàn)?zāi)Ｐ秃瘮?shù)；φ(m,dobs)為反演目標(biāo)函數(shù)。

反問(wèn)題就是對(duì)目標(biāo)函數(shù)求極值。本文采用非線(xiàn)性共軛梯度法利用迭代過(guò)程來(lái)實(shí)現(xiàn)目標(biāo)函數(shù)極值點(diǎn)求取。需要對(duì)目標(biāo)函數(shù)(32)求一階梯度，其表達(dá)式為

將式(32)擴(kuò)展為

1)令i=1，選擇初始模型mi，同時(shí)求解目標(biāo)函數(shù)的梯度gi=-?φ(mi,dobs)。

2)令ui=M-1i gi.

3)選擇1個(gè)αi，使得φm(m+αi gi)極小。

4)選擇下一次更新模型mi+1=mi+αi gi以及目標(biāo)函數(shù)的梯度gi+1=-?φm(mi+1,dobs)。

5)當(dāng)|gi+1|的值做夠小或者達(dá)到設(shè)定的參數(shù)，即停止執(zhí)行；否則執(zhí)行步驟6)。

6)令βi+1=

7)令ui+1=M-1i+1gi+1+βi+1ui，同時(shí)，令i=i+1后返回步驟3)。

從以上計(jì)算過(guò)程可以知道，非線(xiàn)性共軛梯度算法只需要計(jì)算目標(biāo)函數(shù)以及其梯度，不需要過(guò)多的矩陣的分解等運(yùn)算，只需幾個(gè)存儲(chǔ)有限向量，因此，該算法適用于大規(guī)模問(wèn)題的求解。

3 廣域電磁2.5維正演算例分析

3.1 算法驗(yàn)證

為了驗(yàn)證2.5維廣域電磁法程序的正確性，設(shè)計(jì)了如圖1所示的地電模型，具體參數(shù)描述為：2層地電模型電阻率分別為100 Ω·m 與10 Ω·m，第一層厚度為512 m，在坐標(biāo)原點(diǎn)(0，0，0)m處設(shè)置1 個(gè)沿y軸正方向激發(fā)源，源的長(zhǎng)度為1 km，沿著x方向設(shè)置3個(gè)觀(guān)測(cè)點(diǎn)，觀(guān)測(cè)點(diǎn)的坐標(biāo)分別為Site-A(1 020，0，0) m，Site-B(2 020，0，0) m，Site-C(5 020，0，0)m。分別計(jì)算3 個(gè)觀(guān)測(cè)點(diǎn)的Cagniard電阻率及Ey，Hx和Hz這3 個(gè)分量廣域視電阻率，計(jì)算結(jié)果分別如圖2～4所示。從圖2～4可知：2.5D視電阻率曲線(xiàn)與1D 結(jié)果曲線(xiàn)形態(tài)吻合[24]，Site-A曲線(xiàn)中的大部分頻點(diǎn)位于近區(qū)與過(guò)渡帶，而Cagniard視電阻率曲線(xiàn)在低頻呈現(xiàn)連續(xù)上升；廣域視電阻率曲線(xiàn)呈現(xiàn)出D 型曲線(xiàn)特征，與實(shí)際模型吻合。Site-B廣域視電阻率表現(xiàn)出明顯的D型曲線(xiàn)特征，Cagniard 視電阻率在4 Hz 左右出現(xiàn)過(guò)渡帶低阻；Site-C 總體表現(xiàn)出遠(yuǎn)區(qū)曲線(xiàn)特征。觀(guān)測(cè)點(diǎn)Site-A和Site-B的Hx分量的視電阻率曲線(xiàn)在低頻差別較大，原因是Hx分量在近區(qū)與過(guò)渡區(qū)含有的地下介質(zhì)電阻率信息少，在求取廣域視電阻率過(guò)程中，求解非線(xiàn)性方程小的誤差將引起較大的視電阻率差異。以上結(jié)果分析表明，本文開(kāi)發(fā)的2.5維正演程序正確且有效。

圖1 2層介質(zhì)模型Fig.1 Two layer media model

圖2 Site-A觀(guān)測(cè)點(diǎn)2.5D與1D計(jì)算結(jié)果對(duì)比Fig.2 Comparison of 2.5D and 1D calculation results at observation point Site-A

圖3 Site-B 觀(guān)測(cè)點(diǎn)2.5D與D計(jì)算結(jié)果對(duì)比Fig.3 Comparison of 2.5D and 1D calculation results at observation point Site-B

圖4 Site-C 觀(guān)測(cè)點(diǎn)2.5D與1D計(jì)算結(jié)果Fig.4 Comparison of 2.5D and 1D calculation results at observation point Site-C

3.2 靜態(tài)效應(yīng)分析

為了模擬地表不均勻體對(duì)廣域視電阻率的影響，設(shè)計(jì)如圖5所示的地電模型，其中在均勻大地的電阻率為100 Ω·m 中放置了6 個(gè)矩形塊狀異常體，其中4個(gè)上頂埋深為16 m，電阻率從左到右分別為10，1 000，10 和1 000 Ω·m；此外，在地下存在2個(gè)規(guī)模較大的塊狀異常體，塊狀異常體的頂部距離地表埋深為128 m，電導(dǎo)率為10 Ω·m。分別計(jì)算x與y方向的電性場(chǎng)源激勵(lì)下的Cagniard 電阻率及電場(chǎng)分量定義的廣域視電阻率響應(yīng)，計(jì)算結(jié)果如圖6所示。

圖5 靜態(tài)效應(yīng)地電模型Fig.5 Geoelectric model of static effect

圖6 靜態(tài)效應(yīng)地電模型計(jì)算結(jié)果Fig.6 Calculation results of static effect geoelectric model

由圖6可知：在垂直走向方向Idx場(chǎng)源激勵(lì)下，淺部異常體呈現(xiàn)出條帶狀，對(duì)視電阻率剖面影響較大；而沿走向方向Idy場(chǎng)源激勵(lì)下，淺部異常體產(chǎn)生的靜態(tài)效應(yīng)比垂直走向方向激勵(lì)源產(chǎn)生的靜態(tài)效應(yīng)小，在沿走向方向場(chǎng)源激勵(lì)下，淺部異常體呈現(xiàn)出局部圈閉的異?，F(xiàn)象；在與背景電導(dǎo)率存在10 倍差異的條件下，低阻的異?，F(xiàn)象比高阻的明顯；淺部異常體產(chǎn)生的靜態(tài)效應(yīng)對(duì)Cagniard電阻率影響遠(yuǎn)比單分量定義的廣域電阻率的影響大，尤其是當(dāng)場(chǎng)源方向垂直走向方向時(shí)，較深的異常體嚴(yán)重畸變。

3.3 廣域電磁法地形影響

地形能夠使電磁法數(shù)據(jù)產(chǎn)生嚴(yán)重畸變。為了進(jìn)一步分析廣域電磁法和傳統(tǒng)可控源電磁法在視參數(shù)定義受地形影響程度，設(shè)計(jì)均勻半空間存在地形的模型進(jìn)行試算研究，模型參數(shù)如圖7所示，在均勻大地(100 Ω·m)上存在1 個(gè)截面為梯形的地表地形，分別計(jì)算2個(gè)方向場(chǎng)源激勵(lì)下Cagniard電阻及電場(chǎng)分量定義廣域視電阻率地形響應(yīng)，計(jì)算結(jié)果分別如圖8和圖9所示。

圖7 地形模型Fig.7 Terrain model

圖8 x方向場(chǎng)源的計(jì)算結(jié)果Fig.8 Calculation results of x direction horizental electric dipole source

由圖8和圖9可見(jiàn)：地形對(duì)x方向場(chǎng)源的激勵(lì)明顯比對(duì)y方向場(chǎng)源的激勵(lì)大，且對(duì)Cagniard電阻的影響比對(duì)廣域電阻率的影響大，近區(qū)與過(guò)渡帶數(shù)據(jù)畸變大。

圖9 y方向場(chǎng)源的計(jì)算計(jì)算結(jié)果Fig.9 Calculation results of y direction horizental electric dipole source

3.4 廣域電磁2.5維反演算例

3.4.1 單個(gè)異常體反演算例

建立如圖10所示模型，分別開(kāi)展不帶地形、帶地形(圖10中虛線(xiàn))條件下存在單個(gè)異常體的反演。將正演計(jì)算視電阻率加入1%隨機(jī)噪聲作為反演輸入數(shù)據(jù)。

圖10 單個(gè)異常體反演模型示意圖Fig.10 Schematic diagram of a single anomaly body inversion model

圖11所示為單個(gè)異常體x方向電偶源的Ex分量定義廣域視電阻率經(jīng)12 次反演迭代得到的計(jì)算結(jié)果。由圖11可知：帶地形與不帶地形條件下均能夠給出異常體大致位置，異常體的邊界平滑過(guò)渡，反演結(jié)果準(zhǔn)確。圖12所示為采用Cagniard 視電阻率作為反演數(shù)據(jù)得到的反演結(jié)果，同樣顯示了較高的準(zhǔn)確度。以上結(jié)果表明，無(wú)論是廣域視電阻率的反演還是Cagniard電阻率的反演，深部電阻率的收斂都較慢，低阻體向深部且偏向場(chǎng)源方向延伸。

圖11 單個(gè)異常體x方向電偶源Ex分量定義廣域視電阻率反演結(jié)果Fig.11 Inversion results of Ex component of wide field apparent resistivity for single anomaly body in x-directional horizontal electric dipole source

圖12 單個(gè)異常體x方向電偶源Cagniard分量定義視電阻率反演結(jié)果Fig.12 Inversion results of Cagniard apparent resistivity for single anomaly body in x-directional horizontal electric dipole source

3.4.2 組合異常體反演算例

建立如圖13所示模型，分別開(kāi)展不帶地形、帶地形(圖13中虛線(xiàn))條件下存2 個(gè)異常體的反演。將反演數(shù)據(jù)加入1%隨機(jī)噪聲以正演計(jì)算視電阻率。平滑模型約束下，應(yīng)用非線(xiàn)性共軛梯度求解的反演結(jié)果如圖14所示。組合異常體反演結(jié)果表明，x方向電偶源Ex分量定義廣域視電阻率反演得到的低阻異常體擬合程度明顯比高阻的高，高阻反演得到的電阻率差異明顯。同樣地，無(wú)論是廣域視電阻率反演還是Cagniard電阻率反演，深部電阻率的收斂速率都較慢，表現(xiàn)為低阻體向深部且偏向場(chǎng)源方向延伸。

圖13 組合異常體反演模型示意圖Fig.13 Schematic diagram of a combination anomaly body inversion model

圖14 2個(gè)異常體x方向電偶源Ex分量廣域視電阻率反演結(jié)果Fig.14 Inversion results of Ex component of wide field apparent resistivity for combination anomaly body in x-directional horizontal electric dipole source

4 結(jié)論

1)在二維地電條件下，淺部異常體引起廣域電磁的靜態(tài)效應(yīng)?？傮w而言，Ex形式定義的廣域電阻率受靜態(tài)效應(yīng)影響比Cagniard 電阻率的影響小，且赤道裝置得到的廣域視電阻率受靜態(tài)效應(yīng)的影響比軸向裝置得到的廣域視電阻率所受的影響小。另外，地形對(duì)x方向場(chǎng)源激勵(lì)明顯比y方向場(chǎng)源激勵(lì)大，且對(duì)Cagniard電阻影響比廣域電阻率的影響大，近區(qū)與過(guò)渡帶數(shù)據(jù)畸變大。

2)本文開(kāi)發(fā)的廣域電磁法2.5維反演結(jié)果能夠很好地反映異常所在的位置，反演結(jié)果總體準(zhǔn)確，不足之處在于底部邊界收斂速率較慢。

3)為了進(jìn)一步測(cè)試本文算法的可靠性，下一步將開(kāi)展實(shí)際數(shù)據(jù)測(cè)試，以便對(duì)本文開(kāi)發(fā)的算法提出相應(yīng)的改進(jìn)措施，完善算法。