王新宇, 嚴(yán)良俊, 毛玉蓉, 黃 鑫, 謝興兵, 周 磊

1.油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué))，武漢 430100 2.非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心，武漢 430100

0 引言

隨著地球物理電磁法正反演計算技術(shù)和軟硬件的飛速發(fā)展[1-2]，長偏移距瞬變電磁法(long-offset transient electromagnetic method, LOTEM)因其效率高、勘探深度大、應(yīng)用范圍廣等特點，已成功應(yīng)用于資源環(huán)境調(diào)查、剩余油和儲層壓裂動態(tài)監(jiān)測及地?zé)峥辈榈阮I(lǐng)域[3-8]，應(yīng)用效果良好。但地形地貌的復(fù)雜多樣給LOTEM資料處理解釋帶來較大的困難[9]，若不考慮起伏地形影響，很容易導(dǎo)致LOTEM資料解釋錯誤。因此，研究LOTEM在起伏地形下瞬變電磁場的響應(yīng)特征對實際資料處理解釋具有重要指導(dǎo)意義。

電磁法勘探逐漸由粗放向精細發(fā)展，但電磁法數(shù)據(jù)的合理解釋離不開有效的三維正反演算法。對此，已有大量學(xué)者對瞬變電磁正演算法進行研究[10-21]。Commer等[15]提出了一種并行時域有限差分算法，用于LOTEM電磁場響應(yīng)的計算；B?rner 等[16]提出一種有理Krylov方法，并基于矢量有限元在頻率域求解Maxwell方程，通過頻時轉(zhuǎn)換得到時間域電磁響應(yīng)；Um等[17]基于矢量有限元法，采用后退歐拉格式離散時間步長，并提出自適應(yīng)步長加倍方法，采用直接求解器求解線性方程組，當(dāng)步長不變時，僅需進行一次LU分解，大大加快了瞬變電磁三維正演速度；B?rner等[18]提出使用高階有理Krylov近似，結(jié)合直接求解器有效提高了地面瞬變電磁正演效率；Liu等[19]基于擬態(tài)有限體積法離散時間域電磁場方程，采用有理Krylov子空間方法計算時間域電磁場響應(yīng)，實現(xiàn)地空瞬變電磁快速三維正演，并對比了該方法與隱式后退歐拉算法的計算效率；殷長春等[20]基于法向電流連續(xù)的自適應(yīng)后驗誤差估計有效提高了瞬變電磁計算精度，并研究了海底起伏地形對時間域海洋電磁的影響；Liu等[21]基于有限體積法研究了LOTEM在各向異性地層中的電磁場響應(yīng)特征。

近年來，勘探目標(biāo)區(qū)域起伏地形等因素對電磁數(shù)據(jù)的影響規(guī)律越來越受到國內(nèi)外學(xué)者的重視。Liu等[22]采用邊界元法對航空瞬變電磁不同典型地形下的電磁場響應(yīng)特征進行了數(shù)值模擬；Mitsuhata[23]基于有限元法，采用偽Delta函數(shù)分配偶極子源電流，并對起伏地形下LOTEM和可控源電磁法的電磁場響應(yīng)特征進行對比分析；Zhdanov等[24]提出了非均勻背景電導(dǎo)率下復(fù)雜地電結(jié)構(gòu)中的電磁法三維正演積分方程新算法，該算法在正演中引入地形等非均勻背景模型；強建科等[25]實現(xiàn)了三維起伏地形條件下直流電阻率法的有限單元數(shù)值模擬算法；董浩等[26]研究了起伏地形條件下交錯網(wǎng)格有限差分法的大地電磁測深三維正演；Li等[27]基于矢量有限元法研究了復(fù)雜形狀回線源瞬變電磁響應(yīng)特征，并將其應(yīng)用于復(fù)雜地形下任意形狀回線源瞬變電磁三維正演；Lin等[28]采用交錯網(wǎng)格有限差分法對可控源音頻大地電磁法(CSAMT)在典型地形模型山峰、山谷地形下的電磁場響應(yīng)進行了對比，說明地形對CSAMT電磁場響應(yīng)影響劇烈；Li等[29]研究了復(fù)雜地質(zhì)環(huán)境下電性源瞬變電磁場響應(yīng)特征，但僅僅對磁場分量進行討論，并沒有研究電場分量受地形影響的程度；齊彥福等[30]研究了起伏地形對短偏移距瞬變電磁的影響，通過大量數(shù)值模擬，討論了起伏地形對瞬變電磁響應(yīng)的影響。

