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

井间电磁成像的迭代反演算法魏宝君张庚骥曾文冲摘要提出一种基于逐次逼近解法的迭代反演算法，对井间地层电导率的构造进行成像该方法用一阶Born近似将积分方程线性化，得到对电导率分布的初始估计，在迭代反演中用高阶Born近似对井间地层电导率构造进行更精确的估计应用该方法还可以对井间电导率分布进行二次成像，从而使成像分辨率更为准确数值计算结果表明，这种迭代反演算法与基于Born近似、二阶Born近似和扩展Born近似的反演算法相比提高了成像分辨率，且计算效率相当关键词井间电磁成像，迭代反演算法，Born近似，像素THE ITERATIVE INVERSION ALGORITHM IN CROSSWELL ELECTROMAGNETIC IMAGINGWEI BAO-JUN ZHANG GENG-JI(Department of Resources, University of Petroleum, Dongying, Shandong 257062, China)ZENG WEN-CHONG(Shengli Petroleum Administration Burea, Dongying, Shandong 257060, China)AbstractAn iterative inversion method based on SAM(Successive Approximation Method) is developed to image the Earths conductivity structure between two boreholes. The method employs the first Born approximation to linearize the integral equation formulation, and an initial estimate of the conductivity is obtained. A higher-order Born approximation is applied in an iterative manner to achieve better estimates of the structure.Themethod can also be applied to re-image the Earths conductivit struture, and the resolution is much better. Numerical results indicate that the iterative inversion method improves resolution compared with other inversion methods besed on Born approximation or second Born approximation of extended Born approximation, and has the same computational efficiency.Key wordsCrosswell electromagnetic imaging, Iterative inversion algorithm, Bron approximation, Pixel.1引言井间电磁成像的未知量多，像素的数目大，对模型进行严格的正演计算耗费机时，所以一般采用近似成像方法加快反演过程中正演计算的速度已知的近似成像方法有周强等1的Born近似方法Alumbaugh和Morrison2的迭代Born反演方法，Torres-Verdin和Habashy3,4的扩展Born近似反演方法Born近似和迭代Born反演方法只适用于目标电导率与背景电导率比较接近的情况，并且是对小散射体进行低频探测扩展Born近似在正演计算中采用局域非线性方法，是对Born近似的改进，可用于高频和地层电导率对比度较大的情况，但计算精度受到限制Alumbaugh和Morrison在文献5中将Born近似和扩展Born近似结合进行反演他们首先用Born近似估计异常电导率的初始分布，并用局域非线性方法计算散射场，然后用扩展Born近似反演得到对异常电导率的精确估计该过程反复迭代，直至模型的散射场与实测数据吻合到预先给定的精度为止本文参考上述3种近似成像方法，提出了一种基于逐次逼近解法的迭代反演算法该算法既用到了一阶，Born近似和二阶Born近似，也用到了高次Born近似，近似程度较高该算法也适用于较高频率的探测和目标电导率与背景电导率对比度较高的情况，还可以进行二次成像，从而大大提高了成像质量2基于逐次逼近解法的迭代反演算法井间二维问题模型与文献5相同发射电流随时间的变化关系为eit，为发射源的圆频率，在柱坐标系中描述地层磁矢势的积分方程为(1)式中，A为地层中的磁矢势，Ab为磁矢势的背景值，r为接收点的位置坐标，rs为发射源的位置坐标，为电导率异常体所在的区域，r为内任意一点的位置坐标，为二维Green函数，k2ri, k2bib, 为地层磁导率，为地层电导率，b为电导率的均匀背景值(2)式中a为发射线圈半径，nT为发射线圈匝数，I为发射线圈的电流强度(3)其中，J1为1阶Bessel函数，为积分变量反演过程中将井间媒质分成N个相同的矩形像素j，j1，2，N每个像素地层电导率均匀，设为j矩形像素的形状见文献6，其面积为S，节点的局部编号为k1，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列，像素总数目NN1N2用i表示接收点的位置序号，也是像素节点的纵向位置序号，i1，2，N11，用l表示像素节点的横向位置序号，l1，2，N21令对于i2，3，N1，j(i1)N21，N，有(9a)G(i,l,j,k)G(1,l,j(i1)N2,k),(9b)对于i2，3，N11，j1，2，(i1)N2，若mod(j，N2)0，则(10a)(10b)其中，“mod”表示取余当k1时，k3；当k2时，k4；当k3时，k1；当k4时，k2若mod(j，N2)0，则(11a)(11b)在第2次成像过程中，根据像素的划分特点，上述关系式需重新确定4正则化最小二乘解法在每一次迭代成像过程中将(7)式的各方程表示为离散线性方程形式(12)式中piBs(r,rs)表示测量的散射磁场，i表示发射源位置rs与接收位置r的第i次组合，是目标函数，s(n)ij是灵敏度函数，由下式给出(13)(12)式写成矩阵形式为PS.O,(14)式中P为M维复数列向量，S为MN维复数矩阵，O为N维实数列向量解上述线性方程组，得到每个像素电导率的第n次近似(14)式一般是病态的，需用正则化方法求解令下列误差函数取极小值(O)|S.OP|2|O|2,(15)式中是正则化因子(15)式的极小等价于下面线性方程组Re(SS)I.ORe(SP),(16)式中“Re”表示取实部，上角“”表示共轭转置，I为单位矩阵由于矩阵S和列向量P的元素都是复数，在(16)式中需有取实部的符号(推导过程见附录B)在求解(16)式时，为便于使取初始值，将矩阵Re(SS)的对角线元素作归一化处理5数值计算本文迭代反演算法的有效性可通过下面几个例子说明在所有例子中，井间距离均为100m，垂向测量范围也是100m，数值模拟的实测数据均含0.5%的噪声发射装置和接收装置的垂向测量间隔为5m，发射-接收之间共有441种组合，两井间共有400个像素，每个像素的尺寸为5m5m为避免计算过程中奇异现象的出现，第1次成像中未对最靠近发射井的一列像素进行成像这样(14)式就是含441个方程和380个未知量的超定方程组第1个模型中两井间有一个高电导率异常体，电导率为0.1S/m，尺寸为10m10m，背景电导率0.01S/m，发射频率50kHz由图1(b)可以确定出电导率异常体的分布区域及第2次成像的电导率初始值电导率异常区域确定为横向范围4060m，纵向范围4060m，每个像素的尺寸为2m2m，迭代初始值确定为0.03S/m图1井间电磁成像结果(a)地层模型；(b)一次成像结果；(c)二次成像结果Fig.1Results of crosshole electromagnetic imaging第2个模型中有一个低电导率异常体，其电导率为0.2S/m，尺寸是10m10m，背景电导率1.0S/m，发射频率1kHz由图2(b)确定出电导率异常区域的横向范围为4060m，纵向范围为4060m，每个像素的尺寸为2m2m，第2次成像的迭代初始值确定为0.8S/m图2井间电磁成像结果(a)地层模型；(b)一次成像结果；(c)二次成像结果Fig.2Results of crosshole electromagnetic imaging第3个模型中在垂直方向有两个高电导率异常体，电导率均为0.5S/m，尺寸均为15m6m，背景电导率0.1S/m，发射频率为10kHz由图3(b)确定出电导率异常区域的横向范围为3565m，纵向范围为3070m，每个像素的尺寸为3m2m，第2次成像的迭代初始值确定为0.35S/m图3井间电磁成像结果(a)地层模型；(b)一次成像结果；(c)二次成像结果Fig.3Results of crosshole electromagnetic imagingBorn近似、迭代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.OP)(S.OP)OTO,(B1)式中“T”表示转置，由于O是实列向量，对它共轭转置等价于转置由于O使得(O)达到极小，故应有式中微分算子是列向量算子，即它作用到(O)以后，形成列向量(B2)列向量算子作用到行向量OT(o1,o2,oN)上形成单位矩阵I它不能作用到列向量O上，为了对含有列向量O的项求导需将该项转置(B3)“*”表示共轭又因为SS与STS*互为共轭，SP与STP*互为共轭，所以(B4)令该式等于零，就得到公式(16)作者简介：魏宝君，男，1969年生，1991年毕业于北京师范大学物理系1997年获石油大学应用地球物理专业硕士学位现从事教学和电法测井理论研究工作作者单位：魏宝君张庚骥石油大学(华东)资源系，山东东营 257062曾文冲胜利石油管理局，山东东营 257060参考文献1Zhou Q, Becker A, Morrison H F. Audio-frequency electromagnetic tomography in 2-D.Geophysics, 1993, 58(4):4824952Alumbaugh D L, Morrison H F. Electromagnetic conductivity imaging with an iterative Born inversion. IEEE Trans. on Geosci. Remote Sensing, 1993,31(4):7587633Torres-Verdin C, Habashy T M. Rapid 2.5-D forward model