基于交叉梯度結構約束的大地電磁主軸各向異性并行三維反演

王堃鵬， 王緒本， 曹輝， 藍星， 段長生， 羅威， 原健龍

1 成都理工大學地球物理學院， 成都 610059 2 油氣藏地質及開發工程國家重點實驗室(成都理工大學)， 成都 610059 3 中國地質調查局成都地質調查中心, 成都 610081 4 四川省冶勘設計集團有限公司， 成都 610051 5 贛中南地質礦產勘查研究院， 南昌 330029

0 引言

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

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

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

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

1 基于二次場的MT主軸各向異性三維正演

本文將大地電磁三維正演進行背景場和二次場分解，取時諧因子為e-iω t，且忽略位移電流，可得到如下關于電場的似穩態偏微分方程：

(1)

(2)

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

(3)

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

(4)

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

(5)

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

(6)

(7)

(8)

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

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

2.1 各向異性目標函數

本文采用LBFGS法進行反演迭代計算，由于存在三個方向的參數并且觀測數據為復數阻抗，本文設目標函數為：

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

+λcgyφcgy+λcgzφcgz,

(9)

本文采用仿射線性參數變換(Egbert and Kelbert，2012)對真實模型參數進行如下轉換：

(10)

(11)

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

(12)

2.2 交叉梯度項計算

將(13)式按照三個方向進一步展開有：

(14)

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

(15)

其中，mx與mvx分別為X方向的某個電阻率參數和約束模型參數．對于某個具體的tx_x(i,j,z)計算，這里采用圖1來說明．圖1展示了一個從三維離散體抽取的二維離散面，該面是Y軸與Z軸組成的平面，以中心差分為準則，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中的內部單元循環，可以完成tx_x向量所有元素的計算．由于mvx是輸入的約束模型參數，在反演中固定，實際上最終的tx_x向量是由式(17)的矩陣累加形成，最終拓展到整個模型單元可以表述為：

tx_x=Wx_cg_xmx.

(19)

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

ty_x=Wy_cg_xmx,

tz_x=Wz_cg_xmx.

(20)

(21)

(22)

(23)

2.3 各向異性梯度計算

對于LBFGS來說，梯度計算是最重要的步驟之一，根據目標函數(12)，以及式(21)(22)(23)，對轉換參數求導得到以下梯度表達式：

(24)

式中，數據項的目標函數具體形式如下：

(25)

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

關于式(24)后面的模型項和交叉梯度項的梯度，具體如下：

(26)

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

2.4 LBFGS反演框架

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

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

(3)尋找一個相對合適的更新步長，以最小化φ(mi+αiui)；

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

(27)

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

3 合成算例分析

3.1 各向同性模型無約束與約束反演對比

本節首先使用2.2節的交叉梯度約束項進行各向同性反演試驗，以驗證本文的約束方案首先能應用于各向同性介質的結構約束反演.2.2節給出的是主軸各向異性約束反演，但該方案同樣可應用于各向同性反演，僅僅只需要將其他兩個方向的交叉梯度約束項去掉即可.

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

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

圖2 測點分布及模型俯視圖Fig.2 The measuring sites distribution and top view

圖3 第一列為真實模型，第二列為無交叉梯度約束反演結果，第三列為交叉梯度約束反演結果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

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

3.2 各向同性與主軸各向異性反演對比

3.2.1 各向同性與各向異性無約束反演對比

為了說明在主軸各向異性明顯時，各向同性反演的不足，這里首先使用無約束的各向同性與各向異性大地電磁三維反演程序對圖6模型進行反演，最終的反演結果如圖7所示．圖7第一列是各向同性反演，第二、三、四列為各向異性反演．從圖7的結果可以明顯看出，當地層介質存在明顯主軸各向異性時，圖7第一列顯示出了極強的假構造分布，對于實際情況而言這種結果會嚴重誤導構造解釋．

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