鑒于LOTEM在起伏地形下電磁場響應(yīng)規(guī)律研究較少，本文基于矢量有限元法，通過Blender軟件對起伏地質(zhì)模型進行建模，采用開源網(wǎng)格剖分代碼Tetgen[31]離散四面體網(wǎng)格，實現(xiàn)LOTEM三維正演模擬。并通過對層狀模型與典型地電模型數(shù)值模擬結(jié)果的對比分析，驗證算法的正確性和有效性。通過復(fù)雜不同三維模型算例，從理論上探討地形對LOTEM電磁場的影響。

1 電磁場時域矢量有限元

1.1 電磁場矢量有限元控制方程

長偏移距瞬變電磁法滿足的電場擴散方程為

(1)

式中：μ為磁導(dǎo)率；E(r,t)為t時刻r位置的電場矢量；Js(r,t)為t時刻r位置的長導(dǎo)線源電流密度；σ為電導(dǎo)率。

將模擬區(qū)域離散為互不相交的四面體，并引入矢量插值基函數(shù)，將自由度賦予在網(wǎng)格離散后的棱邊上，在任意四面體單元e內(nèi)，電場可展開為

(2)

(3)

當(dāng)方程(1)的E為近似解時，令方程(1)的殘差矢量為

(4)

對于第e個四面體單元，確保其加權(quán)殘差積分為0：

?V(e)N(e)·G(e)dV=0。

(5)

式中：V(e)為第e個單元的體積；G(e)為第e個單元的殘差矢量。利用Galerkin方法，并確保所有四面體加權(quán)殘差總和為0，得到有限元控制方程：

(6)

式中：A為質(zhì)量矩陣；E(t)為棱邊上待求的電場值；B為剛度矩陣；S為外加激勵源項。第e個四面體單元的質(zhì)量矩陣A(e)、剛度矩陣B(e)的第(i,j)個元素和外加激勵源項S(e)的第i個元素可表示為：

(7)

(8)

(9)

式中,σ(e)、μ(e)分別為第e個單元的電導(dǎo)率、磁導(dǎo)率。對于激勵源，將長導(dǎo)線源分割為多段有限長度的電偶極子，將其置于四面體的棱邊上，且電流密度方向平行于棱邊，每個電偶極子的電流密度可表示為[32]

Js(r,t)=δ(r-rs)pI(t)dl。

(10)

式中：δ為脈沖函數(shù)；rs和dl分別為電偶極子的空間位置和長度；p為電流方向；I(t)為t時刻的電流強度。利用脈沖函數(shù)的篩選性質(zhì)，S可表示為

(11)

1.2 初值問題及邊界條件

LOTEM一般采用下階躍激勵。下階躍激勵的初始電場由以下兩部分組成[15]：

(12)

(13)

設(shè)點源電流強度為I，位于rs=(xs,ys,zs)處，利用微分形式的歐姆定律j(r)=σE(r)，根據(jù)電場連續(xù)性得到[33]

·j(r)=Iδ(r-rs)。

(14)

式中，j(r)為位置r處的電流密度。聯(lián)立式(13)和式(14)可得到點源場滿足的三維Poisson方程：

·(σφ(r))=-Iδ(r-rs)。

(15)

對直流電場和時域電磁場采用相同的網(wǎng)格來保證直流電場與矢量電場邊界的一致性，且直流場問題與時域電磁場問題均采用總場方法求解。由于求解區(qū)域較大，在計算過程中對外邊界Γ均施加Dirichlet邊界條件：

φ|Γ=0；(n×E)|Γ=0。

(16)

式中：Γ為計算區(qū)域的外邊界；n為Γ的單位外法向向量。

2 LOTEM三維正演實現(xiàn)

2.1 非結(jié)構(gòu)化網(wǎng)格建模

電磁場有限元前處理是有限元數(shù)值模擬的重要部分，正確且合理的有限元模型以及布局合理的網(wǎng)格密度、質(zhì)量可以有效提高數(shù)值精度與計算速度。本文采用Blender開源軟件建立復(fù)雜地質(zhì)模型。首先將高程數(shù)據(jù)導(dǎo)入Blender軟件進行起伏地形三維建模，生成三維模型文件*.ply。該文件格式簡單，僅僅需要頂點元素列表與面元素列表即可實現(xiàn)空間模型表面網(wǎng)格離散。應(yīng)用開源代碼Tetgen識別模型文件并進行非結(jié)構(gòu)四面體網(wǎng)格離散(圖1)[31]。對于起伏地形的四面體網(wǎng)格剖分，除發(fā)射源、測點與異常體局部加密外，對研究區(qū)域地形也要采用相對較小的網(wǎng)格，以保證數(shù)值模擬精度。

圖1 Tetgen起伏地形網(wǎng)格剖分示意圖

2.2 線性方程組建立及求解

為提高數(shù)值精度，采用二階后退歐拉差分格式對時間進行離散[15]：

