用多分量感应资料快速重建层状地层的纵横向电阻率和水平界面深度

汪宏年，商庆龙，朱天竹，李舟波

1.吉林大学物理学院，长春 130012 2.吉林大学地球探测科学与技术学院，长春 130026

用多分量感应资料快速重建层状地层的纵横向电阻率和水平界面深度

汪宏年1，商庆龙1，朱天竹1，李舟波2

1.吉林大学物理学院，长春 130012 2.吉林大学地球探测科学与技术学院，长春 130026

研究建立用多分量感应测井资料同时快速重建水平层状横向同性介质中横向与纵向电阻率和层界面深度的有效方法。首先，利用电磁场摄动方程、电导率函数与模型参数关系方程以及模式匹配算法得到电磁场并矢Green函数的半解析解，建立多分量感应测井响应的Frèchet导数矩阵的快速算法；在此基础上，借助于规范化处理和奇异值分解技术，给出同时反演水平层状地层中各个地层的纵、横向电阻率以及层界面深度的迭代过程，实现理论合成资料与输入资料的最佳拟合。数值计算证明，该反演算法能够取得较满意的反演效果。

感应测井；Frèchet导数矩阵快速算法；各向异性地层电阻率和层界面的同时反演

0 引言

多分量感应测井作为能够同时测量地层纵横向电阻率的一种新测井技术，在薄交互层、裂缝等各向异性地层评价中发挥着重要作用［1］。与各种常规测井仪器不同，多分量感应测井仪器增加了与井轴方向完全垂直的发射和接收线圈系（共面线圈系），其测量信号能够提供地层纵向电阻率信息；但由于多分量感应测井仪器产生的相应的感应电磁场空间分布不再具有轴对称性，且电场方向与地层界面也不再平行，这样会在地层界面上产生积累面电荷而导致共面线圈系的测井响应更加复杂［2－7］①张烨，汪宏年，陶宏根，等.基于耦合标势与矢势的有限体积法模拟非均匀各向异性地层中多分量感应测井三维响应.地球物理学报，待刊.。数值计算结果显示，共面线圈系的测量结果受井眼和层界面的影响均很大，特别在薄层上或层界面附近往往出现负响应。根据共面线圈系的视电导率σa，xx和σa，yy的大小无法直接估计出地层纵向电阻率的真实变化情况，必须建立一套有效的数据处理技术，以便从测量结果中提取出地层电阻率和层界面。根据笔者目前掌握的资料，多分量感应资料的处理方法主要有2种：一种是Michael等［1］提出的最小二乘法提取纵横向电阻率方法，该方法直接利用均匀介质中电磁场的解析表达式提取地层电阻率参数，由于没有考虑层界面影响，其处理结果只是消除了地层的趋肤效应，在地层厚度较小的交互层上或层界面附近将难以取得满意效果；另一种方法是Zhang等［8］借助迭代反演同时提取层界面和地层电阻率的算法，该算法中没有建立计算目标函数下降方向的快速有效方法，所以其效率往往很低。对于层状地层来说，如果不考虑井眼影响，多分量感应数值模拟仍然属于2.5维问题［2，4］，通过Fourier变化方法可将其转化成3个轴对称向量场的求解问题，借助模式匹配方法可以实现快速正演计算。如果将这种快速数值模拟技术与摄动理论结合起来，建立一套有效的Frèchet导数的快速算法［9－10］，就能够实现多分量感应资料的快速迭代反演。

笔者将研究并建立这种用微分方程反演理论建立起来的多分量感应资料的快速反演方法，将包括摄动方程与Frèchet导数矩阵的快速计算，以及纵横向电阻率和层界面的迭代反演等内容，最后通过数值模拟结果检验在水平层状垂直井眼中的反演结果。

1 理论

为简单起见，用图1表示水平层状各向异性地层。地层由N＋1个各向异性地层组成，dn（n＝2，3，…，N）表示地层界面深度，每个地层的电导率珔σn＝diag（σhn，σhn，σvn）是对角张量，n＝1，2，…，N＋1其中σhn和σvn分别表示地层横向和纵向电导率。对于这种层状地层模型，其电导率空间分布完全由模型向量x＝（σh1，σh2，…，σhn，σv1，σv2，…，σvn，d2，d3，…，dn）确定。建立如图1所示的坐标系xyz，z轴与地层层面垂直，x轴和y轴与地层层面平行。Mp，p＝x，y，z分别是x，y，z正交磁偶极子产生的磁偶极矩，Rp，p＝x，y，z，分别是x，y，z3个不同方向正交接收线圈位置，如图1所示，两者的源距假定为L。p方向的磁偶极子Mp产生的电场Ep和磁场Hp完全由Maxwell方程确定［2－3］：

