基于交叉梯度結(jié)構(gòu)約束的大地電磁主軸各向異性并行三維反演

王堃鵬， 王緒本， 曹輝， 藍(lán)星， 段長(zhǎng)生， 羅威， 原健龍

1 成都理工大學(xué)地球物理學(xué)院， 成都 610059 2 油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))， 成都 610059 3 中國(guó)地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心, 成都 610081 4 四川省冶勘設(shè)計(jì)集團(tuán)有限公司， 成都 610051 5 贛中南地質(zhì)礦產(chǎn)勘查研究院， 南昌 330029

0 引言

大地電磁(MT)是目前廣泛使用的一類重要電磁勘探方法，被廣泛應(yīng)用于油氣礦產(chǎn)勘探及地球動(dòng)力學(xué)領(lǐng)域的研究(Zhang et al.，2014；Heinson et al.，2006；王緒本等，2013)．然而地下介質(zhì)電性各向異性的客觀存在，以及它對(duì)MT觀測(cè)數(shù)據(jù)的影響，使其受到越來越多的關(guān)注(Wannamaker，2005；劉云鶴等，2018)．許多學(xué)者從正演出發(fā)，不斷的研究它對(duì)MT觀測(cè)數(shù)據(jù)的影響規(guī)律．Pek和Verner(1997)使用有限差分法研究了二維MT的各向異性正演. Li和Pek(2008)利用自適應(yīng)有限單元法實(shí)現(xiàn)了一般各向異性的MT二維正演．Kirkby等(2015)嘗試使用各向異性正演研究斷裂分布. 嵇艷鞠等(2016)實(shí)現(xiàn)了一種基于無網(wǎng)格法的各向異性正演．曹曉月等(2018)研究了考慮地形的自適應(yīng)各向異性正演．王宇等(2019)使用有限體積法研究了MT二維各向異性．Guo等(2018)實(shí)現(xiàn)了模塊化的各向異性MT二維正演.Xiao等(2018)利用有限單元法研究了MT三維各向異性．基于這些學(xué)者大量的研究，MT觀測(cè)數(shù)據(jù)受電性各向異性影響的規(guī)律正在被逐漸認(rèn)識(shí)，尤其是對(duì)地質(zhì)解釋可能造成的干擾已經(jīng)受到了廣泛關(guān)注．

同時(shí)，一些學(xué)者嘗試在MT反演中直接反演電性各向異性參數(shù)．Yin(2003)與Pek和Santos(2006)在一維介質(zhì)模型中加入了各向異性參數(shù)反演．Schmoldt和Jones(2013)研究了一種新的二維MT各向異性反演方法．Montahaei和Oskooi(2014)使用神經(jīng)網(wǎng)絡(luò)法實(shí)現(xiàn)了MT二維各向異性反演，并對(duì)南美洲多條MT測(cè)線的實(shí)測(cè)數(shù)據(jù)進(jìn)行了各向同性與各向異性反演研究與分析．Martí(2014)總結(jié)了關(guān)于MT各向異性的正演和反演研究現(xiàn)狀．Cao等(2018)實(shí)現(xiàn)了MT三維主軸各向異性反演．Johansen等(2019)利用MARE2DEM(Key，2016)對(duì)一條海底測(cè)線數(shù)據(jù)實(shí)施了CSEM (Controlled-Source Electromagnetic Method，可控源電磁法)與MT的主軸各向異性二維反演．上述研究在考慮各向異性參數(shù)時(shí)，并沒有充分研究結(jié)構(gòu)約束對(duì)MT各向異性反演的影響．

