考慮關斷時間的地面瞬變電磁三維帶地形反演

齊彥福， 智慶全,3， 李貅*， 景旭, 戚志鵬，孫乃泉， 周建美， 劉文韜

1 長安大學地質工程與測繪學院， 西安 7100542 長安大學地球物理場多參數綜合模擬實驗室， 西安 7100543 中國地質科學院地球物理地球化學勘查研究所， 河北 廊坊 065000

0 引言

隨著經濟的快速發展，我國對礦產資源的需求不斷增加，然而經過數十年的不間斷開采，地形相對平坦、地質條件簡單的淺部礦藏已經開采殆盡，找礦工作正在向地形起伏嚴重的山區進軍.地面瞬變電磁法作為一種低成本、高效率的地球物理勘探方法，已經在礦產資源勘查領域獲得廣泛應用(薛國強等，2007；底青云等，2019).然而，在過去的數十年中受到計算機條件和數據解釋技術的制約，地面瞬變電磁數據解釋工作主要以視電阻率成像和一維反演方法為主(薛國強等，2008；石顯新等,2009；武軍杰等，2013)，此類方法在地形平坦、地下介質成層性明顯的地區可以獲得很好的應用效果，然而在地形起伏嚴重、三維地電結構發育的山區，利用一維處理解釋技術很難獲取準確的地下三維電性分布信息.開發適用于起伏地形的瞬變電磁三維反演算法，已成為瞬變電磁數據精細化解釋的關鍵.

作為反演的基礎，前人在瞬變電磁三維正演算法方面開展了大量工作，開發了一系列三維正演方法(湯井田等，2007)，例如頻時轉換法(殷長春等，2013)、時間域有限差分法(Wang and Hohmann，1993；Commer and Newman，2004)、時間域有限體積法(Haber et al.，2007；Oldenburg et al.，2013)和時間域有限元法(Yin et al.，2016；Cai et al.，2017).頻時轉換法首先計算三維頻率域電磁響應，然后采用頻時轉換技術將其轉換到時間域，完成瞬變電磁響應模擬.時間域有限差分法是一種基于交錯網格和顯式時間差分格式的直接時間域正演算法(孫懷鳳等，2013)，該方法數學概念直觀，表達簡單，但是需要滿足Courant穩定性條件，使得網格尺寸和時間步長均受到嚴格限制.近年來，隨著計算機性能的大幅提升和大型方程組求解技術的快速發展，基于隱式時間差分格式的三維正演算法取得了重大技術突破，其優勢逐漸展現出來.時間域有限體積法(Yang et al.，2014)和時間域有限元法(齊彥福等，2017)通過采用后退歐拉時間離散格式對控制方程進行無條件穩定的隱式時間離散，極大地放寬了對時間步長的限制.與有限體積法所采用的規則六面體網格不同，時間域有限元法可以采用非規則四面體網格進行空間離散.基于四面體網格的靈活性，該方法在處理起伏地形和復雜地電結構時顯示出明顯優勢.Li等(2018)采用時間域有限元方法模擬了復雜形狀發射源地面瞬變電磁響應，結果表明地形起伏和發射線圈形狀的變化對觀測響應影響極大，其導致觀測響應發生復雜畸變.Zeng等(2019)進一步分析了關斷時間對瞬變電磁響應的嚴重影響.

