齊彥福，李 貅，孫乃泉，戚志鵬，周建美，張明晶，邢 濤

1.長安大學地質(zhì)工程與測繪學院，西安 710054 2.長安大學地球物理場多參數(shù)綜合模擬實驗室，西安 710054 3.山東省交通規(guī)劃設(shè)計院，濟南 250031 4.北京探創(chuàng)資源科技有限公司，北京 100071

0 引言

瞬變電磁法是一種建立在電磁感應原理基礎(chǔ)上的時間域電磁探測方法[1]。該方法首先發(fā)射大功率一次電磁信號，然后在斷電間歇觀測導電大地中感應渦流產(chǎn)生的二次電磁場，通過分析二次場的空間和時間分布特征獲得地下電性結(jié)構(gòu)信息[2-3]。根據(jù)發(fā)射源裝置形式的不同，瞬變電磁系統(tǒng)可以分為磁性源瞬變電磁系統(tǒng)和電性源瞬變電磁系統(tǒng)。磁性源瞬變電磁系統(tǒng)采用不接地回線圈作為發(fā)射源，然后利用接收線圈觀測感應電壓，因此磁性源瞬變電磁系統(tǒng)也被稱為回線源瞬變電磁系統(tǒng)。該方法不受接地條件的限制，施工方便，對高阻覆蓋層具有較強的穿透能力[4]。電性源瞬變電磁系統(tǒng)采用接地長導線作為發(fā)射源，然后采用接地電極和接收線圈觀測電場和感應電壓，根據(jù)收發(fā)距的不同可分為長偏移距瞬變電磁[5]和短偏移距瞬變電磁[6]兩種工作形式。其中：長偏移距瞬變電磁系統(tǒng)利用數(shù)千米長的接地導線向地下發(fā)射電流，在3～6倍探測深度的偏移距范圍觀測電磁響應，其探測深度可達幾千米至十幾千米[5]；而短偏移距瞬變電磁系統(tǒng)在小于2倍探測深度的偏移距范圍內(nèi)觀測純二次場[4]。相比于前者，短偏移距瞬變電磁系統(tǒng)的收發(fā)距更小，收發(fā)布置具有更強的靈活性，可根據(jù)測區(qū)的地形條件靈活布置；由于收發(fā)距較小，觀測信號具有更高的信噪比，同時對發(fā)射機的功率和性能要求極大降低，使得野外施工更加方便[4]。

為此，前人針對地形影響開展了大量工作[14-16]。Li等[13]模擬了起伏地表條件下磁性源地面瞬變電磁響應，結(jié)果表明地形對觀測響應影響極大，其導致觀測響應發(fā)生復雜畸變。薛國強等[12]通過分析地形影響的基本規(guī)律，針對磁性源中心回線裝置開展了地形影響快速校正方法研究，并成功應用于理論和實測數(shù)據(jù)中。Tang等[17]利用積分方程法模擬了層狀大地起伏地表模型水平電偶極子源長偏移距瞬變電磁剖面響應，并分析了地形影響特征。然而，與長偏移距瞬變電磁系統(tǒng)相比，短偏移距瞬變電磁系統(tǒng)由于收發(fā)距小，場源影響和地形相互耦合將使得觀測響應更加復雜。此外，當前電磁探測技術(shù)正在向三維觀測模式轉(zhuǎn)變，全面了解三維地形影響的空間分布規(guī)律成為提高瞬變電磁數(shù)據(jù)解釋精度的關(guān)鍵，因此開展系統(tǒng)的短偏移距瞬變電磁地形影響特征研究十分必要。

