頻率域地－井垂直電磁法響應特征數(shù)值模擬

王智，潘和平，駱玉虎

（中國地質(zhì)大學（武漢）地球物理與空間信息學院，湖北 武漢430074）

0 引 言

井間電磁勘探方法廣泛應用于油氣勘探領域，但金屬套管井限制了井間電磁法橫向探測范圍［1］，目前最大探測距離不超過1km。地－井電磁法能達到4km半徑的勘探范圍，但也只能識別油水平面展布，難以精細分辨儲層縱向變化［2－3］。地－井電磁法是另一種具有潛力的方法，國內(nèi)外對地－井電磁法的研究比較少，Augustin［4］、Wilt等［5］曾對地－井電磁法圈定油氣藏進行了數(shù)值模擬研究；魏寶君等［6］進行了地－井電磁系統(tǒng)成像研究，均取得了一定效果。但是，傳統(tǒng)地－井電磁法橫向分辨率低，受近地表強電磁場耦合干擾嚴重，使該法的應用受到限制。張榮鋒等提出固定井源距垂直電磁剖面地井電磁成像方案，但其分辨率和勘探范圍依然有限［7］。李靜和等［7］提出地－井垂直電磁 Walkaway剖面法（地面激發(fā)、接收裝置固定在井孔中，得到以接收點垂直位置為縱坐標，以不同頻率點或者趨膚泛圍為橫坐標的測深剖面，稱為地－井垂直電磁剖面），主要采用積分方程法對模型進行模擬研究。

電磁法數(shù)值模擬主要有3種方法：有限差分法［8－9］、有限單元法［10－11］和積分方程法［12］。有限差分法是用差商代替微商，將待求解的連續(xù)微分方程變換為離散的差分方程，通過求解差分方程得到原微分方程的近似解；有限單元法是一種求解場的變分問題的數(shù)值方法，通過變分原理，等價變換成變分方程；積分方程法是從場參數(shù)（場強或場位）所滿足的微分方程邊值問題出發(fā)，通過某些變換導出有關參數(shù)（場強、位、積累電荷密度或散射電流密度）所滿足的積分方程，用近似計算方法求積分方程的數(shù)值解。

積分方程法只需對異常體進行剖分和求積，不必計算整個區(qū)域的場，不涉及微分方程法中的吸收邊界等復雜問題。該方法存在對復雜模型的適應性差、方法的精確度嚴重依賴于系數(shù)矩陣的精度以及系數(shù)矩陣的精度往往無法保證等缺點。有限元法推導過程相對簡單，對于變分問題自然邊界條件已經(jīng)隱含地得到滿足，只需考慮強邊界條件，在處理復雜的幾何形狀時適應性好。本文針對頻率域地－井垂直電磁法，應用有限元法對不同的模型進行模擬研究，通過模擬均勻介質(zhì)背景下相對高電阻率與低電阻率條件下目標體的異常特征，分析和討論了模型參數(shù)變化、不同埋深、異常體與收發(fā)位置關系變化時的電磁響應特征和規(guī)律，為了解與研究地－井垂直電磁法的探測能力提供參考。

1 頻率域地－井垂直電磁有限元法基本原理

1．1 工作方法

在井中（沿z方向）從淺至深一定范圍內(nèi)按一定間距布設電場接收裝置，在地面井口處布設激發(fā)場源，激發(fā)頻率（1～1000Hz）從高頻至低頻掃描，以位于不同深度的井中測點為縱軸刻度，以頻率為橫軸，可獲得1張垂直電磁剖面［13］。

圖1 研究區(qū)域

1．2 基本原理

假定地下電性結(jié)構(gòu)是二維，將原點取在地面上，平行于走向方向為x軸，垂直方向為y軸，垂直向下方向為z軸，外加源的方向為x方向，此時只存在不隨走向而變的x方向的電場（見圖1）。取外邊界足夠大，所包圍的區(qū)域為Ω，使局部不均勻體的異常場在外邊界上為0。圖1中，圓圈為局部異常體，Γ1為異常體的邊界，r為電磁波的傳播方向，n為外邊界的外法向方向，θ為r與n的夾角。從麥克斯韋方程組出發(fā)得到的有源電磁場的邊值問題為