加入各向異性參數(shù)的反演，從理論而言可以進(jìn)一步減少各向同性反演產(chǎn)生的假異常，但受限于觀測(cè)數(shù)據(jù)的不足，實(shí)際上可能加劇多參數(shù)反演的不確定性．因此，直接反演各向異性介質(zhì)的所有參數(shù)非常困難．現(xiàn)有的常規(guī)思路首先是簡(jiǎn)化各向異性介質(zhì)，比如VTI (Vertical Transversely Isotropy，垂直對(duì)稱軸的橫向各向同性)介質(zhì)、主軸各向異性介質(zhì)等(羅鳴等，2016；Key，2016；彭榮華等，2019)．除此之外，為了進(jìn)一步縮小模型空間，還可以在各向異性反演中強(qiáng)加對(duì)各向異性參數(shù)的約束．Pain等(2003)提出了一種各向異性罰函數(shù)來約束各向異性反演．Meju等(2018，2019)使用交叉梯度法，對(duì)海洋可控源電磁法三維各向異性反演實(shí)施了進(jìn)一步的結(jié)構(gòu)約束．

根據(jù)上述研究現(xiàn)狀，我們將在已有的三維MT各向異性反演基礎(chǔ)上(Cao et al.，2018)，選用交叉梯度項(xiàng)對(duì)各向異性參數(shù)進(jìn)行結(jié)構(gòu)約束，用于約束的結(jié)構(gòu)信息可以來自人工地震、天然地震或重磁等方法．本文的正演將大地電磁場(chǎng)分為一次場(chǎng)和二次場(chǎng)，并將各向異性介質(zhì)簡(jiǎn)化為僅考慮X、Y、Z方向電阻率的主軸各向異性，簡(jiǎn)化后的X、Y、Z三個(gè)方向電阻率可以自由的受輸入結(jié)構(gòu)的約束，本文的結(jié)構(gòu)約束方案同樣能用于聯(lián)合反演，不僅局限于實(shí)現(xiàn)結(jié)構(gòu)約束．

1 基于二次場(chǎng)的MT主軸各向異性三維正演

本文將大地電磁三維正演進(jìn)行背景場(chǎng)和二次場(chǎng)分解，取時(shí)諧因子為e-iω t，且忽略位移電流，可得到如下關(guān)于電場(chǎng)的似穩(wěn)態(tài)偏微分方程：

(1)

(2)

其中，ω為圓頻率，μ0為真空磁導(dǎo)率，σ為介質(zhì)電導(dǎo)率，Js表示源項(xiàng)．將(1)、(2)式做簡(jiǎn)單代換，可得到電場(chǎng)的二階偏微分方程：

(3)

將(3)式減去背景場(chǎng)，可獲得二次場(chǎng)的二階偏微分方程：

(4)

式中σp為背景電導(dǎo)率，Ep為背景電場(chǎng)．這里的背景電場(chǎng)利用一維層狀介質(zhì)的大地電磁方程計(jì)算獲得．為了降低各向異性反演的多解性和不穩(wěn)定性，本文將(4)式中的電導(dǎo)率參數(shù)σ替換為如下主軸各向異性參數(shù)：

(5)

則(4)式在三個(gè)方向的微分方程如下：

(6)

(7)

(8)

針對(duì)(6)—(8)三式，本文采用Yee格式及交錯(cuò)網(wǎng)格有限差分法進(jìn)行離散，最終的方程采用QMR (Quasi Minimal Residual，擬最小殘差)算法求解．關(guān)于主軸各向異性三維正演的更多細(xì)節(jié)可參考作者已發(fā)表的論文(Cao et al., 2018)，本節(jié)不再贅述．

2 基于交叉梯度約束的MT三維主軸各向異性反演

2.1 各向異性目標(biāo)函數(shù)

本文采用LBFGS法進(jìn)行反演迭代計(jì)算，由于存在三個(gè)方向的參數(shù)并且觀測(cè)數(shù)據(jù)為復(fù)數(shù)阻抗，本文設(shè)目標(biāo)函數(shù)為：

φ=φd+λxφm_x+λyφm_y+λzφm_z+λcgxφcgx

+λcgyφcgy+λcgzφcgz,

(9)

本文采用仿射線性參數(shù)變換(Egbert and Kelbert，2012)對(duì)真實(shí)模型參數(shù)進(jìn)行如下轉(zhuǎn)換：

(10)

(11)

在參數(shù)變換策略下，目標(biāo)函數(shù)(9)修改為如下形式：

(12)