作為分析地形影響的重要工具，三維數(shù)值模擬技術(shù)近年來發(fā)展迅速。目前，三維瞬變電磁正演方法主要分為三大類：1)頻時轉(zhuǎn)換算法。該方法首先采用有限差分、有限元、有限體積或積分方程等數(shù)值算法模擬較寬頻帶內(nèi)的頻率域電磁響應，然后利用頻時轉(zhuǎn)換技術(shù)[18]將頻率域響應轉(zhuǎn)換到時間域。2)基于顯式差分格式的時間域直接求解方法。時域有限差分算法是此類方法的代表，該算法采用向前差分的時間離散格式，結(jié)合交錯網(wǎng)格技術(shù)建立顯式時間迭代方程，通過交替求解電磁場的方式實現(xiàn)瞬變電磁正演[19-20]。然而該方法需要滿足Courant穩(wěn)定性條件，對網(wǎng)格尺寸和時間步長的限制比較嚴格。3)基于隱式差分格式的時間域直接求解方法。該方法采用后退歐拉離散格式，可以獲得無條件穩(wěn)定的隱式差分方程，極大地放寬了對時間步長和網(wǎng)格尺寸的限制，具有更高的穩(wěn)定性[21]。根據(jù)空間離散方式的不同可分為時間域有限體積法[22]和時間域有限元法[23-24]。與前者相比，時間域有限元法可以采用非規(guī)則四面體網(wǎng)格進行模型剖分，因此其在處理起伏地形和復雜地電結(jié)構(gòu)時更具優(yōu)勢。

本文基于非結(jié)構(gòu)時間域有限元算法開展電性源短偏移距瞬變電磁地形影響研究。利用非結(jié)構(gòu)四面體網(wǎng)格的靈活性精細刻畫起伏地形，基于后退歐拉離散格式的穩(wěn)定性合理設(shè)定時間步長，提高正演效率。首先簡要介紹非結(jié)構(gòu)時間域有限元正演理論，并對我們所編寫的三維正演程序進行算法精度驗證，然后著重分析地形影響特征，通過模擬大量典型起伏地表模型的三維電磁響應分析起伏地形對短偏移距瞬變電磁觀測響應的影響規(guī)律。

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

我們從時間域麥克斯韋方程組出發(fā)，通過消去磁場并忽略位移電流項建立電場擴散方程：

(1)

式中：μ為磁導率；e(r,t)為電場；r為空間位置(x,y,z)矢量；t為時間；σ為電導率；js(r,t)為源電流密度。首先采用非結(jié)構(gòu)四面體網(wǎng)格對整個計算區(qū)域進行剖分，利用四面體網(wǎng)格的靈活性精細刻畫起伏地形和復雜地電結(jié)構(gòu)。然后利用矢量插值基函數(shù)[25]對電場進行空間離散，則任意四面體單元內(nèi)的電場可表示為

(2)

(3)

式中：A為質(zhì)量矩陣；C為剛度矩陣；S為電流源項。對于第k個單元：

式中：V(k)為第k個單元的體積；dV為體積分。接下來，我們采用二階后退歐拉離散格式對式(3)進行時間離散，得到無條件穩(wěn)定的隱式差分方程：

圍繞三聚氰胺奶粉事件這一案例，設(shè)計了多層級漸進性的探究性課題供學生小組探究（如圖1所示）。多層級漸進性探究型教學模式具有5個特征，分別是“知識由淺入深”“能力培養(yǎng)由低到高”“主動性由弱到強”“探究點的力度逐漸變大”“探究點的難度逐漸增加”，這5個特征以一種循序漸進的方式把學生引入到探究式學習中[5]。

(3A+2ΔtC)e(m+2)=A(4e(m+1)-

e(m))-2ΔtS(m+2)。

(5)

式中：e(m)為m時刻的電場；S(m)為m時刻的源項；Δt為時間步長。為了模擬隨地表起伏的復雜場源，我們將接地導線分解成若干段首尾相連的短導線，每段導線近似為一個電偶極子。每個電偶極子可以表示為

js(r,t)=δ(r-rs)I(t)fdl。

(6)

式中：δ為脈沖函數(shù)；rs為電偶極子的位置；I(t)為電流強度；f為電流方向；dl為電偶極子長度。

本文采用Dirichlet邊界條件[24]，假設(shè)計算區(qū)域外邊界的切向電場分量為0，即

(7)

2 數(shù)值精度驗證

