復雜地形結構下瞬變電磁三維正演模擬?以大西洋洋中脊TAG熱液區為例

聶佐夫，陶春輝,2,3*，沈金松，朱忠民

(1.上海交通大學 海洋學院，上海 200240；2.自然資源部第二海洋研究所，浙江 杭州 310012；3.自然資源部海底科學重點實驗室，浙江 杭州 310012；4.中國石油大學 (北京) 地球物理學院，北京 102249)

1 引言

自1972年，多國政府和大學科學家合作開展美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration,NOAA） TAG（Trans-Atlantic Geotraverse）項目以來，位于26°N大西洋洋中脊的TAG熱液區一直是研究熱液活動的焦點[1?2]。熱液系統循環過程中形成的多金屬硫化物，是一種極具經濟價值的礦產資源，近年已成為海底資源研究的一大熱點[3]。目前，已發現的海底硫化物礦體的規模大小差異較大，由幾十米到幾百米不等，因此需要一種行之有效的探測方法來確定海底多金屬硫化物堆積體的位置、深度、規模，進而有效地評價其資源量及其經濟價值，為后續開采工作打下基礎。根據前人在有熱液活動的TAG丘體[4]的鉆探結果以及在已停止熱液活動的區域[5]的鉆探取芯結果分析可知，硫化物礦體圍巖為導電性差的新鮮或蝕變玄武巖，且其導電性與巖石孔隙中的海水鹽度相關。而硫化物礦體中存在較高含量的金屬元素，它們與圍巖存在顯著的導電性差異，這為海底電磁方法探測熱液多金屬硫化提供了物質基礎。

經過半個世紀的發展，海洋電磁在理論和裝備上都有了長足的進步。20世紀80、90年代，在海洋瞬變電磁研究中比較有代表性的是Edwards研究團隊，1986年Edwards和Chare[6]計算了上覆導電海水的均勻半空間電偶極子瞬變電磁二次場響應，這一項研究成為了海洋瞬變電磁法勘探的理論基礎。次年，Cheesman等[7]、Everect等[8]在此基礎上總結出水平電偶極子和磁偶極子兩種設備十分適合用來開展瞬變電磁近底勘探的結論，隨后進一步開展了海試工作，論證了其在海底電導率勘探上的可行性。最近幾年，海洋電磁法探測海底熱液多金屬硫化物堆積體已成為研究的焦點。2012年，Swidinsky等[9]結合實際海底礦床勘探建立模型，推導同點裝置一維層狀介質正演公式，結果分析表明，在一定半徑線圈及電流條件下，同點探測系統是確定導電礦床厚度及其上覆蓋層厚度的有效工具。2015年，德國GEOMAR設計研制了一套瞬變電磁重疊回線裝置MARTEMIS，并在地中海海域已知硫化物堆積體進行了初步測試，在已知的礦體附近，探測到了瞬變電磁異常的存在[10]。2010年，經自然資源部第二海洋研究所組織、湖南五維地質科技有限公司和北京先驅高技術開發公司共同研發出了海洋瞬變電磁重疊回線裝置，并于南大西洋13.2°S熱液區開展試驗，探測到明顯異常[11]。這一裝置在2012年的中國大洋26航次于大西洋洋中脊TAG熱液區進行了探測，得到熱液區內TAG丘體的瞬變電磁響應，并利用視電阻率反演結果，將海底多金屬硫化物礦體近似看作兩層空間中的T型異常體[12?13]。最近幾年在正演方面，國內對于有限元方法的研究也取得了一定進展。2016年，李瑞雪等[14]采用全空間矢量有限元法，對硫化物礦體進行矩形單元剖分，從頻率域出發，計算出深海熱液硫化物礦體的瞬變電磁響應，且與物理模擬結果具有良好的一致性。2017年，Zhang等[15]直接在時間域求解，實現了基于自適應有限元法的航空瞬變電磁三維正演模擬，并分析了地形因素對航空瞬變電磁響應的影響。2019年，殷長春等[16]實現了時間域海洋電磁三維正演模擬，改善了網格質量，提高了計算精度，并分析了在地形隆起與下凹時，海底有無塊狀導體時的近底電磁響應。2020年，彭榮華等[17]首次考慮海洋可控源電磁探測中的各向異性因素，利用三維正演模擬方法研究海底硫化物堆積體響應，發現各向異性因素會引起海洋可控源電磁結果失真從而不利于準確的數據解釋。