2.2 交叉梯度項(xiàng)計(jì)算

將(13)式按照三個(gè)方向進(jìn)一步展開有：

(14)

參考閆政文等(2020)，在式(14)中，tx_x、ty_x、tz_x三個(gè)向量的任意元素可以表述為：

(15)

其中，mx與mvx分別為X方向的某個(gè)電阻率參數(shù)和約束模型參數(shù)．對(duì)于某個(gè)具體的tx_x(i,j,z)計(jì)算，這里采用圖1來說明．圖1展示了一個(gè)從三維離散體抽取的二維離散面，該面是Y軸與Z軸組成的平面，以中心差分為準(zhǔn)則，tx_x(i,j,z)可以描述為：

圖1 從三維離散體抽取的二維離散面Fig.1 An extracted 2D discrete surface from the 3D discrete body

將式(16)整理并以矩陣形式表述如下：

(17)

其中，各元素具體形式如下：

a33=-a11-a22-a44-a55.

(18)

在圖1中的內(nèi)部單元循環(huán)，可以完成tx_x向量所有元素的計(jì)算．由于mvx是輸入的約束模型參數(shù)，在反演中固定，實(shí)際上最終的tx_x向量是由式(17)的矩陣?yán)奂有纬?，最終拓展到整個(gè)模型單元可以表述為：

tx_x=Wx_cg_xmx.

(19)

同樣的，ty_x、tz_x可以按照上述方式在另外兩個(gè)方向離散，最終以矩陣形式表述為：

ty_x=Wy_cg_xmx,

tz_x=Wz_cg_xmx.

(20)

(21)

(22)

(23)

2.3 各向異性梯度計(jì)算

對(duì)于LBFGS來說，梯度計(jì)算是最重要的步驟之一，根據(jù)目標(biāo)函數(shù)(12)，以及式(21)(22)(23)，對(duì)轉(zhuǎn)換參數(shù)求導(dǎo)得到以下梯度表達(dá)式：

(24)

式中，數(shù)據(jù)項(xiàng)的目標(biāo)函數(shù)具體形式如下：

(25)

其中，Re代表取實(shí)部，Jx、Jy、Jz是關(guān)于真實(shí)電阻率參數(shù)的靈敏度矩陣．對(duì)于LBFGS來說，該部分的求解技巧與NLCG (Non-Linear Conjugate Gradient，非線性共軛梯度)一致，即可通過一系列“擬正演”(Newman and Alumbaugh，2000)避免直接求解靈敏度矩陣，提高計(jì)算效率和節(jié)省大量?jī)?nèi)存．

關(guān)于式(24)后面的模型項(xiàng)和交叉梯度項(xiàng)的梯度，具體如下：

(26)

從式(25)和式(26)可以看到，要完成基于交叉梯度結(jié)構(gòu)約束的MT主軸各向異性反演，梯度計(jì)算比較復(fù)雜．

2.4 LBFGS反演框架

本文的LBFGS反演流程，具體如下：

(1)選擇一個(gè)反演初始模型mi，設(shè)定三個(gè)方向的正則化因子和交叉梯度權(quán)重因子；

(3)尋找一個(gè)相對(duì)合適的更新步長(zhǎng)，以最小化φ(mi+αiui)；

(4)計(jì)算更新模型mi+1=mi+αiui，判斷反演是否滿足停止條件，反之則開始第(5)步；

(27)

其中，si-1為模型差向量，yi-1為梯度差向量.

3 合成算例分析

3.1 各向同性模型無約束與約束反演對(duì)比

本節(jié)首先使用2.2節(jié)的交叉梯度約束項(xiàng)進(jìn)行各向同性反演試驗(yàn)，以驗(yàn)證本文的約束方案首先能應(yīng)用于各向同性介質(zhì)的結(jié)構(gòu)約束反演.2.2節(jié)給出的是主軸各向異性約束反演，但該方案同樣可應(yīng)用于各向同性反演，僅僅只需要將其他兩個(gè)方向的交叉梯度約束項(xiàng)去掉即可.