(17)

式中：E(k)(t)為k時刻的電場解向量；Δt為k時刻前兩道時間步長。將式(17)代入式(6)可得

(3A+2ΔtB)E(k)(t)=
A[4E(k-1)(t)-E(k-2)(t)]-2ΔtS(k)。

(18)

最終形成大型線性方程組，可表示為

KE=b。

(19)

本文使用直接求解器Pardiso[34]對式(19)進行直接求解。對于時域電磁場矢量有限元問題，當(dāng)Δt不變時，方程(19)左端系數(shù)矩陣K固定不變，只需對其進行一次LU分解，將右端項回代求解，可有效提高計算效率。本文正演模擬均在處理器為AMD Ryzen 9-5950、內(nèi)存為96 GB的工作站運行。

3 數(shù)值算例

3.1 三維正演算法驗證

圖2 層狀地層下垂直接觸帶復(fù)雜三維模型

為驗證層狀模型的計算精度，將上述模型基底電阻率設(shè)為100 Ω·m，計算層狀地層的電磁場響應(yīng)，結(jié)果如圖5所示。由圖5c可知，相對誤差在局部急劇增大，這是由于數(shù)值方法均會產(chǎn)生一定的直流偏移，雖然該偏移對整體精度影響較小，但變號處電磁場數(shù)量級一般都很小，直流偏移對變號處電磁場的影響較為嚴(yán)重，導(dǎo)致變號處誤差較大，這在LOTEM正演模擬過程中難以避免。但dBz/dt與Ey的數(shù)值解與解析解整體擬合良好，誤差基本被控制在3%以內(nèi)。進一步說明了本文算法有效性，可進行LOTEM起伏地形條件下的電磁場數(shù)值響應(yīng)規(guī)律分析。

3.2 山脊地形LOTEM三維正演

為研究起伏地形下LOTEM電磁場的響應(yīng)特征，建立如圖6所示的水平地表和山脊地形模型。山脊地形基底長和寬均為1 000 m，高度為50 m，頂部長寬均為500 m，發(fā)射源長度為200 m，沿y方向布設(shè)，中心點坐標(biāo)為(-2 000 m，0 m，0 m)，發(fā)射電流1 A；地表測點、測線間距均為100 m，共21×21=441個測點；圍巖電阻率為100 Ω·m，在水平地表下存在一個電阻率為10 Ω·m的異常體，大小為400 m×400 m×200 m，頂面埋深200 m；分別計算水平地表與山脊地形模型在測點處的電磁場響應(yīng)。將長導(dǎo)線分割為100段電偶極子；在發(fā)射源和測點處局部加密網(wǎng)格，水平地表最終生成1 456 233個網(wǎng)格， 233 004個節(jié)點，1 692 179條棱邊；山脊地形最終生成1 464 446個網(wǎng)格，234 444個節(jié)點，1 701 832條棱邊；水平地表模型正演占用內(nèi)存11.4 G，求解耗時2 183 s，山脊地形模型正演占用內(nèi)存11.2 G，求解耗時2 214 s。

a. 加密區(qū)域；b. 擴邊區(qū)域。

圖4 不同方法的數(shù)值結(jié)果對比

圖7為水平地表和山脊地形模型在測點處不同時刻電場分量的等值線圖。對于有耗介質(zhì)，電磁波在高阻介質(zhì)中衰減快，在低阻介質(zhì)中衰減慢。但隨著時間的推移，電磁波在有耗介質(zhì)中傳播，電場整體逐漸衰減。圖7中10 μs到10 ms之間，感應(yīng)電流主要集中在發(fā)射源附近，靠近發(fā)射源一側(cè)的電場較大，不滿足平面電磁波擴散。隨著時間的增加，感應(yīng)電流擴散范圍逐漸增大，在100 ms后，地表測點電場分布較為均勻，近似滿足平面電磁波擴散。對于水平地表模型的電場響應(yīng)(a1—f1)，異常體的范圍被很好地反映出來。由于低阻異常體埋深較淺，對于山脊地形模型的電場響應(yīng)(a2—f2)，雖然Ey可以反映出異常區(qū)域的存在，但受地形的影響，在地形周圍會出現(xiàn)假異常，且異常響應(yīng)范圍明顯大于水平地表模型的異常響應(yīng)范圍?？梢?，起伏地形會影響LOTEM的觀測數(shù)據(jù)對地下目標(biāo)體的判別，若不考慮地形因素，勢必會給實際資料解釋帶來困難。

圖5 基底電阻率為100 Ω·m時的數(shù)值解與解析解對比

圖6 水平地表(a)和山脊地形(b)模型

a1—f1. 水平地表；a2—f2. 山脊地形。