在三維反演方面，目前三維電磁反演算法主要有高斯牛頓法、非線性共軛梯度法和擬牛頓法.高斯牛頓法通過近似計算Hessian矩陣和目標函數的梯度獲得模型改變量，然后采用線性搜索方式獲得最佳的模型更新步長，實現模型迭代更新.該算法具有超線性收斂特性(Haber et al.，2007)，但是需要計算靈敏度矩陣的乘積來近似Hessian矩陣，導致每次迭代過程的計算量巨大.非線性共軛梯度法首先計算目標函數在當前參考模型下的負梯度，并在梯度的共軛方向搜索最優步長使得目標函數達到極小值，進而實現反演模型更新(翁愛華等，2012).由于該方法在每次反演迭代過程中無需計算靈敏度矩陣，因此計算量明顯減少，然而該算法為線性收斂，反演收斂速度較慢(Newman and Commer，2005).L-BFGS法是一種有限內存的擬牛頓反演算法，該算法通過伴隨正演計算目標函數的梯度，再采用迭代方式近似計算Hessian矩陣的逆，因此每次反演迭代只需要一次正演和一次伴隨正演即可完成模型更新.與非線性共軛梯度法相比，L-BFGS算法在計算時間和模型反演分辨率上更具優勢(劉云鶴和殷長春，2013).Liu等(2019)基于非結構時間域有限元法和L-BFGS算法實現了瞬變電磁三維帶地形反演.

三維反演技術的發展在提高瞬變電磁數據解釋精度的同時，也為改進觀測模式、提高工作效率提供了技術支撐.在傳統的定源回線裝置野外工作中，由于受到中心回線近似解釋方法的限制，數據采集只能在發射線圈中心三分之一區域進行，但是三維反演技術的出現徹底打破了這種限制，可以在發射線圈內外任意位置進行三維觀測，使得瞬變電磁的觀測方式變得更加靈活，同時這種三維觀測方式極大地提高了野外施工效率.然而，觀測模式的改進也帶了一些新的問題.考慮到信號強度隨收發距離的幾何衰減規律，要實現基于三維觀測模式的大深度探測必須加大發射功率.受到當前電子元器件條件的限制，當發射大電流時，無法實現瞬間斷電，需要較長關斷時間，此時不能夠再采用理論階躍波形近似.由于關斷時間對瞬變電磁響應具有嚴重影響(孫懷鳳等，2013；楊海燕等，2019；Zeng等，2019)，因此在瞬變電磁三維反演過程中必須考慮關斷時間.為此本文開展考慮關斷時間的地面瞬變電磁三維帶地形反演算法研究，以提高瞬變電磁數據解釋的精細化程度.我們首先介紹本文所采用的非結構時間域有限元正演理論和瞬變電磁三維反演算法，然后將本文算法應用于理論和實測數據三維反演中以檢驗本文算法的有效性，并分析地形和關斷時間對三維反演效果的影響特征.

1 非結構時間域有限元正演理論

從時間域麥克斯韋方程組出發，并忽略位移電流項，建立電場擴散方程

其中r是空間位置矢量，t是時間，e(r,t)是r處t時刻的電場，js(r,t)是源電流密度，μ和σ分別是磁導率和電導率.采用基于非結構四面體網格的矢量有限元方法對電場擴散方程進行空間離散，則任意四面體單元內的電場均可表示為(齊彥福等，2017)

(2)

(3)

其中M和S分別為質量和剛度矩陣，Js是電流源項.

采用二階后退歐拉離散格式對(3)式進行時間離散，得到無條件穩定的隱式差分方程：

(3M+2ΔtS)em+2(t)=M(4em+1(t)-em(t))

(4)

其中em(t)表示第m時刻的電場，Δt為時間步長.

在瞬變電磁法工作過程中，需要在地表敷設一個很大的發射線圈，當遇到地形起伏時，發射線圈也會隨著地形起伏，使得發射線圈的形狀和有效面積發生變化，為了模擬發射線圈的實際形狀和尺寸，本文采用基于偶極子離散的場源處理方法(Yin et al.，2016)精細模擬實際場源.通過將發射線圈離散成若干段收尾相連的短導線，每段導線當作一個電偶極子進行處理，并考慮發射線圈與四面體單元任意空間相交關系，實現復雜發射場源的模擬.為了模擬斷電時間的影響，本文采用瞬時電流脈沖技術(齊彥福等，2017)，通過在每個時刻的右端源項中加載瞬時電流脈沖強度，實現考慮關斷時間的瞬變電磁響應三維正演.