圖2為測(cè)點(diǎn)分布及模型俯視圖，測(cè)點(diǎn)在X和Y方向間距160 m，共2500個(gè)測(cè)點(diǎn).模型為一個(gè)低阻體(10 Ωm)與高阻體(1000 Ωm)的組合模型，并在空間上接觸面較多.模型的空間分布進(jìn)一步在圖3第一列顯示，背景電阻率為100 Ωm.為了進(jìn)行約束反演研究，這里進(jìn)一步建立了如圖4所示的速度模型，背景值為4000 m·s-1，低速異常體為2000 m·s-1，高速異常體為6000 m·s-1.該速度異常體與電阻率的結(jié)構(gòu)分布一致.

為了比較加入速度結(jié)構(gòu)約束與沒有速度結(jié)構(gòu)約束的反演結(jié)果，使用圖4建立的速度模型作為約束模型輸入交叉梯度項(xiàng)，交叉梯度權(quán)重因子設(shè)為1.0×107，正則化因子為0.1，反演初始模型為50 Ωm均勻半空間模型.經(jīng)過28次迭代后，反演擬合差從7.70下降到0.996，最終的反演結(jié)果如圖3第三列所示.從圖3的結(jié)果可以看到，使用交叉梯度進(jìn)行結(jié)構(gòu)約束后，高阻體的結(jié)構(gòu)相比于無結(jié)構(gòu)約束反演更加清晰，這說明本文基于交叉梯度的結(jié)構(gòu)約束方案，對(duì)各向同性反演是有效的.

圖2 測(cè)點(diǎn)分布及模型俯視圖Fig.2 The measuring sites distribution and top view

圖3 第一列為真實(shí)模型，第二列為無交叉梯度約束反演結(jié)果，第三列為交叉梯度約束反演結(jié)果Fig.3 The first column is the true model．The second column is the inversion results without cross-gradient constraint. The third column is the inversion results with cross-gradient constraint

最后，我們給出了無約束反演和約束反演的步長(zhǎng)曲線及擬合差收斂曲線(圖5)，可以看到本文采用的LBFGS法進(jìn)行無約束與約束反演均非常穩(wěn)定.此外還需要說明的是，關(guān)于本文的交叉梯度因子的選擇，我們經(jīng)過大量的實(shí)驗(yàn)后，有一個(gè)經(jīng)驗(yàn)策略，即約束反演二次迭代的總目標(biāo)函數(shù)值，相較于無約束反演二次迭代的總目標(biāo)函數(shù)值提高70%～80%，約束效果及反演效率相對(duì)較好.

3.2 各向同性與主軸各向異性反演對(duì)比

3.2.1 各向同性與各向異性無約束反演對(duì)比

為了說明在主軸各向異性明顯時(shí)，各向同性反演的不足，這里首先使用無約束的各向同性與各向異性大地電磁三維反演程序?qū)D6模型進(jìn)行反演，最終的反演結(jié)果如圖7所示．圖7第一列是各向同性反演，第二、三、四列為各向異性反演．從圖7的結(jié)果可以明顯看出，當(dāng)?shù)貙咏橘|(zhì)存在明顯主軸各向異性時(shí)，圖7第一列顯示出了極強(qiáng)的假構(gòu)造分布，對(duì)于實(shí)際情況而言這種結(jié)果會(huì)嚴(yán)重誤導(dǎo)構(gòu)造解釋．

圖4 用于交叉梯度約束的速度模型Fig.4 The velocity model for cross-gradient constraint

圖6 第一列為X方向電阻率真實(shí)模型，第二列為Y方向電阻率真實(shí)模型，第三列為Z方向電阻率真實(shí)模型Fig.6 The first column is the true model of X direction resistivity．The second column is the true model of Y direction resistivity. The third column is the true model of Z direction resistivity

圖7 第一列為各向同性反演結(jié)果，第二列為X方向電阻率反演結(jié)果，第三列為Y方向電阻率反演結(jié)果，第四列為Z方向電阻率反演結(jié)果Fig.7 The first column is the isotropic inversion result．The second column is the inversion result of X direction resistivity. The third column is the inversion result of Y direction resistivity. The fourth column is the inversion result of Z direction resistivity