邊值問題的式［1（b）］、式［1（c）］為內(nèi)邊界條件，分別表示在介質(zhì)分界面上電場和磁場的切向分量連續(xù)。式［1（d）］為計算區(qū)域截斷邊界即外邊界條件，是由電磁波的傳播方程u＝u0e－kr求導后與傳播方程兩邊分別乘以k后相加得到的。其中，r為源到邊界上任一點的距離，k為空氣中的波數(shù)時，k0＝和η0分別為空氣中的介電常數(shù)和磁導率；當k為地層中的波數(shù)時和μ分別為地層中的介電常數(shù)和磁導率。對于電磁法，頻率足夠低，可忽略位移電流的影響，地層中的波數(shù)可寫為

與式（1）中的邊值問題等價的變分為

1．3 有限元法

區(qū)域Ω采用矩形網(wǎng)格剖分，利用雙線性插值表示小單元內(nèi)的場值（見圖2），圖2中1、2、3、4為小單元角點的編號，Dy、Dz分別為小單元的邊長。

圖2 網(wǎng)格剖分與小單元編號

單元中的場值u可用形函數(shù)表示為

將式（3）代入式（2），經(jīng)整理可得到方程

插值所用的形函數(shù)以及K1e、K2e、P、K4e、K5e系數(shù)矩陣的數(shù)值請參考文獻［14］。

將K1e、K2e、P、K4e、K5e擴展成全體結(jié)點的矩陣并對式（4）取變分

令式（5）為0可得線性方程組KU＝P，U為電場矢量；P為源項。解方程組，即可得到各節(jié)點的電場值。

1．4 線性方程組求解

BICGSTAB算法的主要思想是基于雙邊Lanczos算法和基于殘差正交子空間的迭代算法［15］。預處理矩陣M采用不完全LU分解。一種不完全LU分解是指把A分解為ˉL和ˉU，它們對應的元素與L和U對應的元素相同，但2個因子分別于A的下、上三角部分有完全相同的非0結(jié)構(gòu)。這種分解就是無填充的不完全LU分解，記為ILU（0）［15］。

1．5 壓縮存儲總體系數(shù)矩陣

用有限元方法作數(shù)值模擬時，由小單元集成后所有節(jié)點的總體系數(shù)矩陣是一個大型稀疏矩陣，包含了大量的0元素，網(wǎng)格節(jié)點越多，0值也越多。結(jié)合稀疏矩陣的特點選用較為常用的CSR存貯格式。

2 頻率域地－井垂直電磁法響應特征模擬

2．1 算法正確性的驗證

在上述理論的基礎上，編制了相應的程序。為了驗證程序的正確性和有效性，本文利用均勻半空間的解析解與數(shù)值解對比進行驗證。

對于均勻半空間，供電電源為線源，在地面的觀測場為

式中，K1（z）為第二類一階變形的虛宗量貝塞爾函數(shù)；k為地層中的波數(shù)；k0為空氣中的波數(shù)；y為偏移距；ε和μ分別為地層中的介電常數(shù)和磁導率。圖3中均勻半空間的電阻率值為100．0Ω·m。圖3（a）與圖3（b）分別為f＝4Hz與f＝1024Hz時計算得到的電場值與解析解結(jié)果的對比，從2幅圖中均可以看出數(shù)值解與解析解的結(jié)果十分接近，頻率越高，邊界條件越能減少邊界的影響。這與文獻［13］的結(jié)論一致。

圖3 均勻半空間下的電場數(shù)值模擬結(jié)果與解析解的對比

根據(jù)上述算法模擬均勻半空間條件下的頻率域地井電磁法響應特征。本文模型將源置于中間網(wǎng)格一個節(jié)點上。選取參數(shù)：橫向剖分時源所在位置附近的網(wǎng)格間距較小，向兩側(cè)逐漸加大，網(wǎng)格數(shù)為50；縱向方向網(wǎng)格，空中以較大的網(wǎng)格間距遞增，向地層中先進行均勻剖分，均勻部分網(wǎng)格間距為20m，網(wǎng)格數(shù)為50，然后逐漸增大。所有模擬介質(zhì)中，空氣為不導電介質(zhì)，電導率定為1．0×1012S·m－1。用上述有限元方法和網(wǎng)格剖分方式對勘探區(qū)均勻半空間中所賦存異常體進行模擬。結(jié)合地－井激電法的分析模式，從異常體性質(zhì)、異常體與收發(fā)位置關系分析異常特征規(guī)律。