其中：Mp＝＾epMpδ（r－rs），＾ep为坐标p上的单位向量；r和rs分别为场点和源点位置。磁导率μ（H／m）为常数；ω（rad／s）为发射信号角频率；因发射频率较低故忽略位移电流影响。

图1 地层模型与坐标系Fig.1 Formation model and coordinate system

1.1 摄动方程与Frèchet导数矩阵的快速计算

利用摄动理论，从方程（1）中可直接推导出电导率微小摄动δ珔σ与电磁场变化关系：

引入磁流源电场并矢Green函数GME（r，rR），磁流源磁场并矢Green函数GMH（r，rR）：

则利用文献［10－11］中的结果，从方程（2）和（3）中可以得到磁场变化δHp的积分解：

其中：rR是接收点位置；上角标“T”表示转置。

对于图1所示的水平层状地层，电导率是分片常数函数，电导率与模型向量x的各个分量间满足如下方程［12］：

其中：δ珔σn＝diag（δσhn，δσhn，δσvn）和δdn分别表示n层上电导率向量和层界面n的微小变化；δ（z）＝H′（z）是Dirac’sδ函数。

因为发射源位于井轴上，利用模式匹配算法［2－3］，在进行磁场正演计算的同时可得到p方向的发射线圈在地层n上的水平和垂直电场强度各个分量的半解析解：

方程（3）中磁流源并矢Green函数珚GME（r，rR）可用同样的方法得到，其3个列向量gq（r；rR），q＝x，y，z3个不同方向正交发射线圈，具有与方程（9）类似的表达式：

由于接收点rR可能不在井轴上，方程（10）右端是无限多个调和分量的和，将方程（9）和（10）代入方程（4）后，得

此外，将方程（5）代入方程（11）并经简单整理，则方程（11）右端项中的被积函数可展开成如下表达式：

所以方程（11）右端的整个积分可简化为

通过方程（7）和（8）得到的半解析式，方程（13）的右端项中关于z的积分可以解析得到，而关于ρ的积分可数值计算且积分与源点和接收点的位置无关，所以可事先计算并储存起来以提高积分效率［9，11］。一旦计算出方程（13）中的各个积分，并代入方程（11）中，就能够得到磁场变化与模型向量摄动方程：

其中：p和q分别表示发射线圈和接收线圈的方向；Fqp（rR；rs）是磁场Hqp（rR；rs）相对于模型向量x的Frèchet导数。取方程（14）中3个主分量的虚部后再除以相应的仪器常数，就得到了多分量感应测井3个主视电导率变化与模型向量摄动方程：

其中，Cpp（zR）是记录点zR处的多分量感应测井视电导率对模型向量的Fréchet导数。由于模型参数δx中包含有电导率和层界面2个不同物理量纲的物理量，为便于反演计算，采用规范化处理技术将所有变量转化为无量纲的纯数：

其中：δq＝（δqσ，δqh）T是模型向量的相对变化量，且δqσ＝（δσh，1／σh，1，…，δσh，N／σh，N，δσv，1／σv，1，…，δσv，N／ σv，N），δqh＝（δh2／（h3－h2），δh3／（h4－h2）／2，…，δhN－1／（hN－hN－2）／2，δhN／（hN－hN－1））；而δypp（zR）＝δσa，pp（zR）／σa，pp（zR）是视电导率的相对变化；Fpp（zR）＝σ－1a，pp（zR）Cpp（zR）Q是规范化Fréchet导数，且

1.2 迭代反演

对于水平层状各向异性地层，在井眼垂直情况下，视电导率σa，xx和σa，yy完全相同且它们与地层的纵横向电阻率有关，而视电导率σa，zz仅反映地层水平电导率变化。所以从多分量感应测井资料中提取地层电阻率可以采用2种不同的处理方法：一种是先对σa，zz进行迭代反演确定出水平电导率，然后再进一步用σa，xx反演地层纵向电导率；另一种方法是将视电导率σa，xx和σa，zz结合在一起，同时反演纵横向电导率和层界面。对于前一种反演方法，由于水平电导率的反演误差会对纵向电导率的反演结果造成较大影响，所以笔者采用后种方法。