為了檢驗本文算法的精度，首先設(shè)計了如圖1所示的100 Ω·m均勻半空間模型，發(fā)射源長度為1 000 m，中心點坐標為(0 m,-500 m, 0 m)，接收點坐標為(0 m, 500 m, 0 m)，發(fā)射電流20 A。利用非規(guī)則四面體網(wǎng)格將整個計算區(qū)域剖分成249 486個單元。為了精細刻畫場源附近電磁場的劇烈變化，同時提高觀測點處的插值精度，對發(fā)射源和接收點處的網(wǎng)格進行了局部加密。核心計算區(qū)域為3 000 m×3 000 m×500 m，為了滿足Dirichlet邊界條件，在每個方向分別設(shè)置20 km的擴邊。采用我們編寫的非結(jié)構(gòu)時間域有限元正演模擬程序模擬了階躍電流波形dBz/dt和ex響應。計算時間被分成了7段，每段100個時間步長，各段時間步長分別為0.1，0.5，2.5，12.5，62.5，312.5和1 560.0 μs。圖2展示了本文三維時間域有限元正演結(jié)果和一維解析解的對比情況，可以看到三維正演結(jié)果與一維解析解吻合得非常好，最大誤差均在5%以內(nèi)，該結(jié)果有效地驗證了本文算法的準確性。

圖1 均勻半空間模型

圖2 均勻半空間模型正演響應和誤差分析

3 典型山峰和山谷地形對電性源瞬變電磁響應影響特征分析

3.1 接收點處地形對觀測響應的影響

本文首先設(shè)計了如圖3所示的3個半空間模型來分析接收點附近地形對觀測響應的影響特征。山峰和山谷地形的最大起伏均為100 m，地下半空間電阻率為100 Ω·m。發(fā)射源沿平行x軸方向布設(shè)，長度為1 000 m，發(fā)射源的中心點坐標為(0 m,-500 m, 0 m)，發(fā)射20 A的階躍電流。采用面積性觀測方式進行測量，點距和線距均為20 m，共71×71=5 041個測點，同時觀測dBz/dt和ex響應。山峰地形的峰頂和山谷地形的谷底分別位于(0 m, 500 m,-100 m)和(0 m, 500 m, 100 m)處。

a.山峰模型；b.水平地表模型；c.山谷模型。

圖4和圖5展示了斷電后整個觀測區(qū)域內(nèi)不同時刻的dBz/dt和ex響應，圖6和圖7展示了不同偏移距的dBz/dt和ex剖面響應。根據(jù)電性源瞬變電磁擴散理論[27]，在水平地表條件下，當導線源中電流突然關(guān)斷后，在發(fā)射源附近的地下介質(zhì)中會瞬間產(chǎn)生一個與導線源同向的感應主電流，并在四周很大范圍內(nèi)形成回流電流，以維持原有磁場；隨著時間的推移，感應主電流開始向下擴散，且主電流的范圍逐漸擴大；根據(jù)畢奧薩伐爾定律，在發(fā)射源兩側(cè)的赤道方向會產(chǎn)生兩個dBz/dt響應中心，而且中心位置隨著主電流下移逐漸向遠離發(fā)射源方向移動(圖 4f—j)。相比之下，ex響應的分布形態(tài)比較簡單，最大響應始終位于發(fā)射源附近(圖5f—j)。

當存在起伏地形時，觀測響應變得十分復雜，尤其是早期觀測信號出現(xiàn)嚴重畸變。對于早期dBz/dt響應，當響應中心未到達山峰地形時，山峰地形對dBz/dt響應表現(xiàn)為“吸引作用”(圖4a)，在峰頂處出現(xiàn)明顯正異常(圖6b)，同時在峰頂遠離發(fā)射源一側(cè)出現(xiàn)了反號現(xiàn)象(圖4a和圖6c)。隨著時間的推移，反號現(xiàn)象迅速消失，但山峰地形對dBz/dt響應的”吸引作用”依然存在(圖4b)。當dBz/dt響應中心到達并經(jīng)過地形時，山峰地形對dBz/dt響應表現(xiàn)為“排斥作用”(圖4c—e)，產(chǎn)生負異常(圖6b)。然而山谷地形恰恰相反，早期dBz/dt響應中心未到達山谷地形時，山谷地形對dBz/dt響應表現(xiàn)為“排斥作用”(圖4k)，在谷底處出現(xiàn)明顯負異常(圖6j)，同時在谷底靠近發(fā)射源一側(cè)出現(xiàn)了反號現(xiàn)象(圖4k和圖6i)。隨著時間的推移，反號現(xiàn)象同樣迅速消失。當dBz/dt響應中心到達并經(jīng)過谷底時，山谷地形對dBz/dt響應表現(xiàn)為“吸引作用”(圖4m—o)，產(chǎn)生正異常(圖6j)。對于ex響應，不論山峰地形還是山谷地形均未出現(xiàn)變號現(xiàn)象。山峰地形對早期ex響應表現(xiàn)為“吸引作用”(圖5a—b)，產(chǎn)生正異常(圖7b)，而對晚期ex響應起到“排斥作用”(圖5c—e)，產(chǎn)生負異常(圖7b)。山谷地形完全相反，山谷地形對早期ex響應表現(xiàn)為“排斥作用”(圖5k—l)，產(chǎn)生負異常(圖7j)，而對晚期ex響應表現(xiàn)為“吸引作用”(圖5m—o)，產(chǎn)生正異常(圖7j)。圖4e、j、o和圖5e、j、o顯示晚期dBz/dt和ex響應受到地形的影響仍然十分嚴重，而且相比于dBz/dt響應，地形對電場ex的影響更明顯。此外，我們還可以發(fā)現(xiàn)不論是dBz/dt響應還是ex響應，地形對早期和晚期信號的影響規(guī)律具有相反特征，dBz/dt響應在早期甚至出現(xiàn)變號現(xiàn)象?？梢钥闯?，接收點處地形對早期和晚期觀測響應均有嚴重影響，如果在數(shù)據(jù)解釋時不考慮地形的影響，將導致解釋結(jié)果存在巨大偏差。