2．2 異常體性質(zhì)對頻率域地－井垂直電磁法的影響

2．2．1 均勻半空間中低電阻率異常體

埋深200m低電阻率異常體位于縱向11～16號節(jié)點之間，埋深400m低電阻率異常體位于縱向21～26號節(jié)點之間，低電阻率異常體電阻率為10Ω·m，圍巖電阻率為100Ω·m。模型各項參數(shù)見圖4，頻率選擇范圍為20～212Hz，圖5至圖7是利用本文所述方法由低電阻率模型得到的正演模擬結(jié)果。

圖4 低電阻率異常體模型

圖5 頻率1Hz正演模擬結(jié)果

圖6 不同頻率下使用均勻半空間場歸一化后正演模擬結(jié)果

圖7 均勻半空間電場歸一化后的地電斷面圖

從圖5（a）可以看出，場在源附近表現(xiàn)為較大的值，在縱向上，場值隨著網(wǎng)格增大而減小，說明隨著離開源的距離的增大，場在衰減，而在低電阻率異常體的11號網(wǎng)格附近，場值有不太明顯的畸變，說明二次場相對于一次場值很小，異常的信息小，卻依然包含在總場當中，為了對低電阻率體引起的電磁異常有較直觀的了解，將不同頻率產(chǎn)生的電場用均勻半空間的電場做歸一化［14］，以突顯純異常引起的電場變化。正演結(jié)果圖如圖5（b）；從圖6中可以看出，歸一化后異常明顯，隨著頻率的增加，異常極值逐漸向深部轉(zhuǎn)移；圖7為歸一化后的電場地電斷面圖，以縱向網(wǎng)格數(shù)為縱軸，頻率的對數(shù)為橫軸?？梢钥吹矫黠@的異常，異常體附近等值線密集，在異常體周圍呈現(xiàn)較低的電場值，低頻反映異常的能力較弱。對于如此微小的異常能夠在圖7中反映出來，說明本文所述的方法對異常體的縱向分辨能力較強。由圖7可以看出，本文所述的方法對不同的埋深的異常體的探測能力依然有效。

2．2．2 均勻半空間高電阻率異常體

高電阻率異常體位于縱向11～16號節(jié)點之間，電阻率為1000Ω·m，圍巖電阻率為100Ω·m。模型各項參數(shù)如圖8所示。

圖8 埋深200m高電阻率異常體模型

圖9為歸一化后的地電斷面圖。從圖9中可以看到明顯的異常帶，在異常體周圍呈現(xiàn)較高的電場值，在高頻段可以看到明顯的異常，隨著頻率的增加，異常極值依然逐漸向深部轉(zhuǎn)移，異常體附近等值線密集。對比圖7（a）與圖9，低電阻率異常體與高電阻率異常體在本方法中均能較好的反應出來，縱向分辨率較高，在高電阻率異常體模型中，能夠較好反應異常體的頻段范圍，與低電阻率體相比較窄，約在27～212Hz。

圖9 均勻半空間電場歸一化后電斷面圖

2．3 異常體與收發(fā)位置關系對頻率域地－井垂直電磁法的影響

選取低電阻率異常體模型，各項模型參數(shù)圖10，藍色虛線分別代表井的不同位置。用上述有限元方法和網(wǎng)格剖分方式對異常體進行數(shù)值模擬。

圖10 埋深200m低電阻率異常體模型

圖11為歸一化后的地電斷面圖，可以看到較為明顯的異常帶，在異常體周圍呈現(xiàn)較低的電場值，在高頻段異常響應較強，隨著頻率的增加，異常極值依然逐漸向深部轉(zhuǎn)移；對比圖11（a）與（b），隨著收發(fā)距的增加，異常體的異常響應越來越弱（歸一化值往1靠近），異常等值線的密集帶逐漸減稀疏，反應異常體的能力減弱，縱向的分辨率變差（見圖11）。因此，在實際應用時要選擇一個較小的收發(fā)距。