圖7第二列和第三列的X與Y方向電阻率反演結(jié)果，沒有出現(xiàn)第一列強(qiáng)烈的假構(gòu)造問題.由于X方向電阻率真實(shí)值(10 Ωm)相較于Y方向電阻率(30 Ωm)更低，而大地電磁對(duì)低阻更靈敏，因此X方向電阻率恢復(fù)情況要好于Y方向電阻率．對(duì)于第四列的Z方向電阻率，從顏色分布上有一定的擾動(dòng)，但絕大多數(shù)值還是與初始模型差不多，這說明大地電磁法來自高空的平面波場(chǎng)源對(duì)地層的垂直電阻率不靈敏(Cao et al., 2018)．

整個(gè)反演的步長(zhǎng)曲線及收斂曲線如圖8所示，各向同性與各向異性反演的最終相對(duì)擬合差分別為0.98，0.93，二者均滿足收斂條件．但從迭代次數(shù)可以發(fā)現(xiàn)，各向同性反演迭代次數(shù)遠(yuǎn)大于各向異性，前者迭代了77次，后者僅僅迭代了11次．因此，盡管圖8的步長(zhǎng)與收斂曲線顯示各向同性與各向異性反演都是穩(wěn)定收斂的，但結(jié)合圖7的反演結(jié)果來看，圖8的曲線再一次說明若地層存在顯著主軸各向異性效應(yīng)時(shí)，簡(jiǎn)單的各向同性反演存在較大的不足．

圖8 (a)為步長(zhǎng)曲線，(b)為rms曲線Fig.8 (a) is the step length curves. (b) is the rms curves

3.2.2 各向同性與各向異性約束反演對(duì)比

在本節(jié)，我們將繼續(xù)使用圖6模型的阻抗數(shù)據(jù)進(jìn)行反演，討論各向同性與各向異性反演在交叉梯度約束下的不同反演結(jié)果．首先建立圖9所示的速度模型，背景速度值為6000 m·s-1，低速異常體速度為2000 m·s-1，該速度異常體與電阻率異常的結(jié)構(gòu)分布一致.在各向同性約束與主軸各向異性約束反演中，所有的交叉梯度權(quán)重因子均為1.1×106，此外我們對(duì)各向同性以及各向異性的X、Y、Z方向的電阻率均使用圖9的速度結(jié)構(gòu)進(jìn)行約束．

圖9 用于交叉梯度約束的速度模型Fig.9 The velocity model for cross-gradient constraint

各向同性與各向異性約束反演的最終結(jié)果如圖10所示，第一列為各向同性的約束反演，從這個(gè)結(jié)果可以看出，交叉梯度項(xiàng)對(duì)最外層邊界的約束相較于圖7第一列有明顯的改善，但沒有改變低阻體內(nèi)部假構(gòu)造的分布，這說明交叉梯度也無法改變各向異性介質(zhì)對(duì)各向同性反演造成的干擾．

圖10第二和第三列分別為X及Y方向電阻率的約束反演結(jié)果，相較于第一列而言仍然比較符合真實(shí)模型的分布.與圖7第二、三列相比，邊界有一定的改善，但效果不明顯，主要原因是模型較為簡(jiǎn)單，我們將在下一節(jié)建立相對(duì)復(fù)雜的模型探討交叉梯度對(duì)主軸各向異性反演的影響．

圖11展示了約束下的各向同性與各向異性反演步長(zhǎng)曲線與收斂曲線，從該結(jié)果可以看到，二者都體現(xiàn)了非常好的穩(wěn)定性，但各向同性反演迭代次數(shù)更多，且很難達(dá)到收斂條件．圖10與圖11的反演結(jié)果說明當(dāng)介質(zhì)存在明顯主軸各向異性時(shí)，對(duì)各向同性反演結(jié)果產(chǎn)生極大的干擾．

