電性源短偏移距瞬變電磁場擴散規(guī)律三維數(shù)值模擬

常江浩，薛國強

(1.中國科學院地質(zhì)與地球物理研究所 中國科學院礦產(chǎn)資源研究重點實驗室，北京 100029；2.中國科學院大學 地球與行星科學學院，北京 100049；3.中國科學院地球科學研究院，北京 100029；4.河北地質(zhì)大學 河北省戰(zhàn)略性關(guān)鍵礦產(chǎn)資源重點實驗室，河北 石家莊 050031；5.長安大學 地球物理場多參數(shù)綜合模擬實驗室(中國地球物理學會重點實驗室),陜西 西安 710054)

0 引 言

隨著經(jīng)濟的快速發(fā)展，中國對礦產(chǎn)資源的需求量逐漸增大。中國成礦地質(zhì)條件良好，資源潛力巨大，但目前礦產(chǎn)勘探開采深度大多小于500 m，受到勘查技術(shù)和裝備方面的限制，深部礦產(chǎn)資源無法被有效識別。因此，加強對深部資源的探測，提高危機礦山的深部和外圍找礦能力，發(fā)現(xiàn)接替資源，是勘探工作的重中之重[1-2]。瞬變電磁法是一種建立在電磁感應原理基礎(chǔ)上的時間域電磁探測方法，具有勘探深度大、不受高阻層屏蔽、工作效率高等優(yōu)點，在深部礦產(chǎn)資源勘探中發(fā)揮著重要作用[3-7]。

瞬變電磁法按照發(fā)射源形式可分為回線源瞬變電磁法和電性源瞬變電磁法。回線源裝置不受接地條件限制，對低阻體分辨能力高，在金屬礦勘探、煤礦水害探查和工程勘察等領(lǐng)域有著廣泛的應用[8-12]。但回線源對高阻體分辨能力較差，其場能量在地層中衰減較快，一般只能達到幾百米的探測深度。電性源采用接地導線源發(fā)射電磁場，其產(chǎn)生的電場不僅具有水平分量，還具有垂直分量，電場水平分量有利于低阻體的探測，垂直分量有利于高阻體的探測[13]。傳統(tǒng)的電性源裝置形式為長偏移距裝置，長偏移距接收信號強度隨收發(fā)距增大將急劇減小，信噪比下降，不利于精細探測。薛國強等針對近源探測的優(yōu)越性，對電性源短偏移距瞬變電磁(SOTEM)響應進行了研究，并提出了電性源短偏移距瞬變電磁法[14-15]。理論和實踐表明，電性源短偏移距瞬變電磁法采用的近源探測形式不僅有較高的信噪比，還能達到較大的探測深度[16-18]。

崔江偉對電性源短偏移距瞬變電磁全程視電阻率求取方法進行了研究，并分析了全程視電阻率對薄層的分辨能力[19]。陳衛(wèi)營等對電性源短偏移距瞬變電磁各分量響應衰減特征進行了分析，并基于一維模型對電場和磁場的敏感區(qū)域進行了研究[20]；之后，又建立了電性源短偏移距瞬變電磁數(shù)據(jù)的一維OCCAM反演方法，并將其應用于三維電性源短偏移距瞬變電磁數(shù)據(jù)的反演[21]。陳康對電性源短偏移距瞬變電磁多分量響應特征進行了研究，并對各個分量的探測能力進行了分析[22]。目前對電性源短偏移距瞬變電磁響應特征的研究主要基于一維模型，且主要討論電場Ex分量和磁場Hz分量的響應特征和探測能力，不利于復雜目標體的精細探查。

本文基于三維時域有限差分法，研究了電性源短偏移距瞬變電磁場在層狀地層中的擴散規(guī)律，并采用三維異常體模型研究了三維異常體對電性源短偏移距瞬變電磁場的影響，以期推進電性源短偏移距瞬變電磁法的基礎(chǔ)理論研究。

1 原 理

1.1 差分方程

