考慮關(guān)斷時(shí)間的地面瞬變電磁三維帶地形反演

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

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院， 西安 7100542 長(zhǎng)安大學(xué)地球物理場(chǎng)多參數(shù)綜合模擬實(shí)驗(yàn)室， 西安 7100543 中國(guó)地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所， 河北 廊坊 065000

0 引言

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

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

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

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

1 非結(jié)構(gòu)時(shí)間域有限元正演理論

從時(shí)間域麥克斯韋方程組出發(fā)，并忽略位移電流項(xiàng)，建立電場(chǎng)擴(kuò)散方程

其中r是空間位置矢量，t是時(shí)間，e(r,t)是r處t時(shí)刻的電場(chǎng)，js(r,t)是源電流密度，μ和σ分別是磁導(dǎo)率和電導(dǎo)率.采用基于非結(jié)構(gòu)四面體網(wǎng)格的矢量有限元方法對(duì)電場(chǎng)擴(kuò)散方程進(jìn)行空間離散，則任意四面體單元內(nèi)的電場(chǎng)均可表示為(齊彥福等，2017)

(2)

(3)

其中M和S分別為質(zhì)量和剛度矩陣，Js是電流源項(xiàng).

采用二階后退歐拉離散格式對(duì)(3)式進(jìn)行時(shí)間離散，得到無(wú)條件穩(wěn)定的隱式差分方程：

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

(4)

其中em(t)表示第m時(shí)刻的電場(chǎng)，Δt為時(shí)間步長(zhǎng).

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

關(guān)于邊界條件，本文采用均勻半空間狄利克雷邊界條件(Yin et al.，2016)，假設(shè)計(jì)算區(qū)域外邊界的切向電場(chǎng)分量為0，即

(5)

關(guān)于初始條件，由于在供電前，沒(méi)有發(fā)射信號(hào)，因此空間任意處的電場(chǎng)均為0，即

e(r,0)=0.

(6)

然后，采用多波前并行求解器MUMPS(Amestoy et al.，2006)對(duì)(4)式進(jìn)行求解即可獲得各個(gè)時(shí)刻空間電場(chǎng)分布，再利用法拉第電磁感應(yīng)定律計(jì)算接收點(diǎn)處的觀測(cè)場(chǎng)值.

2 瞬變電磁三維反演算法

2.1 目標(biāo)函數(shù)建立及正則化方法

首先建立目標(biāo)函數(shù)

(7)

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

(8)

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

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

2.2 L-BFGS反演算法

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

(9)

mn+1=mn+αnqn.

(10)

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

(11)

其中0<β<0.5，β<γ<1，F(xiàn)為正演算子.

2.3 目標(biāo)函數(shù)梯度的計(jì)算

根據(jù)(9)式可知，在L-BFGS算法中近似計(jì)算海森矩陣逆的基礎(chǔ)是獲得目標(biāo)函數(shù)的梯度，本文基于伴隨正演算法計(jì)算目標(biāo)函數(shù)的梯度.首先通過(guò)對(duì)目標(biāo)函數(shù)求導(dǎo)可以獲得目標(biāo)函數(shù)梯度的表達(dá)式

(12)

(13)

由于(13)式中等號(hào)右端第二項(xiàng)很容易計(jì)算，因此計(jì)算目標(biāo)函數(shù)梯度的核心是計(jì)算JTr.我們首先將所有時(shí)刻的正演問(wèn)題簡(jiǎn)寫為

Ae=b，

(14)

則有

e=A-1b.

(15)

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

(16)

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

(17)

則有

JTr=GTA-TQTr，

(18)

JTr=GTv，

(19)

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

ATv=R，

(20)

進(jìn)行計(jì)算.通過(guò)逆向時(shí)間迭代獲得向量v之后即可代入(19)式計(jì)算JTr，進(jìn)而獲得目標(biāo)函數(shù)的梯度，完成模型更新.

3 理論模型反演

3.1 地形對(duì)三維反演結(jié)果的影響特征

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

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

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

圖3 考慮和忽略地形影響的反演結(jié)果對(duì)比圖(a)—(c) 起伏地表模型反演結(jié)果； (d)—(f) 水平地表模型反演結(jié)果.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 數(shù)據(jù)擬合結(jié)果(a)—(c) 含噪聲的理論數(shù)據(jù)； (d)—(f) 起伏地表反演模型的正演數(shù)據(jù)； (g)—(i) 水平地表反演模型的正演數(shù)據(jù).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 考慮和忽略地形影響的反演迭代參數(shù)(a) 均方根誤差； (b) 目標(biāo)函數(shù)； (c) 正則化因子λ.Fig.5 Inversion parameters versus iterations for considering and ignoring topography(a) Root-mean-square error; (b) Objective function; (c) Parameter λ.

3.2 關(guān)斷時(shí)間對(duì)三維反演結(jié)果的影響特征

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

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

圖7 發(fā)射電流波形Fig.7 Transmitting current

圖8 考慮和忽略關(guān)斷時(shí)間影響的反演結(jié)果對(duì)比圖(a)—(c) 考慮關(guān)斷時(shí)間的三維反演結(jié)果； (d)—(f) 理論階躍波形的反演結(jié)果.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測(cè)線數(shù)據(jù)擬合情況Fig.9 Data fitting for survey line x=50 m

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

4 實(shí)測(cè)數(shù)據(jù)反演

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

圖11 發(fā)射源和測(cè)線空間分布Fig.11 Transmitting loop and survey lines over topographic survey area

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

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

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

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

5 結(jié)論

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

致謝特別向各位審稿人和編輯同志對(duì)本文提出的建設(shè)性意見(jiàn)表示感謝.