地球物理学报  2016, Vol. 59 Issue (8): 3110-3120   PDF    
引用本文
潘克家, 汤井田, 杜华坤 等. 2016. 轴对称地层中高分辨率阵列侧向测井信赖域反演法. 地球物理学报, 59(8): 3110-3120, doi: 10.6038/cjg20160833.  
PAN K J, TANG J T, DU H K, et al. 2016. Trust region inversion algorithm of high-resolution array lateral logging in axisymmetric formation. Chinese J. Geophys. (in Chinese), 59(8): 3110-3120, doi: 10.6038/cjg20160833.  
轴对称地层中高分辨率阵列侧向测井信赖域反演法
潘克家1,2 , 汤井田3 , 杜华坤3 , 蔡志杰4     
1. 中南大学数学与统计学院, 长沙 410083;
2. 油气藏地质及开发工程国家重点实验室(成都理工大学), 成都 610059;
3. 中南大学有色金属成矿预测教育部重点实验室, 地球科学与信息物理学院, 长沙 410083;
4. 复旦大学数学科学学院, 上海 200433
收稿日期 2015-05-16, 2016-05-24 收修定稿
基金项目: 国家自然科学基金项目（41474103，41204082），国家863计划项目（2014AA06A602），湖南省自然科学基金项目（2015JJ3148）和成都理工大学油气藏地质及开发工程国家重点实验室资助项目（PLC201305）联合资助.
作者简介: 潘克家,男,1981年生,湖南宁乡人,2004年毕业于中国石油大学,2009年在复旦大学获得理学博士学位,现为中南大学副教授,主要从事地电磁场正反演研究.E-mail:pankejia@hotmail.com
通讯作者: 杜华坤,男,1977年生,云南镇雄人,博士,中南大学讲师,主要从事地球物理数据处理与正反演研究.E-mail:hkdu2008@126.com
摘要: 本文研究轴对称地层中高分辨率阵列侧向测井（HRLA）的多参数信赖域反演方法.首先对前期HRLA的有限元正演方法进行改进，提出基于叠加原理和并行直接稀疏求解器PARDISO的快速正演方案，更适合于反演计算.将HRLA反演问题转化为非线性最小二乘问题，利用信赖域算法求解.为提高反演速度，推导了目标函数对优化参数偏导数的具体计算公式.对典型地层模型，与已有文献中Jacobi预条件共轭梯度法（JCG）计算结果比较，发现PARDIDO比JCG快10倍以上，验证了本文正演程序的正确性与高效性.利用信赖域算法求解了电阻率四参数反演和传统的三参数反演.研究结果表明：并行直接稀疏求解器PARDISO能有效求解此类HRLA正演问题，对6次不同探测深度的测井模拟，所形成的有限元刚度矩阵完全相同，只须进行一次矩阵分解，大大加快了正反演的速度.信赖域算法收敛速度快，且具有全局收敛性. HRLA的信赖域反演结果几乎不依赖于初值的选取，从较差初值出发仍能得到满意的反演结果.另外信赖域算法抗噪能力比较强，即使对测井数据添加信噪比为10dB（甚至5dB）的高斯白噪声，仍能通过反演得到较为准确的地层参数.
关键词: 阵列侧向测井      信赖域算法      反演      全局收敛      PARDISO     
Trust region inversion algorithm of high-resolution array lateral logging in axisymmetric formation
PAN Ke-Jia1,2, TANG Jing-Tian3, DU Hua-Kun3, CAI Zhi-Jie4     
1. School of Mathematics and Statistics, Central South University, Changsha 410083, China;
2. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation(Chengdu University of Technology), Chengdu 610059, China;
3. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education/School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
4. School of Mathematical Sciences, Fudan University, Shanghai 200433, China
Abstract: The high-resolution laterolog array (HRLA) tool can simultaneously provide six curves of logging responses with different depths of investigation, which give abundant information for the inversion of geological formation parameters. This paper proposes a trust region algorithm to recover several parameters in the axisymmetric formation from HRLA data, including the mud-invasion radius, mud resistivity and original formation resistivity. Firstly, we propose an accurate and efficient method to solve the forward problem in order to obtain the logging responses of HRLA tool in the axisymmetric formation, where the superposition principle is used to simplify the calculation of resistivity logging responses, and the parallel direct solver PARDISO is employed to solve the linear system with multiple right-hand sides resulting from finite element approximation. The decomposition of the stiffness matrix is needed only once when using PARDISO to solve these linear systems, which is cheap in terms of work and fast in time. Secondly, the problem of identifying formation parameters is transformed into a nonlinear least squares problem that can be solved by the trust region algorithm. In order to speed up the optimization process, we derive an explicit formula for the partial derivative of the objective function with respect to each resistivity parameter. Finally, for a typical formation model we present resistivity four-parameter and traditional three-parameter inversion results by using different initial guesses.In order to test the validity and effectiveness of the proposed method, we design a typical formation model with borehole radius rh=0.1016 m, invasion radius ri=0.4016 m, permeable formation thickness H=2 m, borehole mud resistivity Rm=1 Ωm, surrounding rock resistivity Rs=10 Ωm and original formation resistivity Rt=500 Ωm. Comparing the forward modeling result with the result in the existing literature, we find that the maximum relative error between the two results is only 0.01%, which verifies the correctness of the forward algorithm program. It takes less than half a second to perform forward modeling with 24881 grid points using PARDISO on a personal laptop. The computational time of the Jacobi conjugate gradient (JCG) method is more than ten times of that of PARDISO, and PARDISO provides a more accurate solution than the JCG. We use the trust region algorithm to invert formation parameters with four different initial values, one of which is far away from the real model. The inversion results show that the method can converge to the true model quickly even with very poor initial guesses. When using noisy data mixed with Gaussian white noise with signal-to-noise ratio of 10dB or even 5dB, the trust inversion algorithm still gives satisfactory results. Trust region algorithm combined with PARDISO is successfully used to multi-parameter inversion of HRLA data. Numerical results show that the trust region inversion algorithm converges fast, and has a global convergence property. Another advantage of the proposed method is the robustness to noise..
Key words: Array lateral logging      Trust region algorithm      Inversion      Global convergence      PARDISO     
1 引言

