起伏地形條件下長(zhǎng)偏移距瞬變電磁三維正演

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

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

0 引言

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

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

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

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

1 電磁場(chǎng)時(shí)域矢量有限元

1.1 電磁場(chǎng)矢量有限元控制方程

長(zhǎng)偏移距瞬變電磁法滿足的電場(chǎng)擴(kuò)散方程為

(1)

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

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

(2)

(3)

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

(4)

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

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

(5)

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

(6)

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

(7)

(8)

(9)

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

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

(10)

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

(11)

1.2 初值問(wèn)題及邊界條件

LOTEM一般采用下階躍激勵(lì)。下階躍激勵(lì)的初始電場(chǎng)由以下兩部分組成[15]：

(12)

(13)

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

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

(14)

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

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

(15)

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

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

(16)

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

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

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

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

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

2.2 線性方程組建立及求解

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

(17)

式中：E(k)(t)為k時(shí)刻的電場(chǎng)解向量；Δt為k時(shí)刻前兩道時(shí)間步長(zhǎng)。將式(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]對(duì)式(19)進(jìn)行直接求解。對(duì)于時(shí)域電磁場(chǎng)矢量有限元問(wèn)題，當(dāng)Δt不變時(shí)，方程(19)左端系數(shù)矩陣K固定不變，只需對(duì)其進(jìn)行一次LU分解，將右端項(xiàng)回代求解，可有效提高計(jì)算效率。本文正演模擬均在處理器為AMD Ryzen 9-5950、內(nèi)存為96 GB的工作站運(yùn)行。

3 數(shù)值算例

3.1 三維正演算法驗(yàn)證

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

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

3.2 山脊地形LOTEM三維正演

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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