尚曉榮1,2，嶽明鑫1,2,3，楊曉冬1,2，吳小平1,2，李勇3
1 中國科學技術大學地球和空間科學學院
2 中國科學技術大學安徽蒙城地球實體國家野外科學觀測研究站
3 中國地質科學院地球實體地球化學勘查研究所自然資源部地球實體電磁法探測技術重點實驗室
第一作者：尚曉榮，碩士研究所學生，2015年畢業于太原理工大學，獲資源勘查工程專業學士學位；現就讀于中國科學技術大學，攻讀地球實體學專業碩士學位，主要研究方向是時頻域可控源電磁法的三維數值模拟和反演解釋。
導讀：
可控源電磁法已在油氣和礦産資源勘探領域廣泛應用，地形起伏對其勘探效果有很大影響，正常可控源電磁資料處了解釋多是假設平坦地形條件，這會導緻反演異常體的位置和形态發生畸變，甚至産生虛假異常。
本文為探究起伏地形對三維可控源電磁場傳播的影響，提高資料處了解釋品質，基于非結構化網格的矢量有限元算法，開展了數值模拟研究，分别讨論了測點處為山脊地形和山谷地形、發射源處于山脊地形、發射源與測點間為山脊地形等情況下頻率域可控源電磁法的響應特征。并結合實際研究了湘西某複雜形态銅鉛鋅多金屬礦三維地質模型的電磁響應特征，驗證了該算法的可靠性和實用性。
本文主要研究成果一是提供了一種可控源電磁資料三維高精度正演方法；二是讨論的不同地形條件下的三維電磁響應特征，為二維和三維可控源電磁法異常地質解釋提供了理論模型，具有實際應用意義。
------内容提綱------
0 引言
1 數值模拟理論
2 數值模拟算例
2.1 層狀模型
2.2 簡單組合異常體模型
2.3 三維起伏地形模型
2.3.1 測點處山脊地形模型
2.3.2 測點處山谷地形模型
2.3.3 發射源處起伏地形模型
2.3.4 發射源與測點間地形起伏模型
3 複雜金屬礦應用執行個體
4 結論
---------------
0 引言
可控源電磁法(Controlled-source electromagnetic method)采用人工場源激發交流電磁波信号(10-3～105Hz)，通過測量目标區交變電磁場達到研究地下媒體電性結構的目的。由于成本低、勘探深度大、抗幹擾能力強等特點，可控源電磁法在油氣和礦産資源勘探、地熱資源開發及環境工程等領域發揮了重要作用。可控源電磁法包括時間域電磁法和頻率域電磁法。時間域電磁法利用不接地回線源或接地線源向地下發射一次脈沖磁場，在一次場間歇期間，利用線圈或接地電極觀測二次渦流場。頻率域電磁法常見的是可控源音頻大地電磁法，主要使用接地水準電偶源作為人工信号發射源，本文研究的即是頻率域可控源電磁法。近年來，為适應“深地”國家重大探測戰略中電磁精細結構探測的需求，可控源電磁儀器的研發及資料采集技術獲得了迅猛發展，資料的反演解釋正逐漸邁入三維實用化階段。與實際需求相比，三維可控源電磁資料的反演技術尚未成熟，計算效率和準确性均有待提高，特别是考慮起伏地形的三維反演尚不多見。正演是反演的基礎，反演解釋的準确性和計算效率很大程度上依賴于正演算法，是以研究可控源電磁資料的三維高精度正演方法十分必要。
三維可控源電磁法常用的數值模拟方法主要包括積分方程法、有限差分法、有限體積法以及有限單元法]等。任政勇等基于有限元方法提出一種新的非結構化網格局部節點加密方法，實作了複雜模型完全非結構化四面體全自動剖分；Puzyrev等基于節點有限元方法，采用複雜疊代求解器和基于稀疏近似逆的有效預調節器對大型稀疏線性方程組進行求解，實作了各向異性媒體中三維可控源電磁場的正演模拟；蔡紅柱等提出一種基于非結構化網格的海洋電磁有限單元正演算法，該算法适用于井中電磁、航空電磁、環境地球實體等非均勻且各向異性媒體的電磁感應基礎研究；湯文武等比較了基于節點有限元與矢量有限元的三維可控源電磁正演，結果表明相同條件下矢量有限元法的計算精度更高，但速度較慢；He等基于棱邊有限元方法實作了各向異性媒體中三維張量可控源電磁法的正演。
正常可控源電磁資料解釋多是基于平坦地形，這會導緻反演異常體的位置與形态常發生畸變，甚至産生虛假異常。随着近年來地質勘查目标和野外實地勘測場景越來越複雜，學者們開始研究起伏地形下的三維電磁正、反演。Nam等采用棱邊有限元方法研究了大地電磁法在三維地形下的響應特征；Lin等分析了可控源音頻大地電磁測深資料的三維地形效應；殷長春等采用自适應非結構有限元法進行了帶地形的時間域海洋電磁法正演。積分方程法、有限差分法和有限體積法大多是基于結構化網格實作的，其中有限單元法對網格剖分比較靈活，非結構化網格更适于模拟真實地形起伏。基于非結構化網格有限元法的二維和三維海洋電磁問題已經取得了一些研究成果。相比之下，起伏地形下的陸地三維可控源電磁數值模拟研究相較滞後，目前關于起伏地形對各個場分量的影響特征的研究尚不多見。
是以，本文基于矢量有限元算法開展了起伏地形下的三維可控源電磁法多分量數值模拟研究。采用非結構化四面體網格對計算區域進行剖分，以準确模拟起伏地形，并基于棱邊基函數實作電磁場的插值和求解。首先，建構層狀媒體模型，将其三維矢量有限元計算結果與一維解析解進行對比，驗證算法的正确性；其次，建構簡單異常體模型，驗證算法對探測低阻異常體的有效性；然後，分别讨論了測點處為山脊地形和山谷地形、發射源處于山脊地形、發射源與測點間為山脊地形等情況下頻率域可控源電磁法的響應特征；最後，研究了複雜三維地質模型的電磁響應特征，驗證該算法的實用性。
1 數值模拟理論
從麥克斯韋方程組出發，假設時諧因子為e-iωt的平面電磁波入射到均勻媒體中，則計算區域電磁場滿足
式中：E為電場強度；H為磁場強度；ω為角頻率；μ為媒體磁導率；σ為電導率；ε為相對介電常數；J為外加電流密度。
由式(1)和式(2)可得整個計算區域内可控源電磁法滿足的電場波動方程
對式(4)采用基于棱邊的非結構四面體矢量有限元法進行求解，電場矢量場可表示為六個矢量棱邊插值基函數的線性組合。非結構四面體單元如圖1所示。
圖1 非結構四面體單元ｅ示意圖
ｅ1～ｅ6為棱邊編号
對起伏地形模型的正演方程組(式(3))求解時，右端項J不僅包含地下電阻率異常體引起的散射電流密度，還包括起伏地形産生的散射電流密度。使用非結構有限元法剖分模型，同時對地形起伏區域進行網格加密，以更好地對地形進行拟合。
将單元矩陣進行組裝，可得到大型稀疏複線性方程組
KE=S (12)
式中：K=Ae-iωWe-ω2Ce，為正演系數矩陣；S=iωSe，反映源的作用。這裡的E即為待求電場強度矢量。通常，K是一大型對稱稀疏矩陣，采用CSR(Compressed Sparse Row，稀疏矩陣按行壓縮)格式進行存儲以減小存儲空間。
最後，使用直接求解器Pardiso求解控制方程，可得到單元網格棱邊上的電場強度和磁感應強度分量。
2 數值模拟算例
2.1 層狀模型
為驗證算法的正确性，建構一個三層地層的層狀模型。對模型的電、磁響應特征進行分析，并根據模拟結果與一維解析解的相對誤差(計算公式為(數值解-解析解)/解析解)分析本文算法的精度。
層狀模型的空間尺寸為100km×100km×100km，如圖2a所示。發射源位于y=-10km處，測點位于原點，發射頻率為0.1～5000Hz，共10個頻點：0.1、0.5、1、5、10、50、100、500、1000、5000Hz。用非結構四面體網格剖分生成的最終網格分布見圖2b。采用Key提出的一維可控源電磁正演算法與本文基于非結構矢量有限元算法進行模拟，并對比響應曲線，結果見圖3。由圖可見，本文算法模拟結果與一維解析解吻合較好，電場的分量Ex實部和虛部的數值解與解析解的相對誤差的絕對值分别小于1%和2%，精度較高，驗證了本文算法的正确性。
圖2 層狀模型幾何結構(a)和網格剖分示意圖(b)
圖3 層狀模型Ex分量實部和虛部數值解與解析解對比(左)及其相對誤差(右)
在保證精度的前提下，分别測試了層狀地層模型不同剖分網格單元數情況下數值模拟所需要的時間，并統計了數值解與解析解的誤差，結果見表1。
表1 層狀模型數值模拟時間和精度統計
從表1可以看出，網格數增加約1倍時，數值模拟耗時大約增加2倍；網格數大于33萬後，誤差下降較慢。是以，在綜合考慮耗時與精度的情況下，本文選擇33萬網格數進行後文模型剖分。
2.2 簡單組合異常體模型
為了分析本文算法對異常體電阻率的敏感性，建構兩組簡單組合異常體模型，背景為均勻半空間。模型包含4個相同的低阻(10Ω·m)/高阻(1000Ω·m)異常體，異常體尺寸為500m×500m×200m(圖4)。模型空間為100km×100km×100km；發射源點位于y=-3km處；測量點位于x軸的-1～1km、y軸-1～1km所圍區域的地面，x、y方向的點距均為200m；發射頻率同前文層狀模型。對接收點附近進行了網格加密，在遠離觀測點處，網格剖分比較稀疏，模型最終網格剖分如圖5所示。
圖4 簡單組合異常體模型
圖5 組合異常體模型網格剖分示意圖
左：x=0剖面；右：z=100m剖面。黃色方塊是異常體
分别對低阻異常體和高阻異常體兩個模型進行模拟，結果見圖6。從圖中可以看出，低阻異常體産生的電場強度異常較高阻異常體明顯，且不同頻率時的相對異常(絕對值)也不同。
圖6 簡單組合異常體模型的Ex振幅曲線(左)及其相對異常曲線(右)
為了進一步分析背景模型、包含高阻異常體和低阻異常體模型的電場特征，繪制了0.5Hz的Ex振幅圖(圖7)。可以看出，異常體的賦存區域在有、無異常體的情況下差異較大；對比低阻和高阻異常體的情況可知，頻率域可控源電磁法對于低阻異常體更敏感。
圖7 0.5Hz時組合異常體模型地表Ex分量振幅圖
(a)不存在異常體(背景模型)；(b)存在低阻異常體；(c)存在高阻異常體；(d)圖b與圖a的相對異常；(e)圖c與圖a的相對異常
2.3 三維起伏地形模型
2.3.1 測點處山脊地形模型
建立一個含有山脊地形的三維簡單異常體模型(圖8)。山脊頂部長、寬均為800m，底部長、寬均為1600m，高度為400m，其中心點位于(0，0，200m)。假設大地背景電阻率為100Ω·m，異常體為不規則形狀(圖8中的藍色區域)，厚度為200m，長、寬均為500m，電阻率為1Ω·m，中心點位于(0，0，750m)。采用非結構化四面體單元對模型進行網格剖分，對發射源和接收點處進行網格加密，模型網格剖分結果見圖8。采用線源，源中點水準位置為(3km，-3km)，發射頻率同前文層狀模型采用的頻點；接收點位于地面-1～1km(x)、-1～1km(y)區域内，x、y方向測點間距均為200m。
圖8山脊地形低阻異常體模型網格剖分示意圖
(a)水準地形的低阻異常體模型(M01)；(b)含山脊地形的低阻異常體模型(M02)；(c)含山脊地形的無異常體模型(M03)
圖9為圖8所示三個模型在測點(0，0，0)處的電場強度分量(Ex)和磁感應強度分量(By)振幅計算結果。由圖可見，模型M02與模型M03的正演Ex和By振幅曲線較一緻，而M01模型與二者相差較大，說明山脊地形會削弱異常體的存在，同時也說明正演時如果忽略山脊地形，會對正演結果産生較大影響。對比圖9b與圖9d可以發現，山脊地形對Ex振幅的影響大于對By振幅的影響。
圖9山脊地形模型Ex和By振幅正演曲線及相對異常
(a)Ex振幅；(b)Ex振幅相對異常；(c)By振幅；(d)By振幅相對異常
圖10為圖8中三個模型在頻率為1Hz時沿x方向(y=0)的正演分量Ex和By的振幅及其相對異常。可以看出，山脊根部位置的正演Ex和By振幅均高于水準地形下的場分量幅值，山脊頂部位置的正演Ex和By振幅均低于水準地形下的場分量幅值，這是由于在地形高處(即山脊頂部)，供電電流是發散的，會呈現低電流密度，而在地形相對低處(即山脊根部)，供電電流是收斂的，呈現高電流密度。此現象說明山脊地形會削弱異常體所産生的異常。
圖10 1Hz時山脊地形模型Ex和By振幅曲線
(a)Ex振幅；(b)Ex振幅相對異常；(c)By振幅；(d)By振幅相對異常
2.3.2測點處山谷地形模型
建立一個含有山谷地形的簡單低阻異常體模型(圖11)。山谷頂部長、寬均為1600m，底部長、寬均為800m，高度為400m，其中心點位于(0，0，200m)。異常體參數與前文的山脊模型低阻體相同。采用與山脊地形模型相同的加密剖分網格，在相同位置布設發射源、接收點，采用同樣的發射頻率。模型網格剖分結果見圖11。
圖11 山谷地形模型網格剖分示意圖
(a)低阻異常體模型(M04)；(b)無異常體模型(M05)
圖12為水準地形的低阻異常體模型(M01)和圖11中的山谷地形模型M04和M05的Ex和By正演振幅曲線及相對異常曲線。由圖可見，M01、M04和M05三個模型的電、磁場強度振幅曲線差異較大；與山脊地形的By正演振幅曲線(圖9c和圖9d)相比，山谷地形低阻體模型的差異更大，可見頻率域可控源電磁法的電場強度、磁感應強度分量受山谷地形的影響更大。
圖12 山谷地形模型Ex和By振幅正演曲線及相對異常
(a)Ex振幅；(b)Ex分量振幅相對異常；(c)By振幅；(d)By振幅相對異常
圖13為山谷模型在頻率為1Hz時沿x方向(y=0)的正演Ex和By振幅及振幅相對異常。可以看出，在山谷兩側Ex和By振幅均低于水準地形下的場分量幅值，在山谷底部則均相反。這一特征與山脊地形的計算結果相似，這是由于在地形高處，地表下供電電流發散，呈現低電流密度；在地形相對低處(即山谷底部)，供電電流收斂，呈現高電流密度。此現象說明山谷地形會增強異常體産生的異常。
圖13 1Hz時山谷地形模型Ex和By振幅沿x方向變化曲線
(a)Ex振幅；(b)Ex振幅相對異常；(c)By振幅；(d)By振幅相對異常
2.3.3 發射源處起伏地形模型
建構一個發射源處為起伏地形的模型(圖14)。源的中點位于起伏地形區域的中心位置(0，-8000m，-395m)，山脊頂部長、寬均為1000m，底部長、寬均為1600m，高度為400m。其他參數設定均與測點處山脊地形模型參數相同。發射源處地形平坦模型(M06)和地形起伏模型(M07)的網格剖分結果見圖14。
圖14 發射源處地形平坦模型(M06，上)和地形起伏模型(M07，下)網格剖分示意圖
圖15為模型M06和M07在測點(0，0，0)處的電、磁場分量Ex和By振幅随頻率變化曲線及相對異常。從圖15a可以看出，頻率高于100Hz時，模型M07的Ex振幅高于模型M06；頻率低于100Hz時，情況相反。這是因為頻率越低，探測深度越大，發射源處的地形為山脊地形，相對增加了媒體厚度；而且，對于高阻異常體，在相同的頻率下，電場衰減較快，振幅衰減也更快。高頻信号的探測深度較小，發射源處向下傳導的電流在相同深度接收點處的電場信号強度會更大。另外，從圖15b可以看出，兩個模型的Ex相對振幅差較大，最大值達48%。由圖15b和圖15d可見兩個模型的By分量相對異常與Ex分量有相似特征。
圖15 發射源處地形起伏模型Ex和By振幅正演曲線及相對異常
(a)Ex振幅；(b)Ex分量振幅相對異常；(c)By振幅；(d)By振幅相對異常
2.3.4 發射源與測點間地形起伏模型
建立一個源點與測點間地形起伏模型(圖16)，分析這類地形對可控源電、磁場響應的影響特征。起伏地形山脊頂部長、寬均為800m，底部長、寬均為1600m，高度為400m，中心點為(0，-2000m，0)。發射源中點位于地表(8km，-8km)處。其他參數設定同2.3.1節模型。圖16為發射源和測點間地形平坦模型(M06)和地形起伏模型(M08)網格剖分圖。
圖16 發射源與測點間地形平坦模型(M06，上)和地形起伏模(M08，下)網格剖分示意圖
圖17為模型M06和M08在測點(0，0，0)處的電、磁場分量Ex和By振幅曲線及振幅相對異常曲線。從圖17a和圖17b可見，兩個模型的Ex振幅較接近，差異較小。另外，發射源與起伏地形的距離、地形起伏程度、地形與異常體的距離等因素也會影響電、磁場。圖17c和圖17d是兩個模型By分量振幅曲線和相對異常曲線，可見By分量振幅相對異常曲線與Ex分量有相似規律。
圖17 發射源與測點間地形起伏模型Ex和By振幅正演曲線及相對異常
(a)Ex振幅；(b)Ex分量振幅相對異常；(c)By振幅；(d)By振幅相對異常
從上述分析可知，由于地形的存在，尤其是發射源之下起伏地形和測點處的起伏地形，采用頻率域可控源電磁法對地下異常體進行探測時，電、磁場分量均發生了畸變，出現了虛假異常。是以在實際資料的處了解釋中，需要考慮地表地形對可控源電磁法電、磁場的影響。
3 複雜金屬礦應用執行個體
為了驗證本文算法在實際應用中的可靠性，以湘西某銅鉛鋅多金屬礦床的成礦模式為基礎建立圖18所示地質模型，所包含的主要岩(礦)石的電阻率資訊如表2所示。
圖18 湘西銅鉛鋅多金屬礦床成礦模式圖
表2 湘西銅鉛鋅多金屬礦床主要岩(礦)石電阻率統計
根據表中岩(礦)石電阻率範圍，取其平均值，給出圖19所示的金屬礦床模型二維電阻率填圖結果及其非結構四面體模型網格剖分結果。對原模型四周區域進行延拓，其中空—地邊界沿原模型邊界兩端分别向兩側延伸，而延拓區域的地層電阻率設為與原模型相鄰區域的電阻率一緻。延拓模型空間範圍為40km×40km×40km，包含四個長度為1km的線源，發射電流為1A，源的中點坐标分别為(0，-3000m，105m)、(-3000m，0，115m)、(0，3000m，162m)和(3000m，0，115m)(圖20)。模型中異常體分布區域為-750m≤x≤750m，-750m≤y≤750m，0≤z≤400m，低阻目标體為圖19中填圖顔色為藍紫色的小規模片狀塊體及青綠色的塊狀異常體。測點沿地形起伏均勻分布，地表範圍為-750m≤x≤750m，-750≤y≤750m，測點間距為100m。發射頻率為0.1～5000Hz，頻點分布同2.1節層狀模型。
圖19 湘西銅鉛鋅多金屬礦床模型二維電阻率填圖(上)及其非結構四面體單元剖分局部示意圖(下)
圖20 湘西銅鉛鋅多金屬礦床延拓模型發射源布置示意圖
分别正演計算有、無異常體情況下，金屬礦床模型中心點在x=0、y=-200m處Ex振幅随頻率的變化曲線，結果見圖21a，二者的相對異常見圖21b。由于模型深部電阻率均一，是以低頻時Ex振幅較穩定，且相對異常基本保持為常數；在高頻段，由于地質異常體分布複雜，圖21b中的相對異常曲線出現了明顯的波動。
圖21 金屬礦模型偏移距x=0、y=-200m處Ex振幅随頻率變化曲線(a)及相對異常曲線(b)
選取正演計算的頻率1Hz時的Ex和By振幅資料，繪制圖22所示的等值線圖。可以發現，Ex和By振幅異常區域均與低阻地層的平面位置大緻吻合(圖22中虛線橢圓區域)；同時，Ex和By振幅相對異常圖中相對異常較大的區域(圖22中實線橢圓區域)也與實際低阻異常體的賦存區域存在很強的相關性。
圖22 金屬礦模型頻率1Hz時正演計算的地表位置Ex(上)和By(下)振幅等值線圖及相對異常圖
(a)包含異常體；(b)不包含異常體；(c)圖b與圖a的相對異常
4 結論
本文基于非結構四面體網格和矢量有限元法實作面向複雜地質模型的可控源電磁法三維多分量數值模拟。通過對比層狀模型的數值解與一維解析解，以及組合異常體模型的電、磁場強度響應特征分析，驗證了本文算法的正确性和有效性，說明了模型剖分方案的合理性。通過測點處山脊和山谷地形模型、發射源下山脊地形模型、發射源和測點間山脊地形模型的數值模拟，說明了地形起伏對可控源電、磁場會産生顯著的影響。是以，在地形複雜的區域開展可控源電磁工作必須考慮地形的影響。實際複雜帶地形金屬礦模型的計算結果也進一步驗證了該算法的實用性和可靠性。
原文來源：尚曉榮，嶽明鑫，楊曉冬，吳小平，李勇.起伏地形對三維可控源電磁響應的影響研究.石油地球實體勘探，2022，57(1)：222-236.