圖6 第一列為X方向電阻率真實模型，第二列為Y方向電阻率真實模型，第三列為Z方向電阻率真實模型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 第一列為各向同性反演結果，第二列為X方向電阻率反演結果，第三列為Y方向電阻率反演結果，第四列為Z方向電阻率反演結果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方向電阻率反演結果，沒有出現第一列強烈的假構造問題.由于X方向電阻率真實值(10 Ωm)相較于Y方向電阻率(30 Ωm)更低，而大地電磁對低阻更靈敏，因此X方向電阻率恢復情況要好于Y方向電阻率．對于第四列的Z方向電阻率，從顏色分布上有一定的擾動，但絕大多數值還是與初始模型差不多，這說明大地電磁法來自高空的平面波場源對地層的垂直電阻率不靈敏(Cao et al., 2018)．

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

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

3.2.2 各向同性與各向異性約束反演對比

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

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

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

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

圖11展示了約束下的各向同性與各向異性反演步長曲線與收斂曲線，從該結果可以看到，二者都體現了非常好的穩定性，但各向同性反演迭代次數更多，且很難達到收斂條件．圖10與圖11的反演結果說明當介質存在明顯主軸各向異性時，對各向同性反演結果產生極大的干擾．

圖10 第一列為各向同性反演結果，第二列為X方向電阻率反演結果，第三列為Y方向電阻率反演結果，第四列為Z方向電阻率反演結果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)為步長曲線，(b)為rms曲線Fig.11 (a) is the step length curves. (b) is the rms curves

3.3 復雜各向異性模型反演

在3.2節，我們證明了在主軸各向異性效應顯著時，強制使用各向同性的無約束和約束反演都會帶來嚴重的假構造分布．由于模型較為簡單，各向異性約束反演并沒有特別好的體現出交叉梯度的優勢和特點，本節我們將建立一個形態較為復雜的模型，探討交叉梯度項對大地電磁主軸各向異性反演的約束效果．

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

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

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

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

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

圖13 第一列為X方向真實電阻率，第二列為Y方向真實電阻率，第三列為Z方向真實電阻率，第四列為各向同性反演結果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方向電阻率反演結果，從圖17b和c可以看到，無約束和約束后的結果依然有較為顯著的差異.此外，這里給出了圖17b和c的結果分別與速度模型進行每個單元的交叉梯度項φcgy的計算，最終結果如圖18所示.圖18a及b的φcgy分布可以看到，約束反演能較好的使Y方向電阻率結構朝著約束模型靠近.

圖15 (a)為X方向真實電阻率， (b)為無約束的X方向電阻率反演結果， (c)為使用結構約束的X方向電阻率反演結果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與速度模型的交叉梯度項φcgx分布， (b)為圖15c與速度模型的交叉梯度項φ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方向真實電阻率， (b)為無約束的Y方向電阻率反演結果， (c)為使用結構約束的Y方向電阻率反演結果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方向電阻率反演結果，圖19b和c的結果相較于圖7與圖10中的Z方向電阻率反演結果而言，有更好的擾動結果.但實際上圖19中的無約束和約束結果都只能反映Z方向電阻率淺部的大致形態，背景值及整體分布都沒有正確的還原，該結果再一次說明MT對垂直方向電阻率的靈敏度很弱.我們在圖20a與b中分別給出了無約束及約束后的Z方向電阻率與速度模型每個單元的φcgz分布，這里依然可以看到約束后的模型交叉梯度值更小.

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

圖18 (a)為圖17b與速度模型的交叉梯度項φcgy分布， (b)為圖17c與速度模型的交叉梯度項φ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方向真實電阻率, (b)為無約束的Z方向電阻率反演結果, (c)為使用結構約束的Z方向電阻率反演結果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與速度模型的交叉梯度項φcgz分布， (b)為圖19c與速度模型的交叉梯度項φ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)為步長曲線， (b)為rms曲線Fig.21 (a) is the step length curves. (b) is the rms curves

3.4 并行效率分析

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

表1 使用不同進程數的耗時及加速比Table 1 The consuming time of adopting different processes and speed-up ratio

4 結論

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

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

(2)若地下介質存在顯著的主軸各向異性時，強制使用無約束或約束的各向同性反演，均會造成地下結構出現嚴重的假異常.這種情況的出現，可能會對實際構造的解釋產生極大的干擾，不利于最終的地質建模及解釋.