時域有限差分法是瞬變電磁場三維數(shù)值模擬的主要方法，其采用的網(wǎng)格單元如圖1所示[23]，這樣的網(wǎng)格形式使三維體表面上場分量的連續(xù)性條件自然得到了滿足。在無源區(qū)域，麥克斯韋(Maxwell)方程組的微分形式[23]為

i、j、k分別為x、y、z方向網(wǎng)格剖分的節(jié)點；Ex、Ey、Ez分別為x、y、z方向上的電場分量；Hx、Hy、Hz分別為x、y、z方向上的磁場分量；圖件引自文獻[23]

(1)

(2)

(3)

(4)

式中：E為電場強度；B為磁感應強度；H為磁場強度；σ為介質(zhì)電導率；ε為介電常數(shù)；t為時間；D為位移電流。

受穩(wěn)定性條件限制，空氣中電磁場的有限差分迭代需要采用較小的時間步長影響計算效率。根據(jù)時間域電磁場的傳播特點，在一次場關(guān)斷后的初期，地面波以光速從空氣中傳播到地表各點，此時的場主要以地面波為主，隨著時間的推移，以地面波傳播的場逐漸衰減，至晚期以地面波傳播的場可忽略不計，晚期的場主要來自地層波。本文根據(jù)場在早期和晚期的特征分別采用了不同的有限差分算法，早期直接采用動態(tài)電磁場方程計算以保證地面波的精度，而晚期采用準靜態(tài)電磁場方程計算以提高計算效率。

在準靜態(tài)近似下，式(2)可忽略位移電流項。為了滿足顯式時間步進格式，Wang等引入了虛擬位移電流項[24]，其表達式為

(5)

虛擬介電常數(shù)的引入使計算晚期場時可以放寬對時間步長的要求[24]。將式(5)進行差分離散即可得到無源區(qū)域準靜態(tài)近似下電場的差分方程，具體離散過程可參考文獻[24]。

1.2 接地導線源的加載及邊界條件

瞬變電磁場直接時域數(shù)值模擬中，發(fā)射源的加載方式主要有3種：①利用負階躍電流在均勻半空間產(chǎn)生的解析解作為場的初始值；②將源電流密度加入麥克斯韋方程，一般采用梯形電流波形；③取關(guān)斷前建立的穩(wěn)定場作為場的初始值。第①種方法要求發(fā)射源附近為均勻介質(zhì)，一般用來計算磁偶極子源產(chǎn)生的初始場；第②種方法對模型沒有限制，可以模擬任意波形發(fā)射電流的瞬變電磁響應；第③種方法需要求解穩(wěn)定電流場，為了避免求解磁場初始值，在迭代中一般用磁場的時間導數(shù)代替磁場。本文采用第②種方法加載發(fā)射源，該方法對模型沒有限制，可以模擬任意發(fā)射電流激發(fā)的瞬變電磁響應。

在有源區(qū)域，式(2)可修改[25]為

(6)

式中：Js為源電流密度。

在施加接地導線源時，可將導線源與電場分量的空間位置重合。將式(6)進行差分離散即可得到有源區(qū)域準靜態(tài)近似下電場的差分方程，具體離散過程可參考文獻[25]。

邊界條件是影響晚期場精度的主要因素，有效吸收邊界不僅能提高解的精度，而且能減少網(wǎng)格數(shù)量。本文在截斷邊界采用卷積完全匹配層(CPML)吸收邊界條件[26-27]。卷積完全匹配層吸收邊界不僅具有完全匹配層(PML)的性能，還能對低頻場進行有效吸收[28]。

1.3 算法驗證

采用均勻半空間模型對算法的精度進行驗證，地下介質(zhì)電阻率為100 Ω·m，空氣電阻率為1.0×105Ω·m，在地表布設長度為100 m的電性源，在偏移距800 m處接收。發(fā)射電流波形為梯形波，梯形波寬度為100 ms，上升沿和下降沿寬度為10 μs。驗證結(jié)果如圖2所示，將時域有限差分解與解析解進行對比。從電場衰減曲線[圖2(a)]可以看出，時域有限差分解與解析解在各個時段都非常吻合；從相對誤差曲線[圖2(b)]可以看出，最大相對誤差小于5%。這說明本文采用的有限差分法具有較高的精度。

