全拖曳式深海直流電阻率法三維任意各向異性正演模擬

楊志龍，殷長春，張 博，劉云鶴，任秀艷，惠哲劍

吉林大學地球探測科學與技術(shù)學院，長春 130026

0 引言

近年來，深海沉積礦產(chǎn)資源，尤其是以高品質(zhì)錳結(jié)核為代表的多金屬結(jié)核礦物的發(fā)現(xiàn)及深海可燃冰提取技術(shù)的進步，加快了人們對深海礦產(chǎn)資源的勘探需求。深海埋藏管線的檢測和海底污染物的運移等環(huán)境工程問題也越來越受到重視。國際深海勘探領(lǐng)域“圈地運動”蓬勃興起。由于深海沉積層的電阻率包含了其物理性質(zhì)的豐富信息，因此開展海洋直流電阻率法研究在海洋礦產(chǎn)探測、海底管線檢測和海底污染物運移等方面具有很大的實用價值。海底是非常穩(wěn)定的環(huán)境，上方海水可近似看作電各向同性均勻介質(zhì)，電極周圍環(huán)境干擾小并且與海水之間的阻抗小，易于獲得很大的供電電流，這在地面直流電阻率法中是很難做到的。除此之外，直流電場空間分布受電阻率的變化影響較大，對海底的高阻油氣層和低阻金屬礦產(chǎn)識別能力較強。由于目前三維海洋直流電阻率法發(fā)表的成果較少，因此開展高精度三維海洋直流電阻率正演模擬研究具有重要理論價值。

制約海洋直流電法發(fā)展的一個重要因素是探測方式的效率。20世紀60年代，海洋電阻率法開始用于調(diào)查北冰洋大陸架上的永凍層和礦產(chǎn)資源勘探，但當時的接收裝置對巖性的區(qū)分效果較弱[1]；Hakim等[2]在對海底進行填圖時提出了兩種特殊的布極方式，提高了電阻率法對巖性的識別能力。Rosenberger等[3]設計了一種自由落體式的電阻率海底沉積物原位探測裝置，取得了很好的應用效果。Edwards等于1986年引入了深海海底勘探瞬變電磁方法理論，于1987年對該理論進行了完善，并于1988年對二維模型進行了模擬計算的展示[4-6]。但以上的探測方式大多適用于淺海直流電阻率法勘探；在深海探測中，由于海底洋流活動異常活躍且活動規(guī)律未知，傳統(tǒng)探測方式中放置海底的接收電極往往會造成丟失和測量偏差等問題，嚴重影響數(shù)據(jù)的質(zhì)量。考慮到以上種種原因，基于前人經(jīng)驗，本文引入一種同心掃面深海全拖曳探測方式。即一條船拖著發(fā)射電極靜止在海平面上，以該船為中心，通過手持GPS儀器提前在海平面布設一系列回形測線，計算并存儲測線上每個測點的位置；然后另一條探測船拖曳接收電極，根據(jù)GPS儀器導航，沿測線依次記錄各測點所接收的電位信息；最后，計算得到拖曳平面內(nèi)測點的視電阻率信息，通過畫圖軟件展示在平面圖上，從而可以形象地觀測拖曳平面各個方位的視電阻率。這種探測方式不僅避免了傳統(tǒng)探測方式接收電極丟失的風險和測量偏差的難題，還實現(xiàn)了拖曳平面內(nèi)的多方位探測，通過計算拖曳的平面內(nèi)各個方位的視電阻率可有效識別海底沉積層電性分布特征，特別是電各向異性信息。