圖10 第一列為各向同性反演結(jié)果，第二列為X方向電阻率反演結(jié)果，第三列為Y方向電阻率反演結(jié)果，第四列為Z方向電阻率反演結(jié)果Fig.10 The first column is the isotropic inversion result．The second column is the inversion result of X direction resistivity. The third column is the inversion result of Y direction resistivity. The fourth column is the inversion result of Z direction resistivity

圖11 (a)為步長(zhǎng)曲線，(b)為rms曲線Fig.11 (a) is the step length curves. (b) is the rms curves

3.3 復(fù)雜各向異性模型反演

在3.2節(jié)，我們證明了在主軸各向異性效應(yīng)顯著時(shí)，強(qiáng)制使用各向同性的無約束和約束反演都會(huì)帶來嚴(yán)重的假構(gòu)造分布．由于模型較為簡(jiǎn)單，各向異性約束反演并沒有特別好的體現(xiàn)出交叉梯度的優(yōu)勢(shì)和特點(diǎn)，本節(jié)我們將建立一個(gè)形態(tài)較為復(fù)雜的模型，探討交叉梯度項(xiàng)對(duì)大地電磁主軸各向異性反演的約束效果．

之前的實(shí)驗(yàn)已經(jīng)證明大地電磁對(duì)X及Y方向電阻率最為靈敏，而在實(shí)際地質(zhì)情況中，褶皺是最有可能出現(xiàn)X與Y方向電阻率不同的．由于沉積作用形成的層理，在褶皺內(nèi)部垂直層理和平行層理的電阻率是不同的，若出現(xiàn)如圖12a所示的褶皺，則如圖12b所示極有可能造成X與Y方向的宏觀電阻率出現(xiàn)明顯差異．

這里建立一個(gè)在地表有剝蝕的褶皺構(gòu)造，共兩部分組成．左邊褶皺出露地表，右邊褶皺頂面埋深150 m．設(shè)左右褶皺的X方向電阻率均為15 Ωm，Y方向電阻率40 Ωm，Z方向電阻率50 Ωm，背景電阻率200 Ωm，最終的模型示意圖見圖13第一列至第三列．本節(jié)計(jì)算頻率8個(gè)，分別為100，80，50，20，10，5，0.5，0.1 Hz．正演數(shù)據(jù)添加的噪聲水平，測(cè)點(diǎn)分布等信息與3.2節(jié)一致．在進(jìn)行各向異性反演前，我們?nèi)匀幌冉o出一個(gè)各向同性反演結(jié)果(圖13第四列)，初始模型為150 Ωm均勻半空間.從這個(gè)結(jié)果可以看到，左右褶皺的形態(tài)均沒有正確恢復(fù)，出現(xiàn)了非常嚴(yán)重的假構(gòu)造，因此實(shí)際情況的三維反演可能還需要更多的考慮.

圖12 褶皺示意圖Fig.12 The schematic diagram of folds

現(xiàn)在開始進(jìn)行無約束與約束的主軸各向異性三維反演，讀入的約束模型如圖14所示，背景速度值為6000 m·s-1，低速異常體速度為3000 m·s-1，該速度異常體與電阻率異常的結(jié)構(gòu)分布一致.所有的交叉梯度權(quán)重因子均為1.5×106，三個(gè)方向電阻率均使用圖14的速度模型約束.反演初始模型為150 Ωm，最終的三個(gè)方向電阻率反演結(jié)果如圖15、17、19所示.針對(duì)本節(jié)的各向異性約束反演，我們并沒有嘗試突出某個(gè)方向的參數(shù)，因此每個(gè)方向的正則化因子都是一樣的，交叉梯度權(quán)重因子也是如此．此外，交叉梯度權(quán)重因子的選擇策略與3.1節(jié)中的方案一致.

