井间电磁成像的迭代反演算法.docx

井间电磁成像的迭代反演算法.docx

《井间电磁成像的迭代反演算法.docx》由会员分享，可在线阅读，更多相关《井间电磁成像的迭代反演算法.docx（10页珍藏版）》请在冰豆网上搜索。

井间电磁成像的迭代反演算法

井间电磁成像的迭代反演算法

魏宝君　张庚骥　曾文冲

摘　要　提出一种基于逐次逼近解法的迭代反演算法，对井间地层电导率的构造进行成像．该方法用一阶Born近似将积分方程线性化，得到对电导率分布的初始估计，在迭代反演中用高阶Born近似对井间地层电导率构造进行更精确的估计．应用该方法还可以对井间电导率分布进行二次成像，从而使成像分辨率更为准确．数值计算结果表明，这种迭代反演算法与基于Born近似、二阶Born近似和扩展Born近似的反演算法相比提高了成像分辨率，且计算效率相当．

关键词　井间电磁成像，迭代反演算法，Born近似，像素．

THEITERATIVEINVERSIONALGORITHMINCROSSWELLELECTROMAGNETICIMAGING

WEIBAO-JUNZHANGGENG-JI

（DepartmentofResources,UniversityofPetroleum,Dongying,Shandong257062,China）

ZENGWEN-CHONG

（ShengliPetroleumAdministrationBurea,Dongying,Shandong257060,China）

Abstract　AniterativeinversionmethodbasedonSAM（SuccessiveApproximationMethod）isdevelopedtoimagetheEarth’sconductivitystructurebetweentwoboreholes.ThemethodemploysthefirstBornapproximationtolinearizetheintegralequationformulation,andaninitialestimateoftheconductivityisobtained.Ahigher-orderBornapproximationisappliedinaniterativemannertoachievebetterestimatesofthestructure.Themethodcanalsobeappliedtore-imagetheEarth’sconductivitstruture,andtheresolutionismuchbetter.NumericalresultsindicatethattheiterativeinversionmethodimprovesresolutioncomparedwithotherinversionmethodsbesedonBornapproximationorsecondBornapproximationofextendedBornapproximation,andhasthesamecomputationalefficiency.

Keywords　Crosswellelectromagneticimaging,Iterativeinversionalgorithm,Bronapproximation,Pixel.

1　引　言

　　井间电磁成像的未知量多，像素的数目大，对模型进行严格的正演计算耗费机时，所以一般采用近似成像方法加快反演过程中正演计算的速度．已知的近似成像方法有周强等〔1〕的Born近似方法Alumbaugh和Morrison〔2〕的迭代Born反演方法，Torres-Verdin和Habashy〔3,4〕的扩展Born近似反演方法．Born近似和迭代Born反演方法只适用于目标电导率与背景电导率比较接近的情况，并且是对小散射体进行低频探测．扩展Born近似在正演计算中采用局域非线性方法，是对Born近似的改进，可用于高频和地层电导率对比度较大的情况，但计算精度受到限制．

　　Alumbaugh和Morrison在文献〔5〕中将Born近似和扩展Born近似结合进行反演．他们首先用Born近似估计异常电导率的初始分布，并用局域非线性方法计算散射场，然后用扩展Born近似反演得到对异常电导率的精确估计．该过程反复迭代，直至模型的散射场与实测数据吻合到预先给定的精度为止．

　　本文参考上述3种近似成像方法，提出了一种基于逐次逼近解法的迭代反演算法．该算法既用到了一阶，Born近似和二阶Born近似，也用到了高次Born近似，近似程度较高．该算法也适用于较高频率的探测和目标电导率与背景电导率对比度较高的情况，还可以进行二次成像，从而大大提高了成像质量．

2　基于逐次逼近解法的迭代反演算法

　　井间二维问题模型与文献〔5〕相同．发射电流随时间的变化关系为eiωt，为发射源的圆频率，在柱坐标系中描述地层磁矢势的积分方程为

（1）

式中，A为地层中的磁矢势，Ab为磁矢势的背景值，r为接收点的位置坐标，rs为发射源的位置坐标，为电导率异常体所在的区域，r′为内任意一点的位置坐标，为二维Green函数，k2r＝－iωμσ,k2b＝－iωμσb,为地层磁导率，为地层电导率，b为电导率的均匀背景值．

（2）

式中a为发射线圈半径，nT为发射线圈匝数，I为发射线圈的电流强度．

（3）

其中，

J1为1阶Bessel函数，λ为积分变量．

　　反演过程中将井间媒质分成N个相同的矩形像素Ωj，j＝1，2，…，N．每个像素地层电导率均匀，设为σj．矩形像素的形状见文献〔6〕，其面积为S，节点的局部编号为k＝1，2，3，4．过像素内任意一点且分别与r轴和z轴平行的直线将像素分割为四部分，每一部分的面积为Sk，定义面积坐标

（4）

其中r0、z0为像素中心点的径向和纵向坐标，Δr、Δz为像素的径向和纵向宽度，上标“int”表示取整．像素内的磁矢势由下面的插值函数给出