空氣層基本不導電，因此在地面直流電阻率法中，地面點源電流向地下呈半球面擴散。而在海洋直流電阻率法中，海水中電流源的發(fā)射電流呈球面擴散，且發(fā)射電流在空氣與海水的分界面處會發(fā)生反射，因此海洋直流電法的正演理論與地面直流電法不同，控制方程中需要使用鏡像法。葛為中[7-8]對一維層狀介質(zhì)的點源場進行了系統(tǒng)研究，基本上解決了水中一維直流電場的正演模擬問題；并于1998年推導出了一維層狀介質(zhì)水下點電源的電位和視電阻率理論公式，討論了水底電測深轉(zhuǎn)換函數(shù)之間的關(guān)系，同時指出水下或水底電阻率法由于電極裝置接近沉積層，比水面電阻率法具有更好的探測能力。王自民[9]論述了在水面上進行電阻率測量的方法。黃俊革等[10]利用規(guī)則網(wǎng)格有限單元法對水下電阻率進行模擬，指出水下電極深度和水底地形都會對測量結(jié)果產(chǎn)生影響，但由于規(guī)則網(wǎng)格對地形模擬較差和所選第一類邊界條件的局限性，取得的模擬結(jié)果精度并不理想。翁愛華等[11]對一維海洋電阻率法進行模擬，詳細討論了深海探測中的影響因素，并指出海洋電阻率法對分辨油氣和礦產(chǎn)資源的優(yōu)越性。在電磁法建模方面：劉長勝等[12]計算了海水為均勻半空間的瞬變電磁響應；李慧[13]計算了淺海瞬變電磁垂直偶極子和中心回線裝置的響應曲線；周勝等[14]計算了深海拖曳式瞬變電磁的響應規(guī)律。

正演是反演的基礎(chǔ)，由于介質(zhì)各向異性會帶來更多的計算參數(shù)，因此目前的正反演只考慮了介質(zhì)的各向同性或簡化的電導率各向異性。而深海沉積物是典型的電各向異性分布[15]，因此在當前的正反演中往往隱含著一些不可忽略的錯誤[16]。在三維各向異性正演計算中可以通過加入三維電導率張量來模擬介質(zhì)在空間中的任意電導率各向異性，從而實現(xiàn)仿真模擬與實際地質(zhì)模型的完全匹配，保證正演模擬的可信度，提高計算的精確度。對于一些簡單的各向異性模型，目前已經(jīng)得到了均勻半空間傾斜各向異性模型的解析解[17]以及各向異性層狀介質(zhì)的地電響應計算結(jié)果[18-21]，但三維正演模擬結(jié)果發(fā)表較少。本文首先提出針對海底礦產(chǎn)勘查的全拖曳式直流電阻率勘查技術(shù)。進而從泊松方程出發(fā)，以鏡像法為基礎(chǔ)，結(jié)合精度較高的第三類邊界條件，通過引入電導率張量求解海洋各向異性環(huán)境下直流電法三維正演問題。為了解決源的奇異性問題，本文采取了一次場與二次場分離的方法，一次場從半空間的解析解獲得，而二次場采用有限元法求取。對于海底起伏物性界面，本文選用了可以模擬任意形狀的非結(jié)構(gòu)網(wǎng)格進行模擬。本文算法與一維半解析解的對比驗證了算法的精確性。最后，通過對海底任意各向異性模型的數(shù)值模擬證明本文算法的有效性，并驗證海洋直流電法對海底各向異性目標體的探測能力。

1 正演理論

1.1 電導率張量

對于任意各向異性地下介質(zhì)，電導率張量σ(σ=ρ-1，ρ為電阻率)可用一個含有6個獨立分量的正定對稱張量表示，即

(1)

σ和主軸電導率張量σc之間存在如下關(guān)系：

σ=RTσcR。

(2)

式中，R=RxRyRz為歐拉旋轉(zhuǎn)矩陣，具體形式由Yin[22]給出。

1.2 基于二次場的邊值問題

假設區(qū)域Ω中存在一個點電源，由其激發(fā)產(chǎn)生的電位u=us+up滿足Poisson方程，up和us分別表示一次場電位和由異常體激發(fā)的二次電位，海平面處的電位滿足Newman邊界條件，區(qū)域其他邊界Γ處的電位滿足第三類邊界條件，則二次場滿足的邊值條件表示如下[23]：

·σ·us=

(3)

(4)

σ·us·n=

-qus-σ·up·n-qPup，rΓ。