TAG熱液區擁有豐富的地質和礦體產狀的背景資料，是目前研究程度最高的活動熱液區之一，在TAG丘體開展瞬變電磁三維正演研究，分析空間電磁場的分布特征和接收裝置的響應特征，有利于我們了解儀器狀態、海底地形以及礦體分布對瞬變電磁二次場響應的影響。結合該區域開展的海洋瞬變電磁探測資料，可以檢驗瞬變電磁方法在探測海底多金屬硫化物堆積體方面的可行性與有效性。

同時，由于目前海洋瞬變電磁實際數據處理仍然以一維反演拼接成擬二維剖面為主[18]，但利用一維層狀介質模型的瞬變電磁響應來解釋復雜地形情況下的響應會存在一定假象[19]，因此有必要采用高維的反演方法，并充分考慮地形起伏、探測裝置測線軌跡以及裝置的姿態等因素，從而得到更為合理的電導率剖面。因此，本研究可為未來海底多金屬硫化物堆積體電導率剖面的精確反演工作打下基礎。

2 TAG熱液區地質背景

TAG熱液區位于大西洋洋中脊26°08′N附近，熱液區面積約為 5 km×5 km，水深范圍為 3 430～3 670 m（圖1）。由1處存在熱液活動的TAG丘體和多處已停止熱液活動的多金屬硫化物丘體組成[2]。其中，TAG丘體存在高溫熱液噴口（溫度最高達363℃），已形成直徑為200 m，高為50 m的環形熱液產物堆積體，活動黑煙囪體高約12 m[20]。大洋鉆探計劃對TAG丘體進行了一系列鉆孔采樣[4]，在丘體周圍的5個區域（TAG-1至TAG-5）共計鉆孔17個，最深鉆孔達125 m。通過對鉆孔中取回的巖芯進行一系列地球化學分析，在海山堆積體內發現了高品位的Cu、Zn、Ag和Au多金屬硫化物[21]，并得到 TAG 丘體的內部結構[4,21?22]：丘體表層為硫化物角礫，下伏富含金屬元素的黃鐵礦、二氧化硅與硬石膏角礫，最深處為金屬含量稍低的脈狀礦化硅化玄武巖。對TAG熱液區內多處硫化物丘體采集的硫化物進行物性測量工作，樣品包括硫化物和作為參照的玄武巖，測得的部分結果如表1所示[23?24]。可見，即使硫化物孔隙度低于玄武巖時，也具有更高的電導率，說明硫化物并不符合基于電解質導電的阿爾奇定律[25]。硫化物因其導電形式更類似于半導體從而具有高電導率[26]。

表1 TAG熱液區樣品及部分陸地玄武巖樣品的電導率測量結果[23?24]Table 1 Conductivity measurement results of sulfide samples from the TAG hydrothermal field and various onshore basalt samples[23?24]

圖1 大西洋洋中脊TAG熱液區高分辨率地形圖（a）以及TAG丘體內部地質結構（b）Fig.1 High-resolution bathymetry map of the TAG hydrothermal field at the Mid-Atlantic Ridge (a) and the geological structure of the TAG mound (b)

3 瞬變電磁響應的數值模擬方法

對于瞬變電磁場計算，我們采用了庫倫規范準則下磁矢量勢?電標量勢與自適應有限元相結合的方法[27]，在此僅給出算法相關的基本關系。在各向同性海水和海底地層中，瞬變電磁場滿足如下Maxwell方程：