3.3 復(fù)雜地形LOTEM三維正演

為研究復(fù)雜地形下LOTEM電磁場的響應(yīng)特征，設(shè)計如圖8所示模型。試驗區(qū)最高海拔為132 m，最低海拔為-68 m，發(fā)射源所在海拔為0 m，發(fā)射源長度為1 000 m，沿y方向布設(shè)，中心點坐標(biāo)為(-1 200 m，0 m，0 m)，發(fā)射電流1 A；y方向的測線位于x=0 m處，測線總長度1 000 m，共41個等間隔測點(圖8a)。淺地表第一層覆蓋層的最大厚度為221 m，最小厚度為41 m，電阻率為100 Ω·m；第二層覆蓋層厚度350 m，電阻率為50 Ω·m；基底圍巖電阻率為100 Ω·m；在第二層覆蓋層和基底中嵌入一個電阻率為10 Ω·m或者1 000 Ω·m的階梯狀異常體，該異常體由3個400 m×100 m×200 m的長方體構(gòu)成，異常體的整體高度為400 m，頂部埋深100 m，y方向正負各延伸200 m；x方向長度為300 m(圖8b)。將長導(dǎo)線分割為400段電偶極子，在發(fā)射源和測點處局部加密網(wǎng)格(圖8c)，最終生成2 777 110個四面體網(wǎng)格，442 984個節(jié)點，3 220 869條棱邊。計算區(qū)域整體尺寸300 km×300 km×200 km，空氣層厚度為50 km；該模型正演占用內(nèi)存32.8 G，求解耗時10 283 s。

圖9為高阻和低阻階梯狀異常體在地表測線的電場響應(yīng)結(jié)果，受地形因素的影響，無論異常體為低阻還是高阻，其異常響應(yīng)基本被“埋沒”。

為進一步研究LOTEM對低阻、高阻異常體的響應(yīng)機理，繪制圖10所示地下電場分布及電流密度矢量圖。由于起伏地形和地下目標(biāo)體的相互耦合作用，使測線附近的電場響應(yīng)發(fā)生了嚴(yán)重的畸變；隨著時間的增加，在有耗介質(zhì)中，電場不斷衰減，可以清晰地看出渦流不斷向下傳播的過程。由圖10a、c、e、g可以看出，由于低阻異常體對電流的吸引作用，感應(yīng)電流在低阻異常體中產(chǎn)生明顯的電流通道效應(yīng)，根據(jù)電流密度法向連續(xù)條件，低阻異常體的電導(dǎo)率較高，感應(yīng)的電場值較小。反之，對于高阻異常體的電場響應(yīng)(圖10b、d、f、h)，電流在高阻異常體周圍產(chǎn)生明顯的電流通道效應(yīng)，高阻異常體的電導(dǎo)率較低，感應(yīng)的電場值較大。

a. 試驗區(qū)等高線；b. 模型yz切面；c. 網(wǎng)格剖分。

a. 低阻異常體；b. 高阻異常體。

a、c、e、g. 低阻異常體；b、d、f、h. 高阻異常體。

4 結(jié)論

本文通過Blender軟件和開源代碼Tetgen實現(xiàn)起伏地形復(fù)雜地質(zhì)模型建模及四面體網(wǎng)格剖分，并基于非結(jié)構(gòu)矢量有限元法實現(xiàn)了長偏移距瞬變電磁(LOTEM)三維正演模擬。通過大量數(shù)值模擬得出以下結(jié)論：

1)通過與解析解、時域有限差分解、有理Krylov方法的擬態(tài)有限體積解的精度對比分析，驗證了本文算法正確性和有效性，微弱的直流偏移雖然對整體計算精度影響較小，但會導(dǎo)致電場變號位置的誤差急劇增大，因此，在數(shù)值模擬或者實際資料處理解釋中，變號附近的場值一般不予考慮。

2)針對起伏地形條件下復(fù)雜地電模型LOTEM的電磁響應(yīng)特征進行分析，表明起伏地形對LOTEM電磁響應(yīng)影響較大，會“削弱”甚至“埋沒”異常體的電磁響應(yīng)，從而無法對異常體的位置和電性特征進行有效分辨。

3)通過對地下電場的等值線及電流密度矢量分布圖研究，由于地形與異常體的相互耦合作用，測線附近的電場畸變嚴(yán)重，使得地表的觀測電場難以對地下目標(biāo)體進行有效判別。

整體研究表明，起伏地形對LOTEM電場響應(yīng)影響較大，在實際勘探中地形效應(yīng)不可忽略。針對以上問題，我們下一步將展開地形校正的研究或開發(fā)反演算法，直接將地形因素考慮到LOTEM的三維正反演中。