首先來看圖15所示的X方向電阻率反演結(jié)果，圖15b和c在中深部的邊界形態(tài)有明顯差異，圖15c的形態(tài)更接近真實(shí)電阻率的分布.為了更直觀的體現(xiàn)反演差異，這里用圖15b和c的結(jié)果分別與速度模型進(jìn)行每個(gè)單元的交叉梯度項(xiàng)φcgx的計(jì)算，最終結(jié)果如圖16所示.從圖16a可以看到，無約束反演的邊界與速度模型的邊界有明顯的交叉梯度值，這說明二者在邊界形態(tài)上的變化是不一致的.在圖16b中，基于結(jié)構(gòu)約束后的X方向電阻率與速度模型的交叉梯度值明顯降低，這個(gè)結(jié)果再一次驗(yàn)證了本文的交叉梯度約束是有效的.

圖13 第一列為X方向真實(shí)電阻率，第二列為Y方向真實(shí)電阻率，第三列為Z方向真實(shí)電阻率，第四列為各向同性反演結(jié)果Fig.13 The first column is the true model of X direction resistivity．The second column is the true model of Y direction resistivity. The third column is the true model of Z direction resistivity. The fourth column is the isotropic inversion result

圖14 用于交叉梯度約束的速度模型Fig.14 The velocity model for cross-gradient constraint

圖17展示了Y方向電阻率反演結(jié)果，從圖17b和c可以看到，無約束和約束后的結(jié)果依然有較為顯著的差異.此外，這里給出了圖17b和c的結(jié)果分別與速度模型進(jìn)行每個(gè)單元的交叉梯度項(xiàng)φcgy的計(jì)算，最終結(jié)果如圖18所示.圖18a及b的φcgy分布可以看到，約束反演能較好的使Y方向電阻率結(jié)構(gòu)朝著約束模型靠近.

圖15 (a)為X方向真實(shí)電阻率， (b)為無約束的X方向電阻率反演結(jié)果， (c)為使用結(jié)構(gòu)約束的X方向電阻率反演結(jié)果Fig.15 (a) is the true model of X direction resistivity. (b) is the inversion result of X direction resistivity without constraint. (c) is the inversion result of X direction resistivity with structural constraint

圖16 (a)為圖15b與速度模型的交叉梯度項(xiàng)φcgx分布， (b)為圖15c與速度模型的交叉梯度項(xiàng)φcgx分布Fig.16 (a) is the cross-gradient φcgx distribution of Fig.15b and velocity model. (b) is the cross-gradient φcgx distribution of Fig.15c and velocity model

圖17 (a)為Y方向真實(shí)電阻率， (b)為無約束的Y方向電阻率反演結(jié)果， (c)為使用結(jié)構(gòu)約束的Y方向電阻率反演結(jié)果Fig.17 (a) is the true model of Y direction resistivity. (b) is the inversion result of Y direction resistivity without constraint. (c) is the inversion result of Y direction resistivity with structural constraint

圖19為Z方向電阻率反演結(jié)果，圖19b和c的結(jié)果相較于圖7與圖10中的Z方向電阻率反演結(jié)果而言，有更好的擾動(dòng)結(jié)果.但實(shí)際上圖19中的無約束和約束結(jié)果都只能反映Z方向電阻率淺部的大致形態(tài)，背景值及整體分布都沒有正確的還原，該結(jié)果再一次說明MT對(duì)垂直方向電阻率的靈敏度很弱.我們?cè)趫D20a與b中分別給出了無約束及約束后的Z方向電阻率與速度模型每個(gè)單元的φcgz分布，這里依然可以看到約束后的模型交叉梯度值更小.

在圖21中分別顯示了各向同性、無約束各向異性、約束各向異性反演的步長(zhǎng)曲線及收斂曲線，從步長(zhǎng)曲線可以看到，各向同性反演在最后階段的反演步長(zhǎng)都不為1，這反映出LBFGS在收斂中遇到了擬牛頓方向不穩(wěn)定的情況，說明當(dāng)介質(zhì)具有較強(qiáng)各向異性效應(yīng)時(shí)，常規(guī)的各向同性反演不適用.此外，從rms收斂曲線也可以進(jìn)一步觀察到，各向同性反演很難達(dá)到收斂停止條件.本節(jié)的模型反演實(shí)驗(yàn)再一次展示了各向異性對(duì)反演帶來的影響是非常顯著的，在實(shí)際應(yīng)用中應(yīng)該特別注意，尤其是大構(gòu)造的解譯是否有必要考慮各向異性效應(yīng)，以盡可能減少冗余構(gòu)造的出現(xiàn).