關于邊界條件，本文采用均勻半空間狄利克雷邊界條件(Yin et al.，2016)，假設計算區域外邊界的切向電場分量為0，即

(5)

關于初始條件，由于在供電前，沒有發射信號，因此空間任意處的電場均為0，即

e(r,0)=0.

(6)

然后，采用多波前并行求解器MUMPS(Amestoy et al.，2006)對(4)式進行求解即可獲得各個時刻空間電場分布，再利用法拉第電磁感應定律計算接收點處的觀測場值.

2 瞬變電磁三維反演算法

2.1 目標函數建立及正則化方法

首先建立目標函數

(7)

其中m是M維地電模型參數向量，dobs是N維觀測數據向量，本文中dobs所包含的元素是各個測點不同時間道的dBz/dt響應，Wd和Wm分別是數據和模型方差矩陣，d是當前模型的正演數據，mref是先驗參考模型，λ是正則化因子.(7)式中第二項是正則化項，其主要作用是壓縮反演模型空間，改善反演的穩定性.利用模型方差矩陣Wm可以約束當前單元與相鄰單元間模型參數的空間變化率，從而獲得光滑的反演模型.本文根據Key(2016)提出的二維非結構網格空間約束矩陣建立方法，綜合考慮與當前單元相鄰的所有單元的電性相關性，并將其推廣到三維情況(Liu et al.，2019)，建立如下Wm矩陣

(8)

圖1 模型空間約束Fig.1 Spatial constraint for inversion model

λ的選取至關重要，將直接影響反演的收斂速度和反演效果，本文采用“金屬冷卻法”(Haber，2014)選取每次反演迭代的λ，即首先人為設定一個較大的初值加強模型約束強度來保證反演的穩定性，然后利用當前λ進行反演迭代，當數據擬合差下降非常緩慢或不下降時減小λ，降低模型約束強度，通過不斷重復上述過程實現數據擬合.

2.2 L-BFGS反演算法

本文采用Nocedal(1980)提出的有限內存擬牛頓法(L-BFGS)實現瞬變電磁數據三維反演.擬牛頓算法(BFGS)的核心是近似計算海森矩陣的逆H，為此首先給定一個對稱正定矩陣H0作為初始海森矩陣的逆，本文中H0取為單位對角矩陣，然后利用

(9)

mn+1=mn+αnqn.

(10)

本文采用線性搜索的方式獲得αn，首先給定初始αn=1，然后逐步減小αn，使其滿足Wolf條件(Nocedal and Wright， 2006)

(11)

其中0<β<0.5，β<γ<1，F為正演算子.

2.3 目標函數梯度的計算

根據(9)式可知，在L-BFGS算法中近似計算海森矩陣逆的基礎是獲得目標函數的梯度，本文基于伴隨正演算法計算目標函數的梯度.首先通過對目標函數求導可以獲得目標函數梯度的表達式

(12)

(13)

由于(13)式中等號右端第二項很容易計算，因此計算目標函數梯度的核心是計算JTr.我們首先將所有時刻的正演問題簡寫為

Ae=b，

(14)

則有

e=A-1b.

(15)

由于瞬變電磁法觀測的是dB/dt響應，由(15)式可以得到

(16)

其中Q是插值矩陣.靈敏度矩陣可寫成

(17)

則有

JTr=GTA-TQTr，

(18)

JTr=GTv，

(19)

其中v=A-TQTr.v通過求解伴隨正演方程

ATv=R，

(20)

進行計算.通過逆向時間迭代獲得向量v之后即可代入(19)式計算JTr，進而獲得目標函數的梯度，完成模型更新.

3 理論模型反演

3.1 地形對三維反演結果的影響特征

