起伏地形條件下長偏移距瞬變電磁三維正演

王新宇, 嚴良俊, 毛玉蓉, 黃 鑫, 謝興兵, 周 磊

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

0 引言

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

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

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

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

1 電磁場時域矢量有限元

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

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

(1)

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

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

(2)

(3)

當方程(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個單元的電導率、磁導率。對于激勵源，將長導線源分割為多段有限長度的電偶極子，將其置于四面體的棱邊上，且電流密度方向平行于棱邊，每個電偶極子的電流密度可表示為[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)

設點源電流強度為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開源軟件建立復雜地質(zhì)模型。首先將高程數(shù)據(jù)導入Blender軟件進行起伏地形三維建模，生成三維模型文件*.ply。該文件格式簡單，僅僅需要頂點元素列表與面元素列表即可實現(xiàn)空間模型表面網(wǎ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)進行直接求解。對于時域電磁場矢量有限元問題，當Δt不變時，方程(19)左端系數(shù)矩陣K固定不變，只需對其進行一次LU分解，將右端項回代求解，可有效提高計算效率。本文正演模擬均在處理器為AMD Ryzen 9-5950、內(nèi)存為96 GB的工作站運行。

3 數(shù)值算例

3.1 三維正演算法驗證

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

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

3.2 山脊地形LOTEM三維正演

為研究起伏地形下LOTEM電磁場的響應特征，建立如圖6所示的水平地表和山脊地形模型。山脊地形基底長和寬均為1 000 m，高度為50 m，頂部長寬均為500 m，發(fā)射源長度為200 m，沿y方向布設，中心點坐標為(-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；分別計算水平地表與山脊地形模型在測點處的電磁場響應。將長導線分割為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之間，感應電流主要集中在發(fā)射源附近，靠近發(fā)射源一側(cè)的電場較大，不滿足平面電磁波擴散。隨著時間的增加，感應電流擴散范圍逐漸增大，在100 ms后，地表測點電場分布較為均勻，近似滿足平面電磁波擴散。對于水平地表模型的電場響應(a1—f1)，異常體的范圍被很好地反映出來。由于低阻異常體埋深較淺，對于山脊地形模型的電場響應(a2—f2)，雖然Ey可以反映出異常區(qū)域的存在，但受地形的影響，在地形周圍會出現(xiàn)假異常，且異常響應范圍明顯大于水平地表模型的異常響應范圍?？梢?，起伏地形會影響LOTEM的觀測數(shù)據(jù)對地下目標體的判別，若不考慮地形因素，勢必會給實際資料解釋帶來困難。

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

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

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

3.3 復雜地形LOTEM三維正演

為研究復雜地形下LOTEM電磁場的響應特征，設計如圖8所示模型。試驗區(qū)最高海拔為132 m，最低海拔為-68 m，發(fā)射源所在海拔為0 m，發(fā)射源長度為1 000 m，沿y方向布設，中心點坐標為(-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)。將長導線分割為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為高阻和低阻階梯狀異常體在地表測線的電場響應結(jié)果，受地形因素的影響，無論異常體為低阻還是高阻，其異常響應基本被“埋沒”。

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

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

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

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

4 結(jié)論

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

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

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

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

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