電性源短偏移距瞬變電磁地形影響特征分析

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

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

0 引言

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

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

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

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

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

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

(1)

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

(2)

(3)

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

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

圍繞三聚氰胺奶粉事件這一案例，設計了多層級漸進性的探究性課題供學生小組探究（如圖1所示）。多層級漸進性探究型教學模式具有5個特征，分別是“知識由淺入深”“能力培養由低到高”“主動性由弱到強”“探究點的力度逐漸變大”“探究點的難度逐漸增加”，這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)

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

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

(7)

2 數值精度驗證

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

圖1 均勻半空間模型

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

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

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

本文首先設計了如圖3所示的3個半空間模型來分析接收點附近地形對觀測響應的影響特征。山峰和山谷地形的最大起伏均為100 m，地下半空間電阻率為100 Ω·m。發射源沿平行x軸方向布設，長度為1 000 m，發射源的中心點坐標為(0 m,-500 m, 0 m)，發射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展示了斷電后整個觀測區域內不同時刻的dBz/dt和ex響應，圖6和圖7展示了不同偏移距的dBz/dt和ex剖面響應。根據電性源瞬變電磁擴散理論[27]，在水平地表條件下，當導線源中電流突然關斷后，在發射源附近的地下介質中會瞬間產生一個與導線源同向的感應主電流，并在四周很大范圍內形成回流電流，以維持原有磁場；隨著時間的推移，感應主電流開始向下擴散，且主電流的范圍逐漸擴大；根據畢奧薩伐爾定律，在發射源兩側的赤道方向會產生兩個dBz/dt響應中心，而且中心位置隨著主電流下移逐漸向遠離發射源方向移動(圖 4f—j)。相比之下，ex響應的分布形態比較簡單，最大響應始終位于發射源附近(圖5f—j)。

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

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

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

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

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

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

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

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

3.3 發射源處地形對觀測響應的影響

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

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

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

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

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

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

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

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

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

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

4 復雜地形影響特征分析

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

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

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

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

5 結論

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

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

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

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

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

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

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

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