式中，E為電場強度（單位：V/m）；H為磁場強度（單位：A/m）；J為傳導電流密度（單位：A/m2）；Je為外加電流密度（單位：A/m2）；μ和 σ 分別為介電常數和電導率。

由公式（1）和公式（2）可得到電場強度E的雙旋度方程，公式為

利用亥姆霍茲分解理論，矢量場E可變為

式中，A為矢量勢；?為標量勢，將公式（4）代入公式（3），得到下式：

同時，將公式（4）代入連續性方程?×σE=??×Je可得公式：

為求解方程（5）和方程（6），考慮各向同性介質的Dirichlet邊界條件，對于求解域Ω和閉合邊界Γ，可得：

式中，n為單位法向量。本文采用Galerkin加權余值法導出方程（5）、方程（6）的弱解形式，并用 Delaunay非結構化四面體單元剖分，四面體單元內各邊與節點的矢量位與標量位可表示為[16]

圖2 四面體單元自由度（修改自參考文獻 [17]）Fig.2 Degrees of freedom associated with a tetrahedral element (modified from reference [17])

建立有限元方程如下求解。得到 A ,?：

式中，系數矩陣元素和源項的計算關系可參考文獻 [17,27?28]，而對于時間步的離散使用后向歐拉法[29?30]。相應地，磁場和電場的計算關系為

4 TAG丘體電性參數的三維建模

本文利用高精度多波束測深數據構建TAG丘體的三維地形圖[31]，為得到海底以下的電導率分布信息，參考大洋鉆探及地球化學探測結果，建立如圖3所示的三維地電模型。模型核心部分為400 m×400 m×300 m的三維結構，其中TAG丘體的尺寸約為180 m×180 m×150 m。模型核心部分外加厚 1 km 的邊界層以確保電磁波衰減，并在層邊界添加Dirichlet邊界條件。根據鉆孔數據確定的TAG丘體內部的電導率分布如圖3c和圖3d所示，丘體內部高品位的黃鐵礦電導率為60 S/m，丘體表層的硫化物角礫電導率為30 S/m，丘體深部為電導率稍低的脈狀礦化硅化玄武巖，電導率約為20 S/m。正演模擬中海水電導率設為3 S/m（海水溫度為2～3℃），圍巖玄武巖電導率設為0.1 S/m。

在瞬變電磁正演模擬中，發射源采用圓形重疊回線裝置，半徑為10 m。發射線圈中施加20 A的階躍電流，為觀測丘體上方的瞬變電磁響應，接收關斷后感應電動勢 ? B/?t 隨時間的衰減，接收間距為10 m，測線如圖3a中紅線所示。

圖3 TAG 丘體所在區域三維電導率分布模型Fig.3 3D conductivity model for the TAG mound

模型建立后，采用非結構化網格（自由四面體）對三維模型進行離散化。為了保證正演模擬的精度，應對發射源和地形變化大、電導率產生差異的界面附近使用局部更為精細的網格。為提高運算效率，根據源的位置對網格進行局部剖分，利用局部網格替代網格數量龐大的全局網格，以減少系統矩陣的自由度。采用PARDISO直接求解器來求解方程組，該求解器是目前最快的線性稀疏矩陣求解方法之一[32]。

5 電磁場響應的特征分析

5.1 電磁場時空分布特征

為了解電磁場在該地質模型中的傳播規律，確定有效信號的觀測范圍，筆者正演計算了瞬變響應在電流關斷后的時空分布。在TAG丘體模型西南角離底20 m處建立重疊回線裝置，線圈半徑為10 m，電流為20 A，通過正演得到TEM脈沖響應? Bz/?t在海底地形上的擴散過程如圖4所示。當發射電流關斷后，變化的磁場會在硫化物丘體內部產生感應渦流，它會隨觀測時間不斷向下和向外擴散且幅值衰減形成“煙圈效應”[33–34]，由于硫化物礦體與圍巖的電導率差異巨大，大多數渦流集中于TAG丘體區域。