电法测井是测井方法中一类重要的方法，其目的是通过求取地层真电阻率来确定地层孔隙度、渗透率及含油气饱和度等重要地质参数. 而电阻率测井受钻井泥浆、井眼、围岩等因素的影响，电阻率测井的读数与地层真实电阻率之间存在偏差. 为消除这些影响，恢复地层真实电阻率，电法测井反演成为首选. 传统的双侧向测井仅能提供深、浅两条电阻率测井曲线，测量地层信息少，无法同时反演3个以上地层参数(如侵入带半径、侵入带电阻率、原始地层电阻率等). 1998年，斯伦贝谢(Schlumberger)公司推出高分辨阵列侧向测井仪器(High-Resolution Lateral Array tool，HRLA)，该测井仪能同时提供5～6条不同探测深度的测井响应曲线，具有纵向分辨率高、径向探测信息丰富等优点，为精确反演地层参数提供了丰富的信息(Chen et al.，1998;张建华等，2002). 因此，对阵列侧向测井反演方法的研究，不仅能为测井仪器的研发提供理论依据，且能指导实际测井工程，具有重要的理论和现实意义.

众所周知，正演是反演的基础．电法测井的正演方法主要包括有限单元法(FEM)、有限差分法(FDM)、数值模拟匹配法(NNM). FDM 主要用于自然电位测井以及感应测井的数值模拟计算(Weller et al.，2003;汪功礼等，2003;沈金松，2004;潘克家和谭永基，2009;Pan et al.，2009;杨守文等，2009;Su et al.，2012). NNM方法为一种半数值、半解析的高效算法，最初用于分析二维非均匀介质中的电磁散射，后经Chew、聂在平、张庚骥、汪宏年等人的完善，已成功应用于各种测井响应的数值分析中(Chew et al.，1991;聂在平和陈思渊，1994;张庚骥和汪涵明，1996;袁宁等，1998;Fan et al.，2000;谭茂金等，2007;汪宏年等，2008;Dun et al.，2010). 有限单元法在电法测井数值模拟中的研究和应用较为成熟．1980年，李大潜的《有限元素法在电法测井中的应用》，是国内最早的有限元法在电法中应用的专著．之后，张庚骥等人进行了一系列具体和深入的工作，将FEM 法推向实用(李大潜等，1980;张庚骥，1984)．

阵列侧向测井仪包含25个电极(见图 1)，具有6种不同的探测深度的发射方式，相比双侧向测井仪，结构更加复杂. 而FEM处理复杂边界、适应问题能力比较强. 因此，目前阵列侧向测井的正演主要采用FEM法. 刘振华和胡启(2002)、仵杰等(2008)、范宜仁等(2009)对于包含井眼、侵入带、围岩和目的层的二维轴对称地层模型，利用FEM法研究了不同层厚、侵入带、井眼条件下阵列侧向测井响应. 邓少贵和李智强(2009)、邓少贵(2010a)、Deng等(2011)采用裂缝组平板裂缝模型以及改进的前线解法，利用三维FEM进行不同裂缝参数条件下的高分辨率阵列侧向测井响应数值模拟. 邓少贵等(2010b)采用三维FEM对倾斜井进行了阵列侧向测井响应的数值模拟，并利用Levenberg-Marquardt(LM)方法进行阵列侧向测井的四参数反演. 潘克家等(2013)建立了高分辨率阵列侧向测井的等值面边值问题数学模型，并采用Jacobi预条件共轭梯度(JCG)法求解有限元方程，提高了测井响应的计算速度．祝鹏等(2015)采用三维有限元法模拟了水平井和大斜度井中阵列侧向测井响应. 值得指出的是，前面研究所提到的阵列侧向测井仪稍有不同. 刘振华和胡启(2002)、潘克家等(2013)和祝鹏等(2015)文中考虑1个主电极A0、6对屏蔽电极A1，A2，…，A6和6对监督电极M1，M2，…，M6的阵列侧向测井仪;而其余文献均研究只含2对监督电极M1，M2的阵列侧向测井仪.

图 1 轴对称地层模型示意图 Rm，Rx0，Rt，Rs分别为泥浆、侵入带、原始地层和围岩电阻率;rh，ri分别为井眼半径和侵入带半径，H为原始地层厚度;z∞，r∞为适当大的有界化常数. Fig. 1 chematic diagram of axisymmetric formation Rm，Rx0，Rs and Rt are resistivities of the borehole mud，invaded zone，surrounding rock and original formation，respectively;rh and ri the radiuses of borehole and invaded zone;H the thickness of the formation;z∞ and r∞ suitably large positive constants.

基于正演研究工作，已有不少学者开始研究阵列侧向测井的多参数反演问题. Zhang和Zhou(2002)采用降维技术与神经网络算法研究了阵列侧向测井的拟二维实时反演. 文中提到的测井仪由1个发射电极和18个测量电极组成，同时测量8个电压数据和16个相邻测量电极的一阶差分数据. 杨韡(2003)针对具有4种不同探测深度的阵列侧向测井仪(含主电极A0和4对屏蔽电极A1，A2，A3，A4)，采用LM算法求解传统的三参数(侵入带半径、侵入带电阻率和原始地层电阻率)反演问题. 刘振华和张霞(2005)采用LM算法进行了阵列侧向测井的四参数(冲洗带电阻率、原始地层电阻率、冲洗带半径和过渡带半径)反演. 顿月芹和袁建生(2009)基于点源模拟柱状电极的快速NMM正演(Dun et al.，2010)，采用BP神经网络实现了阵列侧向测井的三参数非线性反演. 李智强等(2010)利用差分进化算法，对阵列侧向测井响应的三参数反演进行了初步尝试. LM算法结合了高斯牛顿法和最速下降法的优点，但仍只具有局部收敛性，即初值选择不合适，容易导致算法陷入局部极值导致计算失败. 而智能算法(如差分进化算法)通常具有全局寻优能力，可避开反演初值的影响;但计算量相对较大，寻优过程非常耗时，难以实现快速、实时反演.