圖18 (a)為圖17b與速度模型的交叉梯度項(xiàng)φcgy分布， (b)為圖17c與速度模型的交叉梯度項(xiàng)φcgy分布Fig.18 (a) is the cross-gradient φcgy distribution of Fig.17b and velocity model. (b) is the cross-gradient φcgy distribution of Fig.17c and velocity model

圖19 (a)為Z方向真實(shí)電阻率, (b)為無約束的Z方向電阻率反演結(jié)果, (c)為使用結(jié)構(gòu)約束的Z方向電阻率反演結(jié)果Fig.19 (a) is the true model of Z direction resistivity. (b) is the inversion result of Z direction resistivity without constraint. (c) is the inversion result of Z direction resistivity with structural constraint

圖20 (a)為圖19b與速度模型的交叉梯度項(xiàng)φcgz分布， (b)為圖19c與速度模型的交叉梯度項(xiàng)φcgz分布Fig.20 (a) is the cross-gradient φcgz distribution of Fig.19b and velocity model. (b) is the cross-gradient φcgz distribution of Fig.19c and velocity model

圖21 (a)為步長(zhǎng)曲線， (b)為rms曲線Fig.21 (a) is the step length curves. (b) is the rms curves

3.4 并行效率分析

由于主軸各向異性正反演仍然具備頻點(diǎn)、XY及YX模式相對(duì)獨(dú)立的計(jì)算任務(wù)，因此很容易使用MPI技術(shù)實(shí)現(xiàn)粗粒度并行.這里使用3.2節(jié)的各向異性正演計(jì)算進(jìn)行并行效率分析，運(yùn)行的系統(tǒng)為Centos Linux系統(tǒng)，內(nèi)存32 G，工作站最大可開辟24進(jìn)程.我們分別使用1，2，4進(jìn)程數(shù)，來展示總體耗時(shí)和加速比.表1統(tǒng)計(jì)結(jié)果表明，三維并行程序得到了接近線性的加速比，體現(xiàn)了較好的并行度.

表1 使用不同進(jìn)程數(shù)的耗時(shí)及加速比Table 1 The consuming time of adopting different processes and speed-up ratio

4 結(jié)論

在復(fù)雜工區(qū)開展多方法地球物理勘探，已經(jīng)是目前常用的流程.本文的目的就是希望能為大地電磁反演輸入一種先驗(yàn)的結(jié)構(gòu)信息，以達(dá)到進(jìn)一步提高大地電磁反演分辨率的目的.基于此，本文開展了基于交叉梯度結(jié)構(gòu)約束的大地電磁主軸各向異性三維反演研究.為了實(shí)現(xiàn)結(jié)構(gòu)約束，在各向異性反演目標(biāo)函數(shù)中加入了交叉梯度約束項(xiàng)，可輸入任意結(jié)構(gòu)實(shí)現(xiàn)本文提出的主軸各向異性結(jié)構(gòu)約束.反演采用了LBFGS法，并進(jìn)一步利用MPI技術(shù)實(shí)現(xiàn)了并行加速，得到了較好的加速比.本文的模型試算可以得到如下結(jié)論：

(1)基于交叉梯度的結(jié)構(gòu)約束方案，可以對(duì)各向同性與主軸各向異性模型起到較好的結(jié)構(gòu)約束效果，在一定程度上改善了大地電磁單方法反演的分辨率；

(2)若地下介質(zhì)存在顯著的主軸各向異性時(shí)，強(qiáng)制使用無約束或約束的各向同性反演，均會(huì)造成地下結(jié)構(gòu)出現(xiàn)嚴(yán)重的假異常.這種情況的出現(xiàn)，可能會(huì)對(duì)實(shí)際構(gòu)造的解釋產(chǎn)生極大的干擾，不利于最終的地質(zhì)建模及解釋.