圖2 時域有限差分解與解析解對比

2 電性源短偏移距瞬變電磁場擴散規(guī)律

電性源短偏移距瞬變電磁法采用近源觀測方式，近源處瞬變電磁場的擴散規(guī)律與大偏移距處不同，研究電性源短偏移距瞬變電磁場的擴散規(guī)律對于認識近源處瞬變電磁場與地質(zhì)體的相互作用過程以及最佳觀測場量的選擇具有重要作用。本文采用如圖3所示的層狀地層模型，地下地層被分為3層。第一層和第三層電阻率為100 Ω·m，第二層電阻率為10 Ω·m(低阻層)或1 000 Ω·m(高阻層)，分別代表H型地層模型和K型地層模型。其發(fā)射源長度為200 m，沿x軸方向布設，中點位于坐標原點，發(fā)射電流為1 A。

h1、h2、h3分別為第一、二、三層地層的厚度

采用三維時域有限差分法進行正演模擬計算，圖4～6分別為電場各分量在y=200 m剖面上的等值線分布。從圖4可以看出，電場Ex分量下部出現(xiàn)負值區(qū)域，為返回電流影響。在1 ms時，電場Ex分量上部正極值區(qū)域主要集中在發(fā)射源附近，隨著時間的推移逐漸向下擴散；在10 ms時，H型地層模型的電場Ex分量極值區(qū)域位于低阻層附近，而K型地層模型的電場Ex分量極值區(qū)域已經(jīng)穿過高阻層，說明電場在低阻層中擴散速度較慢，在高阻層中擴散速度較快。從圖5可以看出，電場Ey分量正極值區(qū)域和負極值區(qū)域分別位于發(fā)射源兩端，并隨著時間的推移逐漸向下、向外擴散。從圖6可以看出，電場Ez分量在地層分界面急劇變化，這主要受地層界面處存在的面電荷影響。對于H型地層模型，低阻層內(nèi)部電場Ez分量幅值小于上方和下方地層，而K型地層模型的高阻層內(nèi)部電場Ez分量幅值大于上方和下方地層。以上電場的三維模擬結(jié)果與陳衛(wèi)營等的一維模擬結(jié)果[29]基本吻合，進一步說明了本文采用的三維模擬方法的可靠性。

圖4 不同時刻電場Ex分量等值線分布

圖5 不同時刻電場Ey分量等值線分布

圖6 不同時刻電場Ez分量等值線分布

圖7～9為磁場各分量在y=200 m剖面上的等值線分布。從圖7可以看出，磁場Hx分量正極值區(qū)域和負極值區(qū)域分別位于發(fā)射源兩端，并隨著時間的推移逐漸向下、向外擴散。從圖8可以看出，磁場Hy分量上部為負值區(qū)域，下部為正值區(qū)域，說明磁場Hy分量也會受到返回電流影響。比較H型和K型地層模型的結(jié)果，磁場Hy分量在低阻層中擴散較慢，在高阻層中擴散較快，且低阻層對磁場Hy分量擴散影響較大。圖9為磁場Hz分量的等值線分布，磁場Hz分量極值區(qū)域主要位于發(fā)射源下方，并隨著時間的推移逐漸向下移動。受中間異常層的影響，在10 ms時，H型地層模型的磁場Hz分量極值區(qū)域位于低阻層附近，而K型地層模型的磁場Hz分量極值區(qū)域已經(jīng)穿過高阻層，說明低阻層擴散速度較慢，高阻層擴散速度較快。整體上，低阻層對磁場Hz分量擴散影響較大，高阻層對磁場Hz分量擴散影響較小。

圖7 不同時刻磁場Hx分量等值線分布

圖8 不同時刻磁場Hy分量等值線分布

圖9 不同時刻磁場Hz分量等值線分布

3 三維異常體影響下電性源短偏移距瞬變電磁場平面分布特征