圖4 發射線圈關閉后4個時間步長上TAG丘體地形處? Bz/?t擴散圖Fig.4 Forward model of the diffusion of ?Bz/?t across the bathymetry of the TAG mound area at four timesteps after transmitter turn-off

為驗證瞬變電磁法對高導硫化物的探測能力，本文同時也對純玄武巖基底的海底進行了模擬，即把原模型中具有導電性的硫化物角礫、黃鐵礦和脈狀礦化硅化玄武巖都替換為高阻玄武巖，將其正演響應視為背景場，并將含硫化物的異常場與背景場做歸一化，對? Bz/?t歸一化后取對數后的結果如圖5所示，該結果可以指示瞬變電磁響應對高導體的可探測性。在 1×10?5s、3.162 3×10?4s和 1×10?3s時刻，線圈上方的海水因電導率一致而沒有產生異常，但下方TAG丘體處的高導體導致了明顯的 ? Bz/?t差異。同時由于線圈位于丘體西南側（剖面左側），渦流系統經過高導體后，導致整個剖面右側的高阻圍巖響應也出現了明顯的偏差，而兩次正演模擬在剖面右側的差異很小。

圖5 發射線圈關斷后，由背景場歸一化后的異常場，其中第二列圖為西南?東北方向的豎直剖面Fig.5 Forward model of the normalized ?Bz/?t for a homogeneous seafloor model vs.the model which includes the conductive target after transmitter turned off,graphs at the second column are the SW-NE sections of the TAG mound area

5.2 實測瞬變電磁響應的擬合與影響因素分析

5.2.1 TAG 熱液區瞬變電磁數據采集

海上瞬變電磁測量采用重疊回線裝置[12]，線圈為1.95 m×0.75 m 的矩形多匝回線，線圈匝數 40 匝，等效面積為 58.8 m2，發送電流為 20 A 的雙極性方波，占空比為50%，方波頻率為0.625 Hz。海上作業時，科考船以 1～2 kn（0.5～1 m/s）的速度拖曳瞬變電磁拖體勻速前進，通過調節光纜的長度使拖體保持在距海底上方 30～50 m。

2012年的大洋26航次第2航段對海上瞬變電磁裝置進行了測試，結合TAG熱液區礦體位置對測線進行設計，利用深海瞬變電磁雙拖體裝置采集海底多金屬硫化物礦體的響應數據，共完成了6條測線[13]。其中EM02測線由西南向東北從TAG丘體正上方經過，并經過熱液停止活動的MIR區，測線總長為2.67 km。該測線數據在這兩處已知區域均發現明顯的瞬變電磁異常，具有良好的可靠性，因此本文以EM02測線為例來說明實際數據采集中存在的影響因素。

測線位置及測量得到的瞬變響應如圖6所示，圖6b和圖6c分別對應TAG丘體和MIR區上方的瞬變電磁剖面，該數據已做中值濾波處理，并利用源歸一化來剔除線圈大小與發射電流對響應的影響。在此基礎上，為消除儀器本身對瞬變電磁響應的影響從而突出異常的相對變化，利用實測與模擬的海水背景響應對實測數據進行歸一化處理。可以看出，該測線所測得的電壓幅值存在一定的震蕩，在已知的硫化物礦體上方的瞬變電磁響應有明顯的隆起，且TAG丘體處的響應呈“雙峰”狀。為探究其原因，有必要對數據采集過程中可能影響實測數據的各個因素進行研究。

圖6 TAG丘體與MIR區周邊多波束地形圖及大洋26航次瞬變電磁EM02測線（a）、TAG丘體瞬變電磁歸一化響應平剖圖（b）和MIR區瞬變電磁歸一化響應平剖圖（c）Fig.6 The multibeam bathymetry of the area contains the TAG mound and the MIR zone with the EM02 survey line from COMRA 26th Cruise (a),normalized section of TEM response at the TAG mound (b),and normalized section of TEM response at the MIR zone (c)