（5）

根据

（1）式和（5）式得到井间磁矢势的1阶到n阶近似

（6a）

（6b）

……

（6c）

其中A（0）j（rk,rs）由

（2）式得到．由（6）式得到散射磁场z分量的1阶到n阶近似

（7a）

（7b）

……

（7c）

其中

（8）

式中J0为0阶Bessel函数．

　　根据（6）式和（7）式进行迭代反演的具体方法为：

由（7a）式，用B

（1）s（r,rs）近似测量值Bs（r,rs），应用正则化最小二乘法得到各像素电导率的1阶近似σ

（1）j，并由（6a）式得到各像素节点磁矢势的1阶近似A

（1）j（rk,rs）．将A

（1）j（rk,rs）代入（7b）式，用B

（2）s（r,rs）近似Bs（r,rs），得到各像素电导率的2阶近似σ

（2）j．将σ

（2）j依次回代到（6a）式和（6b）式，得到各像素节点磁矢势的2阶近似A

（2）j（rk,rs）．根据A

（2）j（rk,rs），用B（3）s（r,rs）近似Bs（r,rs），得到σ（3）j．重复上述过程，逐次迭代，当n达到某一数值后，在n阶近似的基础上迭代，直至得到满足要求的σ（n）j．

　　利用本文迭代反演方法还可以对井间电导率分布进行二次成像，即第1次对整个井间范围成像后，确定出电导率异常体的大致区域及该区域内电导率的大致分布，然后对电导率异常体区域进行第2次成像．由于第2次成像范围小，在成像中可以减小像素尺寸，在每一个像素内对磁矢势进行插值使计算更精确．另外所选电导率初始值与第1次成像时选背景电导率为初始值相比更接近于真实值，使该反演方法的精确度大大提高，从而提高了成像分辨率．

3　积分的简化

可先算出，反演中调用这些数据即可．上述两个积分都是三重积分，经化简可变为一重积分（附录A）．

　　上述两个积分要根据不同的接收位置或节点位置对所有像素进行计算，运算量大．根据第1次反演过程中像素的划分特点，通过简化方法可避免对积分进行重复计算，使运算量大大减少．第1次成像过程中，井间像素的垂向划分以接收点的纵坐标为边界，且每个像素的尺寸完全相同．设井间成像区域在纵向被划分为N1层，在横向被划分为N2列，像素总数目N＝N1N2．用i表示接收点的位置序号，也是像素节点的纵向位置序号，i＝1，2，…，N1＋1，用l表示像素节点的横向位置序号，l＝1，2，…，N2＋1．令

对于i＝2，3，…，N1，j＝（i－1）N2＋1，…，N，有

（9a）

G（i,l,j,k）＝G（1,l,j－（i－1）N2,k）,

（9b）

对于i＝2，3，…，N1＋1，j＝1，2，…，（i－1）N2，若mod（j，N2）≠0，则

（10a）

（10b）

其中，“mod”表示取余．当k＝1时，k′＝3；当k＝2时，k′＝4；当k＝3时，k′＝1；当k＝4时，k′＝2．若mod（j，N2）＝0，则

（11a）

（11b）

在第2次成像过程中，根据像素的划分特点，上述关系式需重新确定．

4　正则化最小二乘解法

　　在每一次迭代成像过程中将（7）式的各方程表示为离散线性方程形式

（12）

式中pi＝Bs（r,rs）表示测量的散射磁场，i表示发射源位置rs与接收位置r的第i次组合，

是目标函数，s（n）ij是灵敏度函数，由下式给出

（13）

（12）式写成矩阵形式为

P＝S.O,

（14）

式中P为M维复数列向量，S为M×N维复数矩阵，O为N维实数列向量．解上述线性方程组，得到每个像素电导率的第n次近似．

　　（14）式一般是病态的，需用正则化方法求解．令下列误差函数取极小值

Φμ（O）＝||S.O－P||2＋μ||O||2,

（15）

式中是正则化因子．（15）式的极小等价于下面线性方程组

〔Re（S＋S）＋μI〕.O＝Re（S＋P）,

（16）

式中“Re”表示取实部，上角“＋”表示共轭转置，I为单位矩阵．由于矩阵S和列向量P的元素都是复数，在（16）式中需有取实部的符号（推导过程见附录B）．在求解（16）式时，为便于使取初始值，将矩阵Re（S＋S）的对角线元素作归一化处理．

5　数值计算

　　本文迭代反演算法的有效性可通过下面几个例子说明．在所有例子中，井间距离均为100m，垂向测量范围也是100m，数值模拟的实测数据均含0.5%的噪声．发射装置和接收装置的垂向测量间隔为5m，发射-接收之间共有441种组合，两井间共有400个像素，每个像素的尺寸为5m×5m．为避免计算过程中奇异现象的出现，第1次成像中未对最靠近发射井的一列像素进行成像．这样（14）式就是含441个方程和380个未知量的超定方程组．