為了檢驗本文算法的有效性，我們首先設計了如圖2所示三維起伏地表模型，地形最大起伏高差約為150 m，背景圍巖電阻率為100 Ωm，在觀測區域的中心位置埋藏一個電阻率為10 Ωm的良導梯形體，梯形體頂面邊長為200 m，底面邊長為400 m，梯形體的高為100 m，頂面埋深約為200 m.觀測系統采用定源回線裝置，發射線框邊長為650 m，發射電流為10 A的階躍波，基于傳統觀測模式在發射線框內部進行測量，為了完成整個區域的觀測，共布置了9個發射線圈，每個線框內均勻分布8×8=64個測點，共576個測點，點距和線距均為50 m，每個測點觀測從10 μs到10 ms的31道對數等間隔dBz/dt響應數據.我們在理論數據中加入5%的隨機噪聲作為觀測數據，然后利用本文開發的三維反演算法進行反演.本文中的所有反演算例均在PC機上進行，CPU主頻3.20 GHz，共6個計算核心，內存64 GB.

為了研究地形對三維反演結果的影響，本文分別采用起伏地形和水平地表兩種模型完成反演，初始模型和參考模型均采用100 Ωm的半空間模型，初始正則化因子λ均設定為10，然后采用“金屬冷卻法”獲得每次反演迭代的λ.反演網格的核心區域范圍如圖3所示，然后在x、y、z三個方向分別設定100000 m的擴邊區域，使其滿足狄利克雷邊界條件.起伏地形反演網格包含478958個四面體單元，而水平地表反演網格包含435644個四面體單元.圖3展示了經過13次迭代的反演結果.從圖中可以看出帶地形反演結果(圖3a—c)與真實模型(圖2)吻合較好，能夠清晰反映出地下的電性分布情況，較為準確地顯示了良導目標體的位置和尺寸.然而水平地表模型的反演結果(圖3d—f)雖然也可以獲得地下目標的大體分布情況，然而為了強制擬合地形數據，在反演結果中出現了大量假異常，而且這些虛假異常不僅局限于地表，在深部同樣存在，虛假異常的出現將給后續的地質解釋造成很多困難.該反演結果有效驗證了本文算法的可靠性，同時也說明了起伏地形對三維反演結果會產生巨大影響.圖4和圖5展示了兩組反演結果的數據擬合情況.可以看出，考慮地形的三維瞬變電磁反演可以很好地擬合觀測數據，而且反演收斂速度很快，反演共耗時37.39 h，所占內存為22.5 GB.然而，在不考慮地形情況下，通過產生一些局部的虛假異常同樣可以達到數據擬合的目的，但是在數據擬合到一定程度后無法找到合理的三維模型實現進一步的數據擬合，而且由于不合理反演模型的試探次數較多，反演共耗時128.36 h，所占內存為20.8 GB.

圖2 起伏地表梯形良導體模型(a)和(b)分別是不同角度的三維透視圖.Fig.2 Topographic model with a trapezoidal conductor(a) and (b) are perspective views from different directions.

圖3 考慮和忽略地形影響的反演結果對比圖(a)—(c) 起伏地表模型反演結果； (d)—(f) 水平地表模型反演結果.Fig.3 Comparison of inversion results with and without topography(a)—(c) are inversion results for a topographic earth; (d)—(f) are for a flat earth.

圖4 數據擬合結果(a)—(c) 含噪聲的理論數據； (d)—(f) 起伏地表反演模型的正演數據； (g)—(i) 水平地表反演模型的正演數據.Fig.4 Data fitting for invesion results(a)—(c) Synthetic data; (d)—(f) Predicted data from the inverted model with topography; (g)—(i) Predicted data from the inverted model without topography.

圖5 考慮和忽略地形影響的反演迭代參數(a) 均方根誤差； (b) 目標函數； (c) 正則化因子λ.Fig.5 Inversion parameters versus iterations for considering and ignoring topography(a) Root-mean-square error; (b) Objective function; (c) Parameter λ.

3.2 關斷時間對三維反演結果的影響特征