信赖域算法是非线性优化中一类非常重要的全局优化算法，近几十年受到了非线性优化界的广泛重视(袁亚湘和孙文瑜，1997). 目前，信赖域方法已和传统的线搜索方法并列成为非线性规划的两类最主要数值方法. 近年来，已有不少学者将信赖域算法应用各种地球物理反演问题，如一维大地电磁反演(曾奇等，2011)、磁共振测深反演(林婷婷等，2013)和地震定位(田玥和陈晓非，2006). 王彦飞课题组已成功将信赖域算法应用于各种地震数据处理，如地震数据重构(Wang et al.，2011;Cao and Wang，2014)和地震偏移反演成像(Li and Wang，2014).

为进一步提高阵列侧向测井资料的反演速度与反演结果的可靠性，本文对作者前期正演工作(潘克家等，2013)进行改进，更适合于反演计算. 进而引入信赖域算法求解轴对称地层条件下阵列侧向测井的多参数反演问题，设计典型地层模型测试正、反演算法，为下一步实现三维阵列侧向测井数据实时反演奠定基础.

2 阵列侧向测井原理及正演 2.1 阵列侧向测井仪结构及工作原理

高分辨阵列侧向测井仪由1个主电极A0、6对屏蔽电极A1，A2，…，A6和6对监督电极M1，M2，…，M6组成，电极关于主电极对称排列. 图 1给出了阵列侧向测井仪上半部分示意图. 测量时，主电极A0发射单位正电流，监督电极Mi，M′i(1≤i≤6)不发射电流. 只有A0发射电流，其他屏蔽电极为回路电极时构成最浅探测深度的响应RLA0. 此时各个屏蔽电极上的电流由如下等电位条件确定: UA3=UA4=UA5=UA6，UM3=UM4，UM5=UM6. 当A0向两侧每次增加一对屏蔽电极为发射电流电极时，可依次构成RLA1，RLA2，RLA3，RLA4和RLA5测井响应，它们具有不同的探测深度和较高的分辨能力.

2.2 阵列侧向测井数学模型

不失一般性，假设上下围岩电导率相等. 考虑轴对称地层模型，对称轴通过电极中心线. 以对称轴为z轴，电极中心为坐标原点，建立如图 1所示坐标系. 则对包含井眼、侵入带、围岩和原始地层的轴对称地层模型，电位函数u(r，z)在截断区域Ω内满足如下椭圆方程边值问题(潘克家等，2013):

(1)

(2)

(3)

(4)

其中，边界Γ1表示无穷远截断边界(r=r∞和z=z∞)，边界Γ2包括中心对称轴r=0、对称面z=0以及电极系表面上的绝缘环;Γ0i表示第i个电极表面，Ci为待定常数，Ii为第i个电极发射的电流;γ为地层不同区域的交界面，“+”和“-”分别表示函数在交界面两侧的取值，n为边界外法线;r0为电极系半径，R为分块常数的电阻率;Γ0i表示主电极A0表面，Γ0i，Γ06+i(1≤i≤6)分别对应屏蔽电极Ai与监督电极Mi的表面，由监督电极不发射电流可知Ii=0(7≤i≤12).边界条件表示绝缘边界、对称边界及截断边界上电流为零;边界条件表示第i个电极为等势体，且发射电流为Ii;内边界条件表示不同地层交界面上电位及电流均连续.

边界条件与已有文献(潘克家等，2013)中不同，前文中在无穷远截断边界上施加齐次混合边界条件

(5)

此条件虽比齐次Dirichlet边界条件u=0(假定电位为零)，及齐次Neumann边界条件(假定电流为零)精确，能确保正演精度的同时缩小计算区域. 然而，反演中采取边界条件并不合适，因为边界条件导致有限元离散形成的刚度矩阵与发射源到边界的距离有关，故对不同电极发射电流时需分别进行刚度矩阵的组装与分解. 若在无穷远截断边界上施加条件，则对不同的电流发射源(对应阵列侧向测井的多次测量)所形成的刚度矩阵完全相同，采用直接法如PARDISO求解时刚度矩阵只需分解一次，能大大提高正演速度.

2.3 阵列侧向测井有限元正演

记基本解uj(0≤j≤6)为屏蔽电极Aj发射单位正电流，其他电极均不发射电流时，边值问题(1)—(4)的解. 记Ij(k)(1≤k≤6，0≤j≤6)为第k次测量时，屏蔽电极Aj发射的电流值，而u(k)(1≤k≤6)表示第k次测量的电位函数. 则由叠加原理知

(6)

其中，I0(k)=0.5，Ij(k)(1≤j≤6)待定，要求满足各次测量时相应的等电位条件及电流代数和为零的条件.

对第k次测量，设等电位条件为u(k)(Pi(k))=u(k)(Qi(k))(i=1，2，…，5)，其中Pi(k)和Qi(k)分别表示某一屏蔽电极或监督电极. 将(6)代入等电位条件并整理可得

(7)

再加上电流代数和为零的条件

(8)

得到一个关于Ij(k)(j=1，2，…，6)的六元线性方程组，求解即得屏蔽电极所发射的电流值.

利用叠加原理，成功避开了屏蔽电极上电流的不确定性;通过计算电流发射情况非常简单的基本解uj(0≤j≤6)，即可同时得到6次不同探测深度的测井曲线. 而对基本解uj的计算，虽然每次电流发射位置不同，但有限元离散形成的刚度矩阵完全相同，故对7个基本解的计算，实际上只需进行一次矩阵分解，进而计算6条测井曲线只须进行一次矩阵分解，大大加快了HRLA正演速度. 边值问题离散后最终得到的有限元矩阵为对称正定的9对角矩阵，本文采取PARDISO快速求解. 具体的有限元分析过程可参考(潘克家等，2013).