设Y0＝（）T是多分量感应测井仪器测量出的一系列视电导率σa，xx和σa，zz组成的向量，反演目的就是找出某个模型向量x0，使其理论模拟记录Yc能够最佳拟合于Y0。由于测井记录是模型向量x的非线性函数，利用线性化方程（16）可以得到反演结果的迭代解［12］：

其中：x（0）是模型向量的初值；x（k）是第k次迭代反演结果；“I”是单位矩阵；“k”是预定的迭代总次数，实际数值实验表明，通常情况下迭代8次就足够了。而向量δq（k）是方程（16）的广义解：

其中：Δy（k）＝diag（Y0－1）（Y0－Y（k））是理论测井记录与输入资料之间的相对误差；F（k）＝U∑VT＝V∑－1UT是（16）式中规范化处理后的Fréchet导数矩阵的SVD分解。

2 数值结果

为检验上述反演方法的有效性，对理论模型数值模拟结果进行了反演。在反演过程中，仪器的源距L＝1m，源的发射频率f＝20kHz，即仅考虑单发单收这种最简单情形。

图2 无噪声多分量感应测井资料的反演结果Fig.2 The inversion results from noiseless multicomponent induction logging data from 15layer formation

图3 含1%噪声多分量感应测井资料的反演结果Fig.3 The inversion results from 1%noise multicomponent induction logging data from 15layer formation

模型由15个各向异性层状地层组成，从上到下各个地层的纵向电阻率分别是：3、30、4、50、80、15、40、25、110、20、50、120、35、10、6Ω·m，横向电阻率：1、5、2、10、20、5、15、8、25、6、12、30、8、4、2Ω·m。13个层厚：1.65、1.50、1.85、2.30、1.75、2.20、1.85、2.50、1.95、1.80、2.65、2.10、1.85m。对于该模型首先通过正演软件计算出视电导率曲线σa，xx和σa，zz，并作为反演计算的输入数据（采样间距0.05m）。反演前，假定各个地层纵横向电阻率的初值是其真值的1／3，而各个地层界面的初值由公式h（0）n＝hn＋（hn＋1－hn－1）确定，由此得到的地层厚度的初始误差约为20%。

图2给出了迭代8次以后的反演结果。不难看出，所有结果均吻合得很好。但纵向电阻率的反演误差比横向电阻率大，说明纵向电阻率较难反演；此外，与常规的感应测井反演结果不同，在多分量的反演结果中似乎低电阻层上的反演误差更大，说明纵向电导率大的地层其参数更难计算。

为了进一步检验反演方法的抗噪能力，图3给出了在输入数据中加入1%白噪声后的反演结果。图3a的结果显示，反演模型的合成数据与输入资料仍然吻合得较好；而图3b和3c表明，反演结果与地层真值也很接近，但与无噪声情况相比，反演误差有了明显增加，此外，纵向电阻率的反演误差也增加得更为明显。这里需要特别指出的是，反演结果的抗噪能力不仅与反演算法有关，还与模型参数密切相关：当地层厚度较大时，式（16）中磁场相对于该层参数的Frechet导数将比较大，这时其相应的抗噪能力也会更强。

3 结论

1）在薄交互层上，多分量感应测井资料中共面视电导率虽然存在负响应，但由于这种响应与地层参数间具有确定性关系，所以利用迭代反演算法仍然可以快速重构各向异性电阻率和层界面位置。

2）数值结果显示，地层横向电阻率比纵向电阻率更容易反演，其反演结果的抗噪能力也比纵向电阻率更强。

3）与常规的感应测井反演结果不同，在多分量反演结果中似乎低电阻层上的反演误差更大，说明纵向电导率大的地层其参数更难计算。

（References）：

［1］ Michael Z，David K，Ertan P.Foundation of the Tensor Induction Well Logging［J］.Perophysics，2001，42（6）：588－610.

［2］ 汪宏年，陶宏根，姚敬金，等.用模式匹配算法研究层状各向异性倾斜地层中多分量感应测井响应［J］.地球物理学报，2008，51（5）：1591－1599.

Wang Hong－nian，Tao Hong－gen，Yao Jing－jin，et al.Study on the Response of Multicomonent Induction Logging Tool in Deviated and Layered Anisotropic Formations by Using Numerical Mode Matching Method［J］.Chinese Journal of Geophysics，2008，51（5）：1591－1599.

［3］ 汪宏年，胡平，陶宏根.水平层状非均质横向同性地层中阵列多分量感应测井响应的快速计算［J］.地球物理学报，2012，55（2）：717－726.