(5)

其中：

(6)

(7)

而

B=rT·σ-1·r；B=r·σ-1·r；

式中：n為邊界法向量；σ0是圍巖電導率張量；r是網(wǎng)格節(jié)點到源的距離；r是網(wǎng)格節(jié)點到鏡像法中源對應像位置的距離。

電流強度為I的一次電位的解析解可表示為[19]

(8)

在Hilbert函數(shù)空間H1(Ω)內(nèi)取測試函數(shù)N，利用Galerkin有限單元法[24]可得

(9)

我們使用開源代碼TetGen[25]將區(qū)域Ω剖分成Delaunay四面體網(wǎng)格，利用格林恒等式化簡式(9)后代入邊界條件，并對二次場在各個單元中進行插值，得到如下線性方程組：

KUs=b。

(10)

式中：K是全局系數(shù)矩陣；b是右端源項；Us是待求的二次電位列向量。使用并行直接求解器MUMPS求解大型稀疏線性方程組可得到二次場電位，最后與一次場電位相加即可得到網(wǎng)格節(jié)點的總電位。在獲得測點的總電位后，考慮到點源在介質(zhì)中，所以利用鏡像法計算公式來計算相應的視電阻率，計算公式為

(11)

式中：d表示真實點源到測點的距離；d′表示點源的鏡像到測點的距離；u表示測點的總電位。

1.3 非結(jié)構(gòu)網(wǎng)格局部加密技術(shù)

直接由Delaunay refinement triangulation生成的非結(jié)構(gòu)化網(wǎng)格可以通過引入加密技術(shù)來提高數(shù)值解的精度[25-26]。目前的加密技術(shù)大體分為全局加密技術(shù)和局部加密技術(shù)。局部加密技術(shù)相對全局加密來說會大大減少計算量，本文采用了局部加密技術(shù)中的局部節(jié)點加密[27]和局部面積加密。為了提高點源附近的計算精度，我們選擇在源附近均勻地增加節(jié)點以對源附近的網(wǎng)格進行加密；二次場由于異常體和圍巖接觸面處的電阻率差異而產(chǎn)生，因此我們在網(wǎng)格生成器中選擇限制異常體表面的三角形單元大小，以加密網(wǎng)格、提高二次場的計算精度。

1.4 精度驗證

本文算例運行平臺為：Intel(R) Xeon(R) CPU E5-2667 v3 @ 3.20 GHz 3.20 GHz。我們采用任意各向異性的兩層模型(圖1)來驗證本文算法精度。第一層海水的電阻率ρ0=0.3 Ω·m，第二層海底沉積物的主軸電阻率及3個歐拉角為：ρx=1 Ω·m，ρy=1 Ω·m，ρz=4 Ω·m，α=0°，β=0°，γ=0°。考慮到現(xiàn)有的層狀介質(zhì)模型結(jié)果均采用單極-單極裝置，首先假設兩個發(fā)射點源A、B。坐標分別位于(-0.5，0，0)和(0.5，0，0)處，觀測各接收點的電位差并計算視電阻率。圖2給出了本文算法和Yin等[28]半解析解的對比結(jié)果。

圖1 兩層海洋各向異性模型Fig.1 Two-layered anisotropic model for sea

由圖2a可知：當發(fā)收距較小時，海水電阻率得到很好的反映；隨著發(fā)收距的增大，海底沉積層的電阻率得到越來越明顯的反映。由于受海水的影響較大，圖中收發(fā)距較大時得到的電阻率是海水和海底沉積層的綜合響應。另外，由圖2可以看出，本文算法的結(jié)果與半解析解吻合很好，最大相對誤差為2%，說明本文算法具有較高的精度。

2 算例

2.1 模型建立