2.4 稀疏直接求解器与压缩存储格式

(1) 并行直接求解器

PARDISO是Basel大学针对大规模稀疏线性方程组开发的高效并行直接求解器. 该求解器采用BLAS软件包，实现算法中线性代数操作的并行计算. 大量数值测试表明，PARDISO是目前最快的线性稀疏矩阵求解方法之一(Gould et al.，2007). Intel公司获得授权后，Intel Math Kernel Library(Intel MKL)提供了PARDISO的优化版本，效率高，稳定性好.

PARDISO求解器基于LU分解，融合了向左和向右算法的优点，采用超节点消去树进行填充元优化，降低算法的时空消耗. 算法的具体步骤如下:

1) 矩阵重排与符号分解: 针对稀疏矩阵A的结构，设计合适的行列交换矩阵，对原矩阵进行交换与重排，使得新矩阵Ã分解后含有尽量少的非零元素.

2) 矩阵LU分解: 对新矩阵Ã进行LU分解.

3) 方程求解与迭代加精: 利用LU分解结果，求解方程;如对结果精度有进一步要求，使用迭代法进一步提高精度.

4) 迭代结束，释放求解器所占内存.

(2) 压缩稀疏行格式

PARDISO求解器采用目前广泛流行的压缩稀疏行格式(Compressed Sparse Row，CSR)对稀疏矩阵进行压缩存储，大大降低了矩阵元素的访问时间和空间存储成本. 此种格式按行依次存储矩阵A的非零元. 设nnz为n阶对称方阵A上三角部分非零元的个数，则这种存储结构由以下三个数组构成:

1) 双精度数组a——A的上三角部分的非零元素，长度为nnz;

2)整型数组Ja——A中元素的列号，长度为nnz;

3) 整型数组Ia——A中每行第一个非零元在数组a中的位置，长度为n+1，其中Ia(n+1)=nnz+1.

PARDISO采用的CSR格式的存储效率比潘克家等(2012)中采用的基于地址矩阵的压缩存储效率更高. 另外，CSR格式非常灵活，便于操作，如计算矩阵向量乘法Ax的第i个分量为

(9)

图 2为一个对称CSR格式的存储示例. 阵列侧向测井正演问题(1)—(4)通过有限元离散后，得到的刚度矩阵为稀疏、对称矩阵，可采用对称CSR格式存储.

图 2 CSR格式示例 Fig. 2 An example of CSR format
3 阵列侧向测井反演 3.1 阵列侧向测井反演

参数识别是阵列侧向测井的主要任务(Pan et al.，2014;Peng et al. 2001)，其中参数包括几何参数(如地层厚度、井眼半径、侵入带半径等)和物理参数(地层电导率).

实际测井过程中，一部分或全部监督电极Mi(1≤i≤6)作为测量电极，其上的电位可通过测井仪器直接读出，记为Uj(j=1，2，…，n). 我们希望通过这些测量值反演出各地层区域几何尺寸及其电导率，从而帮助判断地层中石油的储量.

采用最小平方误差准则，可将上述反演问题转化为如下多参数优化问题:

(10)

其中，n为测量电位数据的个数，Uj为第j次电位测量值，为其相应的计算值;电导率参数向量σ=(σ1，σ2，…，σp)T， 几何参数向量r=(r1，r2，…，rq)T.

实际测井中，各个参数的重要性不同，而且有些参数往往比较容易获取，如井眼半径、泥浆电导率等. 反演中，原始地层电阻率至关重要，因为它可用来确定地层孔隙度、 渗透率及含油气饱和度等重要地质参数.

3.2 信赖域算法

考虑如下非线性最小二乘问题:

(11)

其中，x∈Rn(m≥n)，F(x)∶Rn→Rm是连续可微函数. ‖·‖表示向量的2范数.

信赖域算法基本思想: 在当前迭代点xk的附近用一简单函数近似目标函数. 利用近似函数在某一邻域内的极小值点作为下一次迭代点. 对当前迭代点xk，设定xk+1=xk+s.信赖域算法通过求解如下信赖域子问题得到修正步长s:

(12)

其中J(xk)∈Rm×n为当前点xk处的Jacobi矩阵，hk为信赖域半径. 不难求解优化问题的步长：

(13)

其中非负迭代参数

(14)

通过引进非负参数λk，克服了Jacobi矩阵J(xk)几乎奇异或坏条件时牛顿步带来的困难. 事实上，方程(13)即传统的LM步. 因此，从这个意义下讲LM法可看作一类特殊的信赖域方法.

值得指出的是，LM法每次迭代修改λk，从而起到隐式调节hk和限制步长‖sk‖2的作用;而信赖域方法直接显式调节hk，故能更好地控制‖sk‖2的大小，而且方法直观容易理解以及具有很强的逼近背景.

信赖域半径hk的选取是该算法的关键. 为了确定信赖域半径，记Δfk为第k步实际下降量，即

(15)

记Δqk为对应的预测下降量，即

(16)

定义比值

(17)

注意到rk从某种程度上反映了二次函数qk(sk)与原目标函数f(xk+sk)的逼近程度，故可按如下方式调整信赖域半径. 若rk比较小接近零时，说明近似程度差，缩小下一步信赖域半径hk+1;若rk较大接近1且信赖域子问题的解在边界‖s‖2=hk上达到，说明近似程度好，而hk偏小了，下一步的hk+1应放大.

求解优化问题的信赖域算法如下:

步骤1 (初始化)给定初始点x1，初始半径h1，精度要求ε，参数0≤p0＜p1＜p2＜1，并令k=1.

步骤2 (收敛性检测)若‖F(xk)‖2＜ε或‖∇Fj(xk)‖2＜ε(j=1，2，…，m)，算法终止，得问题的近似解xk，否则转步骤3.

步骤3 (子问题求解)求解信赖域子问题得到建议步长sk.

步骤4 (信赖域修正)若rk＜p1，令hk+1=hk/2;若rk ＞p2及‖sk‖2=hk，令hk+1=2hk;否则令hk+1=hk.