　　第1个模型中两井间有一个高电导率异常体，电导率为0.1S/m，尺寸为10m×10m，背景电导率0.01S/m，发射频率50kHz．由图1（b）可以确定出电导率异常体的分布区域及第2次成像的电导率初始值．电导率异常区域确定为横向范围40—60m，纵向范围40—60m，每个像素的尺寸为2m×2m，迭代初始值确定为0.03S/m．

图1　井间电磁成像结果

（a）地层模型；（b）一次成像结果；（c）二次成像结果．

Fig.1　Resultsofcrossholeelectromagneticimaging

　　第2个模型中有一个低电导率异常体，其电导率为0.2S/m，尺寸是10m×10m，背景电导率1.0S/m，发射频率1kHz．由图2（b）确定出电导率异常区域的横向范围为40—60m，纵向范围为40—60m，每个像素的尺寸为2m×2m，第2次成像的迭代初始值确定为0.8S/m．

图2　井间电磁成像结果

（a）地层模型；（b）一次成像结果；（c）二次成像结果．

Fig.2　Resultsofcrossholeelectromagneticimaging

　　第3个模型中在垂直方向有两个高电导率异常体，电导率均为0.5S/m，尺寸均为15m×6m，背景电导率0.1S/m，发射频率为10kHz．由图3（b）确定出电导率异常区域的横向范围为35—65m，纵向范围为30—70m，每个像素的尺寸为3m×2m，第2次成像的迭代初始值确定为0.35S/m．

图3　井间电磁成像结果

（a）地层模型；（b）一次成像结果；（c）二次成像结果．

Fig.3　Resultsofcrossholeelectromagneticimaging

　　Born近似、迭代Born近似和扩展Born近似这3种反演方法是本文迭代反演方法的特殊情况．利用本文反演方法进行井间电磁第1次成像能够较准确地确定电导率异常体的范围，第2次成像可以在此基础上准确地确定出电导率值．该方法与上述3种反演方法相比计算效率相当，但成像效果更好．

附录A　积分的简化

第一个积分可表示为

（A1）

（A2）

（A3）

φ1（r0,Δr,λ）内包含着两种积分，其中一种积分为

（A4）

另一种积分可简化为〔6〕

（A5）

其中，H0和H1分别是0阶和1阶Struve函数．

　　简化（A2）式和（A3）式，经计算得到

（A6）

当

时

（A7）

当

时

（A8）

当

时

（A9）

　　在计算Struve函数时，若|x|≤20.0，采用小宗量渐近展开，若|x|≥20.0，采用大宗量展开（文献〔7〕）．对于

只需将第一个积分中的J1（λr）换成J0（λr）即可．

附录B　公式（16）的推导

　　由（15）式，

Φμ（O）＝（S.O－P）＋（S.O－P）＋μOTO,

（B1）

式中“T”表示转置，由于O是实列向量，对它共轭转置等价于转置．

　　由于O使得Φμ（O）达到极小，故应有

式中微分算子是列向量算子，即

它作用到Φμ（O）以后，形成列向量．

（B2）

列向量算子

作用到行向量OT＝（o1,o2,……,oN）上形成单位矩阵I．它不能作用到列向量O上，为了对含有列向量O的项求导需将该项转置．

（B3）

“*”表示共轭．又因为S＋S与STS*互为共轭，S＋P与STP*互为共轭，所以

（B4）

令该式等于零，就得到公式（16）．

作者简介：

魏宝君，男，1969年生，1991年毕业于北京师范大学物理系．1997年获石油大学应用地球物理专业硕士学位．现从事教学和电法测井理论研究工作．

作者单位：

魏宝君　张庚骥　石油大学（华东）资源系，山东东营257062

　　　　　曾文冲　胜利石油管理局，山东东营257060

参考文献

1　ZhouQ,BeckerA,MorrisonHF.Audio-frequencyelectromagnetictomographyin2-D.Geophysics,1993,58（4）:

482—495

2　AlumbaughDL,MorrisonHF.ElectromagneticconductivityimagingwithaniterativeBorninversion.IEEETrans.onGeosci.RemoteSensing,1993,31（4）:

758—763

3　Torres-VerdinC,HabashyTM.Rapid2.5-Dforwardmodelingandinversionviaanewnonlinearscatteringapproximation.RadioScience,1994,29（4）:

1051—1079

4　Torres-VerdinC,HabashyTM.Atwo-steplinearinversionoftwo-dimensionalelectricalconductivity.IEEETrans.onAntennasPropagat.,1995,43（4）:

405—415

5　AlumbaughDL,MorrisonHF.Theoreticalandpracticalconsiderationsforcrosswellelectromagnetictomographyassumingacylindricalgeometry.Geophysics,1995,60（3）:

846—870

6　张庚骥．电法测井（下册）．北京：

石油工业出版社，1986

　　ZHANGGeng-Ji.ElectricLogging（II）（inChinese）.Beijing:

PetroleumIndustryPress,1986

7　AbramowitzM,StegunIA.HandbookofMathematicalFunctions.NewYork:

DoverPublications,1972.496—498