本文設計了如圖6所示的起伏地表模型，地形最大起伏高差約為380 m.在100 Ωm的半空間中埋藏一個電阻率為10 Ωm的良導L形目標體，目標體位于發射線圈正下方，沿x和y方向的長度均為700 m，寬度和高度均為200 m，頂面埋深約為300 m.觀測裝置采用定源回線裝置，發射線圈邊長為1000 m，發射如圖7所示的電流波形，其中供電時間ton=19 ms，關斷時間tramp=100 μs，峰值電流為20 A.在發射線圈內外同時進行觀測，共觀測16×16=256個測點，點距和線距均為100 m.每個測點觀測從發射電流完全關斷后10 μs到10 ms的31道對數等間隔dBz/dt響應數據.我們在理論數據中加入5%的隨機噪聲作為觀測數據，然后利用本文開發的三維反演算法進行反演.

為了分析關斷時間對三維反演結果的影響，我們采用考慮關斷時間和理論階躍波兩種電流波形對觀測數據進行反演.初始模型和參考模型均采用100 Ωm的半空間模型，初始正則化因子λ均為0.01，然后采用“金屬冷卻法”獲得每次反演迭代的λ.反演網格的核心區域范圍如圖8所示，然后在x、y、z三個方向分別設定100000 m的擴邊區域，使其滿足狄利克雷邊界條件.反演模型網格包含535752個四面體單元.圖8展示了經過30次迭代后的兩種電流波形的反演結果.從圖中可以看出考慮關斷時間的反演結果(圖8a—c)與真實模型(圖6)吻合較好，能夠清晰反映出地下的電性分布情況.然而理論階躍波形的反演結果(圖8d—f)雖然也可以獲得地下目標的大體分布，但是在反演結果中出現了很多假異常，且大部分集中在地表，這是由于關斷時間的變化主要影響瞬變電磁響應斷電后早期的電磁響應，這種影響隨著時間的增加逐漸減弱.虛假異常的出現嚴重影響了瞬變電磁數據三維解釋精度.圖9和圖10展示了兩組反演結果的數據擬合情況.可以看出，考慮關斷時間的三維瞬變電磁反演可以很好地擬合觀測數據，而且反演收斂速度很快，反演共耗時12.56 h.然而，在不考慮關斷時間的情況下數據擬合情況并不理想，尤其是早期數據無法擬合(如圖9所示)，導致錯誤的模型搜索方向，進而在一定程度上影響深部目標體的反演效果.由于不合理反演模型的試探次數較多，反演共耗時16.37 h.兩種模型反演所占內存均為20.3 GB.

圖7 發射電流波形Fig.7 Transmitting current

圖8 考慮和忽略關斷時間影響的反演結果對比圖(a)—(c) 考慮關斷時間的三維反演結果； (d)—(f) 理論階躍波形的反演結果.Fig.8 Comparison of inversion results for considering and ignoring ramp time(a)—(c) Inversion results considering ramp time; (d)—(f) Inversion results ignoring ramp time.

圖9 x=50 m測線數據擬合情況Fig.9 Data fitting for survey line x=50 m

圖10 考慮和忽略關斷時間影響的反演迭代參數(a) 均方根誤差； (b) 目標函數； (c)正則化因子λ.Fig.10 Inversion parameters versus iterations for considering and ignoring ramp time(a) Root-mean-square error; (b) Objective function; (c) Parameter λ.

4 實測數據反演