a—d.山峰模型響應；e—h.水平地表模型響應；i—l.山谷模型響應。

a—d.山峰模型響應；e—h.水平地表模型響應；i—l.山谷模型響應。

上述復雜的地形影響規(guī)律是由觀測點高程變化和地形與周圍介質(zhì)電性差異共同引起的[28]。不同高度接收位置的電磁響應數(shù)值模擬結(jié)果(圖4和圖5)表明，接收點高度的增加將導致早期響應被增強，晚期響應被削弱；而接收點高度降低將產(chǎn)生完全相反的影響。與水平地表模型相比，山峰地形接收點高程明顯高于水平地面，而山谷地形接收點高程低于水平地面，這成為影響地形頂點處電磁響應產(chǎn)生畸變的主要因素。除了影響觀測點高程之外，地形本身還是一個處在地表附近的異常體，其中山峰是空氣中的低阻體，而山谷則是埋藏在地下的高阻體。當電源供電時，在地下空間中會產(chǎn)生與供電電流方向相反的回流電流，由于地形與周圍介質(zhì)間存在電性差，因此在地形兩側(cè)會產(chǎn)生積累電荷。當突然斷電時，地下會產(chǎn)生與供電電流方向相同的主電流，此時地形兩側(cè)的積累電荷也會放電，產(chǎn)生散射電流場。由于山峰地形相當于一個低阻體，地形內(nèi)部散射電流方向與半空間產(chǎn)生的主電流方向相反，因此在山峰地形靠近源一側(cè)產(chǎn)生與半空間同向的感應電壓，使得響應增強，而在遠離源一側(cè)產(chǎn)生與半空間反向的感應電壓，由于早期散射電流場很強導致感應電壓dBz/dt出現(xiàn)變號。同時，根據(jù)電流的連續(xù)性，山頂?shù)乇頃a(chǎn)生散射回流，且與主電流方向相同，因此導致早期電場ex響應增強。同理，對于山谷地形，由于是高阻體，產(chǎn)生的散射電流場與半空間產(chǎn)生的主電流方向相同，因此在山谷地形靠近源一側(cè)產(chǎn)生與半空間反向的感應電壓，產(chǎn)生變號現(xiàn)象，而在遠離源一側(cè)產(chǎn)生與半空間同向的感應電壓，使得響應增強。同時，谷底地表會產(chǎn)生與主電流方向相反的散射回流，導致早期電場ex響應被削弱。到了晚期，積累電荷放電結(jié)束，由其產(chǎn)生的局部畸變也隨之消失，此時接收點高程的影響逐漸占主導地位。

a—e.山峰模型響應；f—j.水平地表模型響應；k—o.山谷模型響應。

a—e.山峰模型響應；f—j.水平地表模型響應；k—o.山谷模型響應。

3.2 發(fā)射源和接收點之間的地形對觀測響應的影響