Wang Hong－nian，Hu Ping，Tao Hong－gen，et al.Fast Algorithm of Responses of Array Multicomponent Induction Logging Tool in Horizontally Stratified Inhomogeneous TI Media［J］.Chinese Journal of Geophysics，2012，55（2）：717－726.

［4］ Wang Hong－nian，Tao Hong－gen，Yao Jing－jin.An Efficient and Reliable Simulation of Multicomponent Induction Logging Response in Horizontally Stratified Inhomogeneous TI Formations by Numerical Mode－Matching Method［J］.Geoscience and Remote Sensing，DOI：10.1109／TGRS.2012.2183135.

［5］ 姚东华，汪宏年，杨守文，等，用传播矩阵法研究层状正交各向异性地层中多分量感应测井响应［J］.地球物理学报，2010，53（12）：3026－3037.

Yao Dong－hua，Wang Hong－nian，Yang Shou－wen，et al.Study on the Responses of Multi－Component Induction Logging Tool in Layered Orthorhombic Anisotropy Formations by Using Propagator Matrix Method［J］.Chinese Journal of Geophysics，2010，53（12）：3026－3037.

［6］ 杨守文，汪宏年，陈桂波，等.倾斜各向异性地层中多分量高频电磁波测井响应的三维时域有限差分（FDTD）算法［J］.地球物理学报，2009，52（3）：833－841.

Yang Shou－wen，Wang Hong－nian，Chen Gui－bo，et al.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［J］.Chinese Journal of Geophysics，2009，52（3）：833－841.

［7］ Wang Hong－nian，So Poman，Yang Shou－wen，et al.Numerical Modeling of Multicomponent Induction Well Logging Tools in the Cylindrically Stratified Anisotropic Media［J］.Geoscience and Remote Sensings，2008，46（4）：1134－1146.

［8］ Zhang Zhiyi，Yu Liming，Mauser K B.Simultaneous Determination of Relative Angles and Anisotropic Resistivity Using Multicomponent Lnduction Logging Data［J］.Geophysics，2004，69（4）：898－908.

［9］ Wang Hong－nian，Tao Hong－gen，Yao Jing－jin，et al.Fast Multiparameter Reconstruction of Multicomponent Induction Well Logging Datum in Deviated Well in a Horizontally Stratified Anisotropic Formation［J］.Geoscience and Remote Sensings，2008，46（5）：1525－1534.

［10］ Wang Hong－nian.Adaptive Regularization Iterative Inversion of Array Multicomponent Induction Well Logging Datum in a Horizontally Stratified Inhomogeneous TI Formation［J］.Geoscience and Remote Sensing，2011，49（11）：4483－4492.

［11］ Wang Hong－nian.Simultaneous Reconstruction of Geometric Parameters and Resistivity from the Multiarray Induction Log in a Horizontally Layered Formation［J］.Geoscience and Remote Sensing，2003，41（1）：185－194.

Simultaneously Fast Reconstruction of Resistivities and Interfaces in Horizontally Stratified TI Formation by Using Multicomponent Induction Well Logging Data

Wang Hong－nian1，Shang Qing－long1，Zhu Tian－zhu1，Li Zhou－bo2

1.College of Physics，Jilin University，Changchun 130012，China 2.College of GeoExploration Science ＆Technology，Jinlin University，Changchun 130026，China

We advance an efficient method to simultaneously reconstruct the horizontal and vertical resistivities and the interface depth each bed in the horizontally stratified transversely isotropic（TI）formation from the multicomponent induction well logging（MCIL）data.We first give the perturbed equation of electromagnetic（EM）fields，relation of conductivity function to the model parameters and the semi－analytic solutions of EM tensor Green function provided by the numerical mode matching method in the horizontally layered TI formation，and further set up a fast algorithm of Frèchet derivatives of multicomponent induction logging response with respect to all the formation resistivities and interfaces.Then，by combination of normalization and singular value decomposition（SVD）technique，we advance an iterative method to simultaneously reconstruct both horizontal and vertical resistivities and interface per bed so that we realize the best fit of the input data with the modeling logs.Numerical tests validate the algorithm.

induction logging；fast algorithm of Frèchet derivatives；reconstruction of resistivitiesand interface at same time

book=2012,ebook=678

P631.81

A

1671－5888（2012） 04－0900－06

2012－04－18

国家自然科学基金项目（40874058）

汪宏年（1962－），男，教授，博士生导师，主要从事非均质与各向异性介质中电磁场数值模拟与反演研究，E－mail：wanghn＠mail.jlu.edu.cn。