我們將本文開發的三維反演算法應用于青海夏日哈木礦區的實測數據解釋中.夏日哈木銅鎳礦區地處柴達木盆地西南緣，是東昆侖成礦帶新發現的一處超大型巖漿熔離型銅鎳硫化物礦床，超基性巖體具有多期次侵入的特征，巖性較為復雜，整體上深部礦體品位相對上部較高.礦區出露地層為古元古代金水口巖群白沙河組，主要巖性有黑云斜長片麻巖、云母二長片麻巖、斜長角閃巖、大理巖等.區域內巖漿巖較為發育，主要形成于古特提斯造山旋回的不同階段，主要由早二疊世閃長巖、晚三疊世花崗巖等組成(杜瑋等，2014).圍巖導電性較差，最高電阻率可達1000 Ωm以上.然而工業礦體的導電性良好，最低電阻率可達3.5 Ωm.研究區內地貌主要為戈壁荒漠及殘山景觀，近山脊基巖裸露、地形較陡峻，溝谷風塵黃土覆蓋.為完成對地下礦體的勘探，在測試區內敷設一個600 m×600 m的發射線框，共布置5條測線，線距80 m，點距50 m，共83個測點(如圖11所示)，其中L2和L3測線與地質勘探線重合.發射如圖7所示的電流波形，其中供電時間ton=20 ms，關斷時間tramp=1 ms，峰值電流為15 A.每個測點觀測從發射電流完全關斷后54 μs到16 ms的31道對數等間隔dBz/dt響應數據.

圖11 發射源和測線空間分布Fig.11 Transmitting loop and survey lines over topographic survey area

為了更加準確地反映實際情況，我們在反演過程中同時考慮了起伏地形和關斷時間的影響.初始模型和參考模型均采用200 Ωm的半空間模型，初始正則化因子λ設定為0.01，然后采用“金屬冷卻法”獲得每次反演迭代的λ.反演網格的核心區域范圍如圖12所示，然后在x、y、z三個方向分別設定100000 m的擴邊區域，使其滿足狄利克雷邊界條件.反演模型網格包含373072個四面體單元.經過38次迭代反演收斂.圖12展示了考慮關斷時間的三維帶地形反演結果透視圖，從圖中可以看出，良導礦體沿x軸方向延伸，大體垂直于測線方向.圖13顯示三維反演展示的低阻區域與礦體的空間位置基本吻合，礦體埋深約為300 m，礦體主要分布在y=7853到y=8253之間，這一結果進一步驗證了本文算法的有效性.從圖14可以看出在反演初始階段數據擬合差和目標函數下降較快，隨著迭代的進行，下降速度逐漸減小，此時需要減小正則化因子λ，減弱對模型的約束強度，實現進一步的數據擬合.反演共耗時20.11 h，占用內存約為17.9 GB.

圖12 三維反演結果透視圖Fig.12 Perspective view of the 3D inversion model

圖13 三維反演結果與地質勘探線剖面圖Fig.13 3D inversion model and geological profile

圖14 反演迭代參數(a) 均方根誤差； (b) 目標函數； (c) 正則化因子λ.Fig.14 Inversion parameters versus iterations(a) Root-mean-square error; (b) Objective function; (c) Parameter λ.

5 結論

本文通過結合非結構時間域有限元法和L-BFGS反演算法，成功實現了考慮關斷時間的地面瞬變電磁三維帶地形反演.利用非結構網格的靈活性精細模擬起伏地形的變化，采用偶極子離散處理技術模擬發射源的實際形狀，基于瞬時電流脈沖技術模擬關斷時間的影響，最后應用L-BFGS反演算法進行模型迭代更新，完成了基于三維觀測模式的地面瞬變電磁數據三維反演.理論和實測數據的三維反演結果表明地形和關斷時間對反演模型影響較大，其中地形起伏導致發射線圈形狀以及測點和源的空間位置關系發生變化，尤其是起伏地形和地下三維異常體相互耦合導致觀測響應產生了十分復雜的畸變，在反演過程中如果利用水平地表模型強制擬合地形數據將導致反演模型中地表出現大量假異常，而且假異常不只是出現在地表，這給地質解釋造成一定困難.而關斷時間主要影響早期觀測數據，如果利用理論階躍波形響應擬合具有一定關斷時長的觀測數據將導致地表出現大量虛假異常，當關斷時間較長時甚至無法找到合理的三維模型實現數據擬合.早期數據無法擬合將導致錯誤的模型搜索方向，進而影響深部目標體的反演效果.開發基于三維觀測模式的三分量反演算法是我們未來的工作重點.

致謝特別向各位審稿人和編輯同志對本文提出的建設性意見表示感謝.