在實際工作過程中，經(jīng)常會遇到在發(fā)射源和測區(qū)之間存在地形起伏的情況，為了分析該種情況下地形對觀測響應的影響，我們在圖6和圖7中展示了圖3模型中不同偏移距的dBz/dt和ex剖面響應，其中y=500 m剖面恰好經(jīng)過地形中心，以y=850 m和y=1200 m兩條測線來分析發(fā)射源和接收點之間的地形對觀測響應的影響?？梢钥闯霾徽撌莇Bz/dt響應還是ex響應，當觀測點距離地形較近時(圖6c、g、k和圖7c、g、k)，觀測響應會受到很大影響，但是隨著偏移距的增加這種影響迅速減弱(圖6d、h、l和圖7d、h、l)，因此當觀測位置距離地形較遠時，發(fā)射源和接收點之間地形的影響可以忽略。

3.3 發(fā)射源處地形對觀測響應的影響

本文設(shè)計了如圖8所示的發(fā)射源附近地表起伏模型，其中山峰和山谷模型的最大地形起伏均為100 m。半空間電阻率為100 Ω·m。發(fā)射源平行于x軸方向，長度為1 000 m，沿起伏地表敷設(shè)，發(fā)射源中心點的水平坐標均為(0 m,-500 m)，發(fā)射20 A的階躍電流。觀測系統(tǒng)與圖3相同。圖9和圖10展示了斷電后整個觀測區(qū)域內(nèi)不同時刻的dBz/dt和ex響應。從圖9可以看出觀測源處地形對dBz/dt響應形態(tài)并沒有明顯影響，但是對早期響應幅值影響嚴重(圖11a)。其中山峰地形使得早期dBz/dt響應增強，產(chǎn)生正異常(圖11a)，而山谷地形的存在削弱了dBz/dt響應，產(chǎn)生負異常(圖11a)。對于ex響應，山峰地形未改變觀測響應的分布形態(tài)，只改變了響應幅值，產(chǎn)生了正異常(圖11b),而山谷地形完全改變了早期ex響應的分布形態(tài)，并產(chǎn)生了明顯的負異常(圖10k、l、圖11b)。這主要是由發(fā)射源和觀測點相對高程變化以及電磁波傳播路徑改變引起的。對于山谷地形，發(fā)射源高程低于水平地表，感應電磁波到達接收點的路徑增長，而且中間介質(zhì)為導電大地，因此高頻信號被吸收，導致早期響應幅值降低；對于山峰地形，雖然電磁波傳播路徑增長了，然而中間介質(zhì)為空氣，對電磁波的吸收更弱，因此早期信號增強。隨著時間的推移，地形影響逐漸減弱。這是由于早期信號主要集中在地表，且高頻信號占主導部分，具有較高的分辨率，此時對地形十分敏感；隨著時間的推移，感應渦流逐漸向地下傳播，高頻成分迅速減弱，低頻成分占主導地位，分辨率降低，地形的影響逐漸減弱。此外，我們還可以發(fā)現(xiàn)相對于dBz/dt響應，ex響應受發(fā)射源處地形的影響更加嚴重。與接收點處地形相比，發(fā)射源處地形的影響主要集中在早期，對晚期信號影響較小。

a.山峰模型；b.水平地表模型；c.山谷模型。

a—e.山峰模型響應；f—j.水平地表模型響應；k—o.山谷模型響應。

a—e.山峰模型響應；f—j.水平地表模型響應；k—o.山谷模型響應。

圖11 發(fā)射源附近地表起伏模型水平位置(0 m, 500 m)處dBz/dt和ex響應

3.4 發(fā)射源和接收點附近同時存在地形對觀測響應的影響

本文進一步考慮了在發(fā)射源和接收點附近同時存在地形的情況(圖12)，地下電阻率為100 Ω·m的半空間。發(fā)射源和接收點處地形最大起伏均為100 m，發(fā)射和接收裝置與圖3和圖8一致。從圖13可以看出，dBz/dt響應的分布形態(tài)與只有接收點附近存在地形時(圖4)基本一致。然而，受到發(fā)射源和接收點兩處地形的影響，ex響應分布形態(tài)產(chǎn)生嚴重畸變，由于兩處地形的相互耦合導致早期ex響應出現(xiàn)了變號現(xiàn)象；隨著時間的推移，晚期ex響應分布形態(tài)與只有接收點附近存在地形時基本一致(圖14和圖5)?？梢姡邮拯c處地形對晚期觀測信號的影響要強于發(fā)射源處地形的影響。