5.2.2 瞬變電磁響應擬合與分析

本節模擬了重疊回線裝置的近底拖曳測量方式，離底高為20 m時得到的瞬變電磁響應以及高導異常點對應位置如圖7所示。高導異常響應與無明顯異常的背景響應在5 ms前沒有明顯區別，因為早期響應大多貢獻自儀器附近的海水。在5 ms后，二者開始產生差異，深部高導體的存在導致了更強的二次場，從而使異常響應衰減變慢。將存在異常的響應所在的位置投點至平面圖可以發現，高導異常出現的位置均在TAG丘體上方，沒有異常的背景響應則采集自遠離丘體的測點上。晚期20 ms處異常最為明顯的響應來自于TAG丘體的山頂處，此處礦體厚度最大。

圖7 瞬變電磁重疊回線系統接收線圈響應（a）和高導異常所在測點位置（b）Fig.7 Transient electromagnetics response for coincident loop system (a) and the locations that consist high conductivity anomalies (b)

我們將正演模擬結果與TAG熱液區實測瞬變電磁數據進行對比。經過預處理后的部分實測數據如圖 8 所示，采集有效時窗范圍為 3～30 ms，30 ms后的數據因信噪比過低需被剔除。黑色曲線為已知TAG丘體的瞬變電磁響應，淺紅色與灰色曲線分別為離底高度20 m時異常最大與最小的正演模擬結果。可知TAG丘體的三維正演結果與該區域的實測數據能夠很好地擬合，且實測數據表明，已有效探測到了高電導率異常。

圖8 TAG 熱液區部分實測瞬變電磁數據Fig.8 Part of the collected TEM data at the TAG hydrothermal field

為研究儀器離底高度變化對海底多金屬硫化物探測效果的影響，把離底高度分別設置為40 m、60 m和80 m來對比其二次場響應特征。為方便比較，繪制電壓平剖圖如圖9所示。綜合TAG丘體地形可以看出，由于地形起伏及儀器拖曳方式改變了線圈和高導礦體之間的耦合關系，使早期響應的異常形態發生畸變。而對于晚期響應，當礦體厚度逐漸變大時，響應的幅值也變大。對比 20 m，40 m、60 m 與 80 m的正演模擬結果，可知儀器離TAG丘體越遠，異常幅值越小，且早期畸變響應也變小。理論上而言，即使是離底高度為80 m時，晚期的瞬變電磁響應也能探測到丘體的高導異常，但實際上，實測數據是存在噪聲干擾的。以中國大洋26航次在TAG熱液區采集到的瞬變電磁數據為例，實測數據的有效時窗范圍為 1～30 ms，30 ms后的數據由于幅值很小，受到的噪聲干擾也就隨之變大，往往被剔除。圖中的方框虛線即代表30 ms時的背景場平剖曲線，即有效信號的范圍在該曲線以上，此時可以發現，離底高度80 m的平剖曲線已幾乎與背景場重合，無法有效反映出TAG丘體的高導異常。

圖9 不同離底高度下的電壓平剖圖（a）和對應測線軌跡（b）Fig.9 Section of TEM response under different survey altitudes (a) and the corresponding survey lines (b)

在實際近底勘探應用中，瞬變電磁重疊回線儀器會因拖曳過程及洋流擾動而晃動。為研究儀器姿態對瞬變電磁響應的影響，在線圈模型建立時引入傾角，以X軸方向為例，當線圈位于TAG丘體正上方20 m及80 m時，分別計算線圈傾角為45°和60°時的接收線圈響應，結果如圖10所示。可以看出由于激發的一次場隨著線圈發生傾斜，導致二次場產生變化，同時接收線圈的傾斜也使得實測數據產生干擾。對于離底高度為20 m的高導礦體響應，這一干擾主要發生在信號的早期（小于 3 ms）及晚期（大于 70 ms）；而對于離底高度為80 m的響應，姿態變化產生的干擾減弱，且集中在30 ms之后。只有當儀器離底高度遠大于儀器探測深度時，所測得的響應可等效為將儀器置于無限空間海水中的響應，此時線圈的姿態變化則不會對所測響應產生影響。因此，線圈的離底高度和海底地質體電導率分布都會使其在姿態不穩時的響應產生干擾。