步骤5 (可接受检测)若rk≤p0，则令xk+1=xk;否则令xk+1=xk+sk，k=k+1， 转步骤2.

本文计算中取p0=0，p1=0.25，p2=0.75. 信赖域算法迭代计算过程中，耗时部分主要集中在Jacobi矩阵的计算.

3.3 目标函数梯度的计算

采用信赖域算法求解优化问题的过程中，需要计算目标函数关于反演参数的梯度

(18)

直接对(10)求偏导数可得

(19)

将方程(6)代入上式，则有

(20)

类似地，可得

(21)

由此可知，如何计算基本解uk对参数的偏导数和将成为求解反演问题的关键.

利用有限元理论，容易证明基本解uk满足的椭圆型方程边值问题等价于如下变分问题(李大潜等，1980): 求u∈V使得

(22)

其中有限元空间

(23)

而表示函数u的梯度.

将式(22)对参数σi形式求导，可得满足如下变分问题:求vi∈V使得

(24)

这种形式求导是合理的(Cai，1998)，因此求解变分问题即可得到偏导数ukσi.

基本解uk对几何参数ri的偏导函数，难以得到其解析表达式或满足的变分问题，可利用如下中心差商来逼近:

其中，步长ε为较小的正常数，ei为第i个单位向量，仅第i个分量为1，其余分量均为0.

4 数值结果

所有计算均在个人笔记本上完成，CPU为Intel(R)Core(TM)i3-3120(2.5 GHz)，内存为4 GB. 采用Fortran语言编制了阵列侧向测井的多参数反演程序，编译环境为Intel Parallel Studio XE 2011，该开发工具套件包含Intel Visual Fortran Compiler XE 12.1和Intel® Math Kernel Library 10.3.

4.1 正演程序验证

反演的速度与精度严重依赖于正演. 为了验证正演程序的正确性，构造与已有文献(潘克家等，2013)中完全相同的正演模型: 井径rh=0.1016 m，侵入带半径ri=0.4016 m;井眼泥浆电阻率Rm=1 Ωm，围岩电阻率Rs=10 Ωm，侵入带电阻率Rx0=50 Ωm，原始地层电阻率Rt=500 Ωm;储层厚度H=2 m.

采用与文献(潘克家等，2013)中相同的网格剖分，电极附近r，z两个方向采用非常密的均匀剖分，远离电极采用对数均匀剖分，最终得到24564个矩形单元，24881个网格节点. 利用本文方法与文献(潘克家等，2013)中的方法分别进行计算. 两种方法得到的测井曲线几乎完全重合，屏蔽电极上的电流值Ij(k)(k，j=1，2，…，6)的最大相对偏差小于0.01%，验证了本文方法及正演程序的正确性.

计算时间方面，本文采用的PARDISO求解器正演一次(计算6条测井曲线)仅耗时0.45 s，残差2范数小于10-13，程序总耗时为0.85 s. 而JCG迭代算法要求残差2范数小于10-8，有限元方程求解耗时6.30 s. 因此，PARDISO求解器不仅计算精度高，且对此类问题速度比JCG法快十多倍，有利于快速反演.这主要是因为PARDISO为直接法，对7次不同源位置的基本解的求解，刚度矩阵只须进行一次分解;而迭代法JCG对7个基本解必须分别独立计算.

4.2 物理四参数反演

本节考虑物理参数即电阻率的反演. 首先利用正演算法对4.1节的理论模型进行计算，得到阵列侧向测井仪主电极A0上的6次测量的电位值. 然后以此6个不同的电位值作为测量数据，利用信赖域算法同时反演泥浆电阻率Rm、围岩电阻率Rs、侵入带电阻率Rx0和原始地层电阻率Rt四个物理参数.

表 1给出了不同初值下的电阻率四参数反演结果、迭代步数、迭代时间和最终目标函数值. 从表中可以看出，从4个不同初值出发，信赖域算法都能较快收敛到电阻率的真值，反演结果的最大相对误差达到10-6的量级. 当初值比较接近真值时，信赖域算法只须迭代很少几步，就能迅速收敛到原问题的真值. 即使初值非常差，如表中初值4(Rm=Rs=Rx0=Rt=100 Ωm)，信赖域算法仍能收敛到问题的准确解.

表 1 不同初值下的电阻率反演结果 Table 1 Resistivity inversion results for different initial values
图 3 信赖域算法目标函数收敛曲线 Fig. 3 Objective function convergence curves of the trust region algorithm

图 3给出了初值1、初值2和初值3情形下，信赖域算法的目标函数迭代收敛曲线. 从图中可看出: 1)信赖域算法对此类反演问题的收敛速度非常快，仅迭代几次，目标函数就能下降十几个数量级;2)初值离真值越近，信赖域算法收敛速度越快. 对初值1，信赖域算法仅迭代4次，目标函数值就从104下降到10-11.

下面考察信赖域算法反演阵列侧向测井数据的抗噪能力，为此向测量数据中添加一定信噪比(SNR)的高斯白噪声. 信噪比SNR定义为(潘克家等，2008)

其中，U和$\tilde{U}$分别为测井仪主电极A0上的真实电位值和添加了噪声的电阻值.

在上例中加入SNR分别为10、20、40和80 dB的高斯白噪声，以较差初值4作为初始猜测值，利用信赖域算法得到的反演结果、迭代步数、最终目标函数值和最大相对误差如表 2所示. 从表中可看出: 1)算法迭代次数(代表反演时间)与不加噪声时表 1中初值4的迭代次数50差不多;2)随着SNR的减小，噪声水平增大，反演结果变差;3)信赖域反演算法抗噪声的能力比较强，即使添加了20 dB的高斯白噪声，仍能得到比较满意的反演结果，最大相对误差为2.10%. 若更改初值进行反演，得到的反演结果、目标函数值和相对误差都与表格2中几乎相同，只有迭代步数有所改变. 这进一步验证了信赖域算法反演结果对初值的不敏感性，为一种全局收敛算法.