為研究海底金屬礦床及圍巖各向異性的電阻率響應特征，我們參考公開發(fā)布的海底礦物電學性質(zhì)和礦藏規(guī)模[29]，主要依據(jù)已有的正演模擬文獻選取海底沉積層電阻率[28]。由于深海多金屬結(jié)核尤其是錳結(jié)核以錳鈷鎳氧化物為主，因此主要參考軟錳礦的電阻率范圍選取異常體電阻率[14]，建立如圖3所示的模型。在10 km10 km10 km的模型空間中，建立z軸垂直向下的笛卡爾坐標系，坐標原點位于海平面中心[30]。海水深度為3 000 m，電阻率為0.3 Ω·m，海底沉積物的主軸電阻率為ρ0x=1 Ω·m，ρ0y=2 Ω·m，ρ0z=4 Ω·m，沉積物中金屬礦埋深30 m，長寬均為1 000 m，厚度為800 m，主軸電阻率為ρ1x=0.5 Ω·m，ρ1y=0.002 5 Ω·m，ρ1z=0.01 Ω·m。假設電流發(fā)射船靜止于海平面，下方拖曳著距海底30 m的發(fā)射電極，接收船拖曳著距海底30 m的接收電極，并在回形測線上沿測點依次完成測量。我們采用單極-單極裝置。

圖2 本文算法與半解析解結(jié)果對比Fig.2 Comparison of the apparent resistivity computed by finite element method of this paper with semi-analytical solution for the model in Fig.1

圖3 三維海洋直流電法模型示意圖Fig.3 3D model for marine DC resistivity modeling

2.2 任意各向異性海底沉積層對視電阻率的影響

為了研究海底地層和異常體各向異性對視電阻率響應的影響特征，本文針對圖3給出的三維模型計算了各種電各向異性組合模型的視電阻率，結(jié)果如圖4所示。為對比，圖4同時展示了各向異性圍巖中存在(圖4e、f、g、h)和不存在(圖4a、b、c、d)各向異性異常體時的海底視電阻率響應。

由圖4a可以看出，視電阻率等值線呈橢圓狀分布，在實際電阻率小的方向被壓縮，在實際電阻率大的方向被拉伸。由于受海水影響較大，視電阻率體現(xiàn)海水與海底沉積層的綜合影響。圖4b描述的是當海底沉積層繞x軸旋轉(zhuǎn)45°時的視電阻率分布，此時y軸方向的實際電阻率變大，x軸方向的電阻率不變，因此與圖4a對比，圖4b中沿y方向的視電阻率拉伸更大。圖4c描述的是海底沉積層繞x軸和z軸旋轉(zhuǎn)90°時的視電阻率分布，由于繞x軸旋轉(zhuǎn)90°，圖中視電阻率等值線相對于圖4a、b長軸達到了最大拉伸；同時由于地層繞z軸旋轉(zhuǎn)90°，因此視電阻率也相應地繞z軸順時針旋轉(zhuǎn)了90°。圖4d描述的是當?shù)貙永@x軸和z軸分別旋轉(zhuǎn)90°和135°時的視電阻率分布，與圖4c相比，視電阻率等值線很好地反映了地層各向異性特征，極性圖方向有效地指示了海底沉積地層的電性主軸方向。由此可見，視電阻率分布可準確地反映出海底沉積層的各向異性分布信息。