圖10 線圈在不同離底高度下，沿X軸傾角變化時的接收線圈響應曲線Fig.10 Transient electromagnetic method response under different loop angle along the X-axis and different altitude

在上文中，我們發現正演模擬的早期響應存在明顯的畸變，這一現象在實測數據中也有體現。為更好地了解這一現象產生的原因，有必要進一步探究地形起伏以及測線布置對瞬變電磁響應的影響。設置3組不同的正演模型：起伏地形+水平測線、水平地形+水平測線以及水平地形+起伏測線。其中，水平地形指的是將原TAG丘體區域地形轉換為水平地形，地下高導體形態與之前模型一致；起伏地形為原TAG丘體區域；水平測線的離底高度為20 m，起伏測線為原模型測線。分別進行正演模擬，得到不同的正演模擬結果（圖11）。

圖11 不同地形及測線條件下的瞬變電磁響應平剖圖Fig.11 Section of TEM response under different topography and survey line conditions

起伏測線條件下，當探測過程不受熱液區崎嶇地形影響時，可以發現其與起伏地形下的平剖圖相比，早期響應更為平緩；此時的畸變信號來自儀器離底高度的變化，因為其改變了重疊回線裝置以下的電導率分布。對起伏地形進行水平測線觀測，可以發現早期響應不再發生畸變，其幅值變化可以更好地反映海底良導體的厚度，但固定深度的觀測方式會使丘體兩側水深更深的礦體響應發生衰減，需要結合晚期響應才可圈定異常范圍。

6 結論

本文采用有限元法實現了海洋瞬變電磁響應的三維正演模擬。通過多波束、鉆探和地質資料建立了復雜地形與結構的TAG丘體三維電導率模型，并計算得到了電磁場的時空分布，正演響應結果與TAG熱液區的實測數據能夠很好地擬合。研究了儀器不同離底高度與姿態條件下的瞬變電磁響應規律，以及地形因素和測線設置對瞬變電磁響應的影響。總結得到以下認識：（1）電流關斷后的磁場在富含多金屬硫化物的TAG丘體內產生渦流，其接收線圈響應特征與高阻體響應特征具有明顯差異，因此瞬變電磁重疊回線裝置可以有效探測海底多金屬硫化物堆積體；（2）瞬變電磁的有效探測深度與所測數據的時窗有關，以TAG熱液區實測瞬變電磁數據的信噪比為例，應保證瞬變電磁重疊回線裝置在拖曳過程中離底高度不高于60 m，才能有效探測到高導體異常；（3）復雜的熱液區地形以及儀器拖曳過程中高度的變化會嚴重影響瞬變電磁早期響應，因此在后續數據處理過程中應結合精確的海底地形數據以及超短基線定位數據，以防止對實測數據的錯誤解釋；或者在保證離底高度不大于探測深度的前提下，采用水平測線的拖曳方式以降低早期響應的擾動，簡化后續處理反演流程；（4）拖曳過程中線圈姿態變化會干擾瞬變電磁的響應，所以有必要記錄探測裝置的姿態，來更好地分析所探測異常的來源。本文驗證了瞬變電磁法在海底硫化物勘探中良好的應用前景，同時考慮了復雜地形、測線軌跡、儀器姿態因素對數據的擾動。在實際反演過程中，若已知并引入這些參數，將能夠為復雜地形與結構條件下的熱液區硫化物堆積體的精確三維反演提供支撐。