表 2 不同SNR下的电阻率反演结果 Table 2 Resistivity inversion results under different noise levels
4.3 传统三参数反演

考虑传统意义上的三参数反演，即同时反演侵入带半径ri、侵入带电阻率Rx0和原始地层电阻率Rt. 这三个参数也是实际测井资料处理中非常重要的参数. 人们可通过地层的电阻率来确定地层孔隙度、渗透率及含油气饱和度等重要地质参数. 而像泥浆电阻率Rm等参数比较容易获取，往往不需要通过反演获得. 同样以阵列侧向测井仪主电极A0上6次测量的不同电位值，利用信赖域算法同时反演ri、Rx0和Rt三个参数.

表 3给出了不同初值下传统三参数的反演结果、迭代步数、迭代时间和最终目标函数值. 从表中可看出，对不同的初值，信赖域算法都能较快收敛到参数的真值，反演结果的最大相对误差达到10-5的量级. 与4.2节电阻率四参数反演类似，传统三参数反演时初值接近真值时，信赖域算法只须迭代几步，就能迅速收敛到问题的真值. 即便初值不好，如表 3中初值4(ri=0.2 m，Rx0=Rt=100 Ωm)，信赖域算法仍能收敛到问题的真解. 这进一步验证了信赖域算法为区域搜索算法，具有全局收敛性，以及很强的初值适应能力.

表 3 不同初值下的三参数反演结果 Table 3 Inversion results of three parameters for different initial values

在上例反演数据中加入SNR分别为5、10、20 dB和40 dB的高斯白噪声，以ri=0.3 m，Rx0=40 Ωm，Rt=1000 Ωm作为反演初值，利用信赖域算法得到的传统三参数反演结果、迭代步数、最终目标函数值和最大相对误差如表 4所示. 从表 4中同样可得出类似结论: 加噪声时迭代次数与不加噪声时迭代次数相当;对传统的三参数反演，信赖域反演算法抗噪声的能力非常强，测量数据即使添加5dB的高斯白噪声，仍能得到非常满意的反演结果，最大相对误差小于4%.

表 4 不同SNR下的三参数反演结果 Table 4 Inversion results of three parameters under different noise levels
5 结论与展望

本文研究轴对称地层条件下高分辨率阵列侧向测井的多参数反演. 首先对阵列侧向测井的有限元正演方案进行改进，使之更适合反演计算. 建立非线性反演模型，在此基础上研究了目标函数对参数偏导数的计算，提出同时反演物理参数和几何参数的信赖域反演算法. 具体结论如下:

(1) 线性叠加原理有效地避开阵列侧向测井正演时屏蔽电极上电流的不确定性，且采用PARDISO直接求解器，对6次不同探测深度的测井模拟，有限元刚度矩阵只须分解一次，大大加快了HRLA的正演速度.

(2) 阵列侧向测井仪可测量6条具有不同探测深度的测井响应，利用主电极上的6个电位值可同时反演多个地层参数.

(3) 信赖域算法收敛速度快、具有全局收敛性，在不添加噪声的前提下，对较大范围内的初值，电阻率四参数反演和传统三参数反演均能得到非常准确的结果.

(4) 信赖域算法抗噪能力比较强，对测量数据添加SNR为10 dB(甚至5 dB)的高斯白噪声，仍能通过反演得到较为满意的地层参数.

本文基于阶跃式侵入模型，利用阵列侧向测井仪主电极A0上的6个电位值反演4个地层参数. 若同时利用监督电极M1，M2，…，M6或屏蔽电极A1，A2，…，A6上的电位测量值，有望反演更加复杂的泥浆侵入模型，甚至得到目的层特定区域的精细电导率成像.

致谢

本文完成于第一作者在美国北卡罗来纳大学夏洛特分校数学与统计系访学期间，作者衷心感谢国家留学基金委的资助. 同时感谢复旦大学“电法测井的数学建模及数值模拟”课题组李大潜院士和谭永基教授的有益指导.