圖4e描述的是電各向異性海底沉積層中存在各向異性異常體時的視電阻率分布，與圖4a對比可以發(fā)現(xiàn)：受各向異性圍巖的影響，異常體邊緣處視電阻率等值線出現(xiàn)畸變，異常體內(nèi)部的視電阻率等值線較好地反映了各向異性異常體的電性分布情況；但由于受圍巖影響較大，異常體對應位置的視電阻率等值線在x軸方向的拉伸與y軸方向的壓縮并不十分明顯。圖4f給出沉積層主軸電阻率繞x軸旋轉(zhuǎn)45°、異常體主軸電阻率不旋轉(zhuǎn)時的視電阻率分布：與圖4b對比可以發(fā)現(xiàn)，異常體邊緣的等值線出現(xiàn)畸變，由此可大致判斷異常體邊界位置；與圖4e對比可知，異常體內(nèi)部的視電阻率等值線在實際電阻率小的方向被拉伸、實際電阻率大的方向被壓縮，這是受圍巖電阻率影響的結(jié)果，可以根據(jù)這個現(xiàn)象推測圍巖主軸電阻率的旋轉(zhuǎn)方向。圖4g描述的是沉積層主軸電阻率繞x軸和z軸各旋轉(zhuǎn)90°、異常體主軸電阻率繞x軸旋轉(zhuǎn)45°時的視電阻率分布，與圖4c對比發(fā)現(xiàn)：由于異常體主軸電阻率繞x軸旋轉(zhuǎn)，導致橢圓等值線的短軸被壓縮；由于受圍巖的影響，異常體邊界處形成了短軸向外凸的花瓣狀曲線。圖4h描述的是沉積層主軸電阻率繞x軸和z軸分別旋轉(zhuǎn)90°和135°、異常體主軸電阻率繞z軸旋轉(zhuǎn)135°時的視電阻率分布，與圖4d對比發(fā)現(xiàn)：異常體內(nèi)部等值線與外部圍巖的視電阻率橢圓等值線正交，很好地區(qū)分異常體和圍巖各向異性的差異；異常體邊界處受圍巖影響，等值線出現(xiàn)畸變。

綜合上述，深海直流拖曳裝置可以有效地反映海底沉積層和良導礦體的電性分布特征，異常體邊界處的視電阻率畸變等值線可以有效地指示異常體邊界位置，視電阻率等值線的分布可有效識別圍巖和異常體的各向異性特征。

2.3 海底地形對視電阻率的影響

深海探測中，由于海底火山運動活躍，海底地形往往高低起伏，對實際探測造成一定影響。為了更好地模擬實際起伏地形的視電阻率響應，我們建立了如圖5所示的不同海底起伏地形模型。圖5a1、b1、c1、d1中的異常體長寬均為1 000 m，厚800 m。

圖5a1為海底無地形起伏時的正演模型，異常體上表面距海底30 m；圖5b1為海底同時存在海峰、海谷地形時的正演模型，海峰、海谷寬度均為2 000 m，海峰高50 m，海谷深50 m，礦體位于坐標原點正下方3 030 m；圖5c1為海底存在海谷地形時的正演模型，海谷寬度為4 000 m，深50 m，礦體位于坐標原點正下方3 080 m；圖5d1為海底存在海峰地形時的正演模型，海峰寬度為4 000 m，高50 m，礦體位于坐標原點正下方2 980 m。圖5a2、b2、c2、d2則為a1、b1、c1、d1對應的視電阻率響應圖。

圖5a2為海底無地形起伏時各向異性模型的視電阻率響應，主要作為分析含地形起伏時視電阻率變化的參考圖件。對比圖5a2、b2、c2、d2可以發(fā)現(xiàn)，海峰、海谷地形造成對應于異常體的視電阻率異常不對稱、邊界被光滑，整體異常由海谷向海峰一側(cè)移動。

3 結(jié)論

本文基于點源電場的泊松方程，利用非結(jié)構(gòu)節(jié)點有限元方法，并采用一次場與二次場分離技術(shù)，成功實現(xiàn)了全拖曳式深海直流電阻率三維數(shù)值模擬，并通過與半解析解對比驗證了算法的精度。通過對任意各向異性模型的數(shù)值模擬實驗得出如下結(jié)論：

1)全拖曳式深海直流電阻率法可有效勘探海底礦產(chǎn)資源。

2)本文提出的算法可以有效模擬各向異性海洋沉積層環(huán)境中直流電阻率法的視電阻率響應。

3)即使在海底各向異性環(huán)境及海底起伏地形條件下，海底異常體的電性特征(包括各向異性特征)仍然可以得到有效探測和識別。

全拖曳式海洋直流電阻率法探測方式靈活、觀測效率高、勘探成本低，可解決傳統(tǒng)深海探測的系列技術(shù)難題。隨著海洋電阻率勘探儀器的發(fā)展，可望該技術(shù)將在海底礦產(chǎn)資源勘查中發(fā)揮越來越重要的作用。