為認識三維異常體影響下電性源短偏移距瞬變電磁場的擴散規(guī)律，建立了如圖10所示的含三維異常體的地電模型。均勻半空間含有一低阻異常體，大小為200 m×200 m×200 m，埋深為100 m。發(fā)射源長度為200 m，沿x軸方向布設，中點位于坐標原點，發(fā)射電流為1 A。異常體位于發(fā)射源赤道位置，距離發(fā)射源400 m。均勻大地電阻率設為100 Ω·m，異常體電阻率為10 Ω·m。

圖(b)括號中的數(shù)據(jù)為電阻率

采用三維時域有限差分法進行正演模擬計算，圖11、12為發(fā)射電流關(guān)斷后不同時刻在地表觀測到的電場Ex、Ey分量的等值線分布(在地表無法觀測電場Ez分量)。從圖11可以看出：在1 ms時，電場Ex分量在異常體位置發(fā)射畸變，等值線凹向發(fā)射源方向，說明低阻異常體使電場擴散速度減慢；在10 ms時，Ex分量在異常體上方出現(xiàn)極小值，在異常體兩側(cè)出現(xiàn)極大值。從圖12可以看出，電場Ey分量顯示出了4個異常區(qū)域，均位于異常體外部且分布在異常體四周。

圖11 異常體影響下電場Ex分量平面分布

圖12 異常體影響下電場Ey分量平面分布

圖13～15為發(fā)射電流關(guān)斷后不同時刻在地表觀測到的磁場各分量的等值線分布。從圖13可以看出，由于受到異常體的吸引，磁場Hx分量在異常體外部產(chǎn)生畸變。從圖14可以看出：在1 ms時，磁場Hy分量的分布受異常體影響較小，隨著時間的推移，其分布受異常體影響逐漸增大；在10 ms時，其等值線出現(xiàn)兩個極小值區(qū)域，一個位于發(fā)射源附近，另一個位于異常體上方，說明低阻異常體會使磁場Hy分量產(chǎn)生負異常。圖15為磁場Hz分量的等值線分布，其等值線在異常體兩側(cè)(靠近發(fā)射源和遠離發(fā)射源位置)畸變較大，而在異常體上方?jīng)]有明顯異常。

圖13 異常體影響下磁場Hx分量平面分布

圖14 異常體影響下磁場Hy分量平面分布

圖15 異常體影響下磁場Hz分量平面分布

以上分析表明，電性源短偏移距瞬變電磁場的5個分量對異常體的靈敏區(qū)域不同。電場Ex分量和磁場Hy分量在異常體上方異常最大；而電場Ey分量和磁場Hx分量異常響應最大值位于異常體外部，且分布在異常體四周；磁場Hz分量在異常體兩側(cè)(靠近發(fā)射源和遠離發(fā)射源位置)異常最大。這說明電性源短偏移距瞬變電磁場的5個分量在地表應選擇不同的區(qū)域進行觀測。

4 結(jié) 語

(1)電場Ex分量的正極值區(qū)域主要集中在發(fā)射源附近，并隨著時間的推移逐漸向下擴散；電場Ey分量正極值區(qū)域和負極值區(qū)域分別位于發(fā)射源兩端，并隨著時間的推移逐漸向下、向外擴散；電場Ez分量在地層分界面產(chǎn)生躍變，說明在地層分界面處存在面電荷。

(2)磁場Hx分量正極值區(qū)域和負極值區(qū)域分別位于發(fā)射源兩端，并隨著時間的推移逐漸向下、向外擴散；磁場Hy分量上部為負值區(qū)域，下部為正值區(qū)域，說明磁場Hy分量也會受到返回電流影響；磁場Hz分量極值區(qū)域主要位于發(fā)射源下方，并隨著時間的推移逐漸向下移動。

(3)對于三維異常體，電場Ex分量和磁場Hy分量的靈敏區(qū)域位于異常體上方；而電場Ey分量和磁場Hx分量的靈敏區(qū)域位于異常體外部，且分布在異常體四周；磁場Hz分量的靈敏區(qū)域位于異常體兩側(cè)(靠近發(fā)射源和遠離發(fā)射源位置)。這說明電性源短偏移距瞬變電磁場的5個分量對三維異常體的靈敏區(qū)域不同，不同的分量應選擇不同的區(qū)域進行觀測。