参考文献
Cai Z J. 1998. Multiple parameters identification problems in resistivity well-logging. Chinese Annals of Mathematics , 19(3): 265–272.
Cao J J, Wang Y F. 2014. Seismic data restoration with a fast L1 norm trust region method. J. Geophys. Eng. , 11(4): 045010.
Chen Y H, Chew W C, Zhang G J. 1998. A novel array laterolog method. The Log Analyst , 39(5): 23–30.
Chew W C, Nie Z, Liu Q H, et al. 1991. An efficient solution for the response of electrical well logging tools in a complex environment. IEEE Transactions on Geoscience and Remote Sensing , 29(2): 308–313.
Deng S G, Fan Y R, Li Z Q, et al. 2011. Response characteristics of array lateral logs and their primary inversion in reservoirs with fracture-induced anisotropy. Petroleum Science , 8(1): 11–16.
Deng S G, Li Z Q. 2009. Simulation of array Laterolog response of fracture in fractured reservoir. Earth Science-Journal of China University of Geosciences (in Chinese) , 34(5): 841–847.
Deng S G, Li Z Q, Chen H. 2010a. The simulation and analysis of array lateral log response of fracture in coalbed methane reservoir. Coal Geology & Exploration (in Chinese) , 38(3): 55–60.
Deng S G, Xu R W, Jiang J L. 2010b. Study on the array Laterolog response of heterogeneous formation in deviated well. Well Logging Technology (in Chinese) , 34(2): 130–134.
Dun Y Q, Yuan J S. 2009. Fast inversion for array lateral electric-logging. Journal of Tsinghua University (Science and Technology) (in Chinese) , 49(11): 1871–1875.
Dun Y Q, Tang Z H, Zhang F, et al. 2010. Fast calculation and characteristic analysis of array lateral-logging responses. International Journal of Applied Electromagnetics and Mechanics , 33: 145–151.
Fan G X, Liu Q H, Blanchard S P. 2000. 3-D numerical mode-matching (NMM) method for resistivity well-logging tools. IEEE Transactions on Antennas and Propagation , 48(10): 1544–1552.
Fan Y R, Jiang J L, Deng S G, et al. 2009. Numerical simulation of high resolution array lateral logging responses. Well Logging Technology (in Chinese) , 33(4): 333–336.
Gould N I M, Scott J A, Hu Y F. 2007. A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations. ACM Transactions on Mathematical Software (TOMS) , 33(2): Article No. 10.
Hu H L, Xiao X, Pan K J, et al. 2014. Finite element modeling of 2. 5D DC resistivity based on locally refined graded mesh. Journal of Central South University (Science and Technology) (in Chinese) , 45(7): 2259–2267.
Li D Q, Zheng S M, Tan Y J, et al. Application of Finite Element Method in Electric Well Logging . (in Chinese) Beijing: Petroleum Industry Press, 1980 .
Li Z H, Wang Y F. 2014. A subspace trust-region method for seismic migration inversion. Optimization Methods and Software , 29(2): 286–296.
Li Z Q, Fan Y R, Deng S G, et al. 2010. Inversion of array laterolog by improved difference evolution. Journal of Jilin University(Earth Science Edition) (in Chinese) , 40(5): 1199–1204.
Lin T T, Hui F, Jiang C D, et al. 2013. Layered multi-exponential inversion method on surface magnetic resonance sounding dataset. Chinese J. Geophys. , 56(8): 2849–2861. doi: 10.6038/cjg20130833.
Liu Z H, Hu Q. 2002. Calculation and characteristics of array Laterlog responses. Journal of Xi'an Petroleum Institute(Natural Science Edition) (in Chinese) , 17(1): 53–57.
Liu Z H, Zhang X. 2005. Multi-parameter inversion of array laterolog responses. Journal of Xi'an Shiyou University(Natural Science Edition) (in Chinese) , 20(1): 30–33.
Nie Z P, Chen S Y. 1994. The numerical analysis for the response of dual lateral logging tool in a complex environment. Acta Electronica Sinica (in Chinese) , 22(6): 30–38.
Pan K J, Chen H, Tan Y J. 2008. Multi-exponential inversion of T2 spectrum in NMR based on differential evolution algorithm. Acta Physica Sinica (in Chinese) , 57(9): 5956–5961.
Pan K J, Tan Y J. 2009. High-efficient numeric simulation of spontaneous potential log in complex beds. Oil Geophysical Prospecting (in Chinese) , 44(3): 371–376.
Pan K J, Tan Y J, Hu H L. 2009. Mathematical model and numerical method for spontaneous potential log in heterogeneous formations. Applied Mathematics and Mechanics , 30(2): 209–219.
Pan K J, Tang J T, Hu H L, et al. 2012. Extrapolation cascadic multigrid method for 2. 5D direct current resistivity modeling. Chinese J. Geophys. , 55(8): 2769–2778. doi: 10.6038/j.issn.0001-5733.2012.08.028.
Pan K J, Wang W J, Tang J T, et al. 2013. Mathematical model and fast finite element modeling of high resolution array lateral logging. Chinese J. Geophys. , 56(9): 3197–3211. doi: 10.6038/cjg20130932.
Pan K J, Liu J L. 2014. A parameter identification problem for spontaneous potential logging in heterogeneous formation. Journal of Inverse and Ill-posed Problems , 22(3): 357–373.
Peng Y J, Tan Y J. 2001. Multi-parameter identification and applications in well-logging. Computational Geosciences , 5(4): 331–343.
Shen J S. 2004. Modeling of the multi-component induction log in anisotropic medium by using finite difference method. Progress in Geophysics (in Chinese) , 19(1): 101–107.
Su B Y, Fujimitsu Y, Xu J L, et al. 2012. A model study of residual oil distribution jointly using crosswell and borehole-surface electric potential methods. Applied Geophysics , 9(1): 19–26.
Tan M J, Zhang G J, Yun H Y, et al. 2007. 3-D numerical mode-matching (NMM) method for resistivity logging responses in nonsymmetric conditions. Chinese J. Geophys. (in Chinese) , 50(3): 939–945.
Tian Y, Chen X F. 2006. Simultaneous inversion of hypocenters and velocity using the quasi-Newton method and trust region method. Chinese J. Geophys. (in Chinese) , 49(3): 845–854.
Wang G L, Zhang G J, Cui F X, et al. 2003. Application of staggered grid finite difference method to the computation of 3-D induction logging response. Chinese J. Geophys. (in Chinese) , 46(4): 561–567.
Wang H N, Tao H G, Yao J J, et al. 2008. Study on the response of a multicomponent induction logging tool in deviated and layered anisotropic formations by using numerical mode matching method. Chinese J. Geophys. (in Chinese) , 51(5): 1591–1599.
Wang Y F, Cao J J, Yang C C. 2011. Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling. Geophys. J. Int. , 187(1): 199–213.
Weller A, Furche M, Sch n J H. 2003. Detection of layer boundaries in wells using multi-electrode resistivity data. Geophys. J. Int. , 153(1): 175–186.
Wu J, Xie W W, Xie X C, et al. 2008. Forward response analysis of array lateral logging tool. Journal of Xi'an Shiyou University (Natural Science Edition) (in Chinese) , 23(1): 73–76.
Yang S W, Wang H N, Chen G B, et al. 2009. The 3-D finite difference time domain (FDTD) algorithm of response of multi-component electromagnetic well logging tool in a deviated and layered anisotropic formation. Chinese J. Geophys. (in Chinese) , 52(3): 833–841.
Yang W. 2003. Forward and inversion of array laterolog logging data. Progress in Exploration Geophysics (in Chinese) , 26(4): 305–308.
Yuan N, Nie Z P, Nie X C. 1998. Numerical analysis of SP log response in complex heterogeneous media. Acta Geophysica Sinica (in Chinese) , 41(S1): 429–436.
Yuan Y X, Sun W Y. Optimization Theories and Methods . (in Chinese) Beijing: Science Press, 1997 .
Zeng Q, Shi X M, Wu Y S, et al. 2011. A study of one-dimensional magnetotelluric sounding inversion using the trust region algorithm. Progress in Geophys. , 26(3): 885–893. doi: 10.3969/j.issn.1004-2903.2011.03.013.
Zhang G J. Electrical Well Logging . (in Chinese) Beijing: Oil Industry Press, 1984 .
Zhang G J, Wang H M. 1996. Solution of the normal resistivity logging with the numerical mode-matching method. Journal of China University of Petroleum (Edition of Natural Sciences) (in Chinese) , 20(2): 23–29.
Zhang J H, Liu Z H, Wu J. Theory and Application of Electrical Well Logging . (in Chinese) Xi'an: Northwest University Press, 2002 .
Zhang Z Y, Zhou Z Q. 2002. Real-time quasi-2-D inversion of array resistivity logging data using neural network. Geophysics , 67(2): 517–524.
Zhu P, Lin C Y, Li Z Q, et al. 2015. Numerical simulation of array laterolog response in horizontal and highly deviated wells. Journal of Jilin University (Earth Science Edition) (in Chinese) , 45(6): 1862–1869.
邓少贵, 李智强. 2009. 裂缝性储层裂缝的阵列侧向测井响应数值模拟. 地球科学(中国地质大学学报) , 34(5): 841–847.
邓少贵, 李智强, 陈华. 2010a. 煤层气储层裂隙阵列侧向测井响应数值模拟与分析. 煤田地质与勘探 , 38(3): 55–60.
邓少贵, 徐悦伟, 蒋建亮. 2010b. 倾斜井非均匀地层的阵列侧向测井响应研究. 测井技术 , 34(2): 130–134.
顿月芹, 袁建生. 2009. 阵列侧向电法测井的快速反演. 清华大学学报(自然科学版) , 49(11): 1871–1875.
范宜仁, 蒋建亮, 邓少贵, 等. 2009. 高分辨率阵列侧向测井响应数值模拟. 测井技术 , 33(4): 333–336.
李大潜, 郑宋穆, 谭永基, 等. 有限元素法在电法测井中的应用. 北京: 石油工业出版社, 1980 .
李智强, 范宜仁, 邓少贵, 等. 2010. 基于改进差分进化算法的阵列 侧向测井反演. 吉林大学学报:地球科学版 , 40(5): 1199–1204.
林婷婷, 慧芳, 蒋川, 等. 2013. 分层多指数磁共振弛豫信号反演方法研究. 地球物理学报 , 56(8): 2849–2861.
刘振华, 胡启. 2002. 阵列侧向测井响应的计算及其特征. 西安石油学院学报(自然科学版) , 17(1): 53–57.
刘振华, 张霞. 2005. 阵列侧向测井响应的多参数反演. 西安石油大学学报(自然科学版) , 20(1): 30–33.
聂在平, 陈思渊. 1994. 复杂介质环境中双侧向测井响应的高效数值分析. 电子学报 , 22(6): 30–38.
潘克家, 陈华, 谭永基. 2008. 基于差分进化算法的核磁共振T2谱多指数反演. 物理学报 , 57(9): 5956–5961.
潘克家, 谭永基. 2009. 复杂地层中自然电位测井的高效数值模拟. 石油地球物理勘探 , 44(3): 371–376.
潘克家, 汤井田, 胡宏伶, 等. 2012. 直流电阻率法2. 5维正演的外推瀑布式多重网格法. 地球物理学报 , 55(8): 2769–2778.
潘克家, 王文娟, 汤井田, 等. 2013. 高分辨率阵列侧向测井的数学模型及有限元快速正演. 地球物理学报 , 56(9): 3197–3211.
沈金松. 2004. 用有限差分法计算各向异性介质中多分量感应测井的响应. 地球物理学进展 , 19(1): 101–107.
谭茂金, 张庚骥, 运华云, 等. 2007. 非轴对称条件下用三维模式匹配法计算电阻率测井响应. 地球物理学报 , 50(3): 939–945.
田玥, 陈晓非. 2006. 利用拟牛顿法和信赖域法联合反演震中分布与一维速度结构. 地球物理学报 , 49(3): 845–854.
汪功礼, 张庚骥, 崔锋修, 等. 2003. 三维感应测井响应计算的交错网格有限差分法. 地球物理学报 , 46(4): 561–567.
汪宏年, 陶宏根, 姚敬金, 等. 2008. 用模式匹配算法研究层状各向异性倾斜地层中多分量感应测井响应. 地球物理学报 , 51(5): 1591–1599.
仵杰, 谢尉尉, 解茜草, 等. 2008. 阵列侧向测井仪器的正演响应分析. 西安石油大学学报(自然科学版) , 23(1): 73–76.
杨守文, 汪宏年, 陈桂波, 等. 2009. 倾斜各向异性地层中多分量电磁波测井响应三维时域有限差分(FDTD)算法. 地球物理学报 , 52(3): 833–841.
杨韡. 2003. 阵列侧向测井的正反演. 勘探地球物理进展 , 26(4): 305–308.
袁宁, 聂在平, 聂小春. 1998. 自然电位测井响应的高效数值模拟. 地球物理学报 , 41(S1): 429–436.
袁亚湘, 孙文瑜. 最优化理论与方法. 北京: 科学出版社, 1997 .
曾奇, 师学明, 武永胜, 等. 2011. 一维大地电磁信赖域反演法研究. 地球物理学进展 , 26(3): 885–893.
张庚骥. 电法测井. 北京:: 石油工业出版社, 1984 .
张庚骥, 汪涵明. 1996. 普通电阻率测井的数值模式匹配解法. 石油大学学报(自然科学版) , 20(2): 23–29.
张建华, 刘振华, 仵杰. 电法测井原理与应用. 西安: 西北大学出版社, 2002 .
祝鹏, 林承焰, 李智强, 等. 2015. 水平井和大斜度井中阵列侧向测井响应数值模拟. 吉林大学学报(地球科学版) , 45(6): 1862–1869.