圖11 歸一化后地電斷面圖

3 結(jié) 論

（1）異常體的響應結(jié)果在圖中都表現(xiàn)較為明顯，說明了正演算法的有效性和正確性，算法中引入壓縮存儲數(shù)組，減少了內(nèi)存占有量，便于迭代算法解方程組時程序的調(diào)用。

（2）隨著頻率的增加，電場的異常極值逐漸向深部移動，較高頻段對異常體的反應較強，在實際應用中，盡可能選擇一個合適的頻率段；隨著收發(fā)距的增加，異常逐漸減弱，異常等值線的密集帶逐漸減稀疏，縱向分辨的能力變差，實際應用時要選擇一個較小的收發(fā)距。

（3）模擬研究為地－井電磁法三維反演奠定了基礎。

（4）頻率域地－井電磁法響應特征規(guī)律復雜，本文只考慮了異常體電阻率的大小、埋深、異常體與收發(fā)位置關系對其異常產(chǎn)生的影響，模型的尺寸、形狀以及各種組合異常體同樣需要考慮；選取的場源及模型都是在二維情況下計算的，為了得到更好的模擬效果，需選擇更適當?shù)膱鲈催M行三維模擬提高模擬精度。

（5）在地井觀測方式下，為了突出異常體的響應，減小總場對二次場的影響，采用的是利用背景場進行歸一化的方法。在實際的工作中很難準確得到背景場的值。為了更有效地在總場中提取微弱的二次場異常是一個更重要的問題，需要更進一步深入的研究。

［1］楊曦，潘和平．井間電磁場數(shù)值模擬及成像技術［J］．物探與化探，2009，33（2）：140－147．

［2］何展翔．圈定油氣藏邊界的井地電法研究［D］．成都：成都理工大學，2006．

［3］王志剛，何展翔，魏文博，等．井地電法的準解析近似三維反演研究［J］．石油地球物理勘探，2007，42（2）：220－225．

［4］Augustin A M，Kennedy W D，Morrison H F，et al．A Theoretical Study of Surface to Borehole Electromagnetic Logging in Cased Holes［J］．Gephysics，1989，54（1）：90－99．

［5］Wilt M J，Alumbaugh D L，Morrison H F，et al．Croswwell Electromagnetic Tomography：System Design Considerations and Field Results［J］．Gephysics，1995，60（3）：871－885．

［6］魏寶君，張庚驥．地面－井孔電磁系統(tǒng)成像方法研究［J］．石油大學學報：自然科學版，1999，23（6）：24－29．

［7］李靜和，何展翔．地井垂直電磁 Walkaway剖面法油藏開發(fā)三維模型電場響應特征［J］．石油地球物理勘探，2012，47（4）：653－666．

［8］周峰，潘和平，文國軍，等．基于有限差分的井中激發(fā)極化三維正演［J］．電波科學學報，2010，25（4）：785－792．

［9］孟慶鑫，潘和平．地－井瞬變電磁響應特征數(shù)值模擬分析［J］．地球物理學報，2012，55（3）：1046－1053．

［10］孫向陽，聶在平，李愛勇，等．用高階疊層矢量有限元法計算隨鉆測井的三維電磁響應［J］．電波科學學報，2009，24（2）：273－279．

［11］劉申芬，李天成，慕洪濤．三維井地電阻率有限元數(shù)值模擬及反演［J］．物探與化探，2009，33（1）：88－90．

［12］李琳，李麗娟，高長征，等．電介質(zhì)半空間兩耦合線天線的時域積分方程法［J］．電波科學學報，2008，23（3）：411－416．

［13］王若，王妙月，底青云．頻率域線源大地電磁法有限元正演模擬［J］．地球物理學報，2006，49（6）：1858－1866．

［14］李杰．頻率域可控源二維有限元正演數(shù)值模擬［D］．長沙：中南大學，2010．

［15］柳建新，蔣鵬飛，童孝忠，等．不完全LU分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用［J］．中南大學學報：自然科學版，2009，40（2）：484－491．