a.發(fā)射源位于山谷，接收點位于山峰；b.水平地表模型；c.發(fā)射源位于山峰，接收點位于山谷。

a—e.發(fā)射源位于山谷，接收點位于山峰；f—j.水平地表模型；k—o.發(fā)射源位于山峰，接收點位于山谷。

a—e.發(fā)射源位于山谷，接收點位于山峰；f—j.水平地表模型；k—o.發(fā)射源位于山峰，接收點位于山谷。

4 復雜地形影響特征分析

本文最后設(shè)計了如圖15所示的復雜起伏地表模型，該模型為實際地形，地形數(shù)據(jù)來自于陜西省南部。地形最大高差約為450 m，地下半空間電阻率為100 Ω·m，電阻率為1 Ω·m、半徑為200 m的良導球體埋藏在半空間中，球體中心水平位置位于(0 m, 600 m)處，頂面埋深約為110 m。發(fā)射源長度為1 000 m，沿平行x軸方向布設(shè)于y=-650 m處的山谷中，發(fā)射源的中心點水平坐標為(0 m,-650 m)，發(fā)射20 A的階躍電流。采用面積性觀測方式觀測dBz/dt和ex響應。從圖16和圖17可以看出：相比于水平地表模型(圖16k—o和圖17k—o)，復雜地形使得觀測響應產(chǎn)生嚴重畸變，早期的dBz/dt和ex響應均出現(xiàn)了變號現(xiàn)象(圖16a，b，f，g；圖17 a，b，f，g)；隨著時間的推移，地形對dBz/dt響應的影響開始減弱，良導球體的影響逐漸凸顯，即使存在地形的影響，但是仍然可以有效分辨地下良導目標(圖16h—j)；然而，復雜地形對ex的早期和晚期響應均產(chǎn)生了嚴重影響，地形和異常體的相互耦合使得觀測響應產(chǎn)生嚴重畸變，導致無法對地下目標體進行有效分辨，這給后續(xù)的數(shù)據(jù)解釋造成了巨大困難(圖17h—j)。這一對比結(jié)果表明地形對ex觀測信號的影響更加嚴重。

a.起伏地表半空間模型；b.起伏地表良導體模型；c.水平地表良導體模型。

a—e.起伏地表半空間模型響應；f—j.起伏地表良導異常體模型響應；k—o.水平地表良導異常體模型響應。

a—e.起伏地表半空間模型響應；f—j.起伏地表良導異常體模型響應；k—o.水平地表良導異常體模型響應。

5 結(jié)論

本文利用非結(jié)構(gòu)時間域有限元算法成功實現(xiàn)了電性源瞬變電磁復雜地形響應模擬。利用非結(jié)構(gòu)四面體網(wǎng)格進行空間離散，精細刻畫了起伏地形；采用后退歐拉離散格式進行時間離散，放寬了對時間步長的限制；最后利用直接求解技術(shù)求解有限元控制方程實現(xiàn)了非結(jié)構(gòu)時間域有限元正演模擬。

通過模擬大量典型起伏地表模型瞬變電磁響應，我們發(fā)現(xiàn)地形對dBz/dt和ex觀測響應具有嚴重影響，且具有以下特征：

1)地形對早期觀測響應的影響尤為突出，使觀測信號產(chǎn)生嚴重畸變，甚至產(chǎn)生變號現(xiàn)象。

2)相對于dBz/dt響應，電場ex響應受地形響應更加嚴重。

3)接收點處地形對晚期觀測信號的影響要強于發(fā)射源處地形的影響。

4)隨著遠離地形，發(fā)射源和接收點之間的地形對觀測響應的影響逐漸變?nèi)?，當距離足夠遠時可以忽略地形影響。

5)復雜地形與地下目標體相互耦合使得觀測響應十分復雜，導致無法對地下目標體進行有效分辨，給地下目標體的準確探測帶來了巨大困難。

鑒于地形對電性源短偏移距瞬變電磁響應的巨大影響，開展地形矯正工作將是我們下一步的工作重點。