引用本文
王磊, 范宜仁, 袁超, 巫振观, 邓少贵, 赵伟娜. 随钻方位电磁波测井反演模型选取及适用性[J]. 石油勘探与开发, 2018,45(5): 914-922.
WANG Lei, FAN Yiren, YUAN Chao, WU Zhenguan, DENG Shaogui, ZHAO Weina. Selection criteria and feasibility of the inversion model for azimuthal electromagnetic logging while drilling (LWD)[J]. Petroleum Exploration & Development, 2018,45(5): 914-922.  
Doi: 10.11698/PED.2018.05.18
Permissions
Copyright©2018, 石油勘探与开发》编辑部
石油勘探与开发》编辑部 所有
随钻方位电磁波测井反演模型选取及适用性
王磊1,2,3, 范宜仁1,2,3, 袁超4, 巫振观1,2,3, 邓少贵1,2,3, 赵伟娜5
1. 中国石油大学（华东）地球科学与技术学院,山东青岛 266580
2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛 266071
3. 中国石油大学CNPC测井重点实验室,山东青岛 266580
4. 中国石油勘探开发研究院,北京 100083
5. 中国科学院力学研究所流固耦合系统力学重点实验室,北京 100190

联系作者简介：巫振观（1991-）,男,福建大田人,现为中国石油大学（华东）在读博士研究生,主要从事电法测井正反演方面研究。地址：山东省青岛市经济技术开发区长江西路66号,中国石油大学（华东）地球科学与技术学院,邮政编码：266580。E-mail: wuzg2014@163.com

第一作者简介：王磊（1989-）,男,山东潍坊人,现为中国石油大学（华东）在站博士后,主要从事电法测井正反演与电测井解释方面研究。地址：山东省青岛市经济技术开发区长江西路66号,中国石油大学（华东）地球科学与技术学院,邮政编码：266580。E-mail: upcwanglei199133@163.com

收稿日期: 2018-01-19
基金项目: 国家自然科学基金项目（41474100,41574118,41674131）; 中央高校基本科研业务费专项资金资助（17CX06041A）; 中国石油大学（华东）研究生创新工程（YCX2017004）
摘要

随钻方位电磁波测井可提供邻近地层界面的准确位置信息而被广泛运用于实时地质导向钻井,但如何选取合适的反演模型和最优化算法则是决定随钻方位电磁波测井资料反演速度和精度的关键。为此,首先基于降维策略将复杂地层随钻方位电磁波测井资料3D反演问题简化为一系列1D问题,然后探讨不同1D反演模型、反演算法的可行性与反演效果,并给出了两者组合的最优选取方法。数值模拟结果表明,1D反演模型的选取取决于靶层的厚度,而反演算法的选取则依赖于反演模型的层数,即对靶层厚度为4.0 m以上的厚层,应选用单界面反演模型和梯度反演算法组合;当靶层厚度为1.0~4.0 m时,仅需将单界面反演模型替换为双界面反演模型即可提供准确的邻近上、下地层界面;对薄层随钻方位电磁波资料,则需采用多界面反演模型和Bayesian反演算法。图8表1参14

关键词: 随钻方位; 电磁波测井; 反演模型; 反演算法; 地质导向; 地层界面
中图分类号:P631.8 文献标志码:A 文章编号:1000-0747(2018)05-0914-09
Selection criteria and feasibility of the inversion model for azimuthal electromagnetic logging while drilling (LWD)
WANG Lei1,2,3, FAN Yiren1,2,3, YUAN Chao4, WU Zhenguan1,2,3, DENG Shaogui1,2,3, ZHAO Weina5
1. School of Geosciences, China University of Petroleum, Qingdao 266580, China
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China
3. CNPC Key Laboratory for Well Logging, China University of Petroleum, Qingdao 266580, China
4. Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China;
5. Key Laboratory for Mechanics in Fluid Solid Coupling Systems, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China
Abstract

Azimuthal electromagnetic logging while drilling (LWD), which is capable of providing accurate position information approaching bed boundary, has been widely applied in real-time geosteering. For the inversion speed and precision of azimuthal electromagnetic LWD data, the key lies in the selection of proper inversion model and corresponding optimization algorithm. In this study, we first simplified the complex three-dimensional (3D) inversion of data into a series of one-dimensional (1D) inversion problems by using the dimensionality reduction scheme. Then, the feasibility and inversion performance of various 1D inversion models and different optimization methods were investigated, and the best combination between the inversion model and inversion algorithm was also obtained. Numerical simulation results show that the selection of 1D inversion model is dominated by the thickness of targeted beds, whereas the determination of inversion algorithm relies on the total layers amount of the inversion model. For the formation with thickness larger than 4.0 m, the single boundary inversion model and gradient optimization method are recommended. When the bed thickness is between 1.0 m and 4.0 m, the two-boundary inversion model instead of the single-boundary one is needed to estimate upper and lower boundaries around the borehole. For the inversion of azimuthal electromagnetic LWD data of thin layers, the multiple- boundary inversion model and the Bayesian algorithm should be employed.

Keyword: azimuth while drilling; electromagnetic logging; inversion model; inversion algorithm; geosteering; bed boundary
0 引言

地层界面位置的实时确定是随钻地质导向的关键之一, 其对准确钻井与油气产能最大化具有重要意义[1, 2]。虽然利用随钻方位电磁波测井响应特征可以定性判断界面位置, 但复杂的井下环境导致测井响应极为复杂[3], 井眼到地层界面的定量评价须通过资料的反演处理方可获得。

大斜度井/水平井电测井资料处理通常为复杂的三维（3D）问题, 反演时需考虑地层的倾斜度、各向异性及钻井液侵入等各种问题[4, 5, 6]。当忽略井眼大小、钻井液侵入且假定地层倾斜方向沿某一方向保持不变时, 电测井资料反演可简化为二维（2D）反演[7, 8]。但是, 电测井3D和2D反演总体仍不成熟, 这主要归因于：①井下电测井测量信息数量和种类少, 导致3D和2D反演结果存在一定不确定性; ②缺乏3D/2.5D正演快速算法; ③反演参数过多, Jacobian矩阵计算量大。目前, 工业界普遍采用基于滑动开窗的降维处理策略对随钻电磁波类测井资料进行快速处理。Yang等忽略地层横向非均质的影响, 基于简化的3层反演模型消除了视电阻率异常“ 犄角” [9]。利用同样的反演方法, Li和Zhou实现了邻近地层界面的实时提取[10]。因此, 将复杂的3D反演简化为1D问题, 可以大大简化反演算法和模型的复杂程度。

随钻方位电磁波资料1D处理的关键在于简化模型的选取和反演算法的确定。目前, 对随钻方位电磁波测井资料的处理仅限于双界面模型, 而双界面模型的适用性和单、多界面模型的可行性仍有待研究。一般而言, 反演算法的选取依赖于反演模型。梯度类算法具有迭代次数少、反演速度快等优点, 被广泛应用于基于双界面模型的随钻方位电磁波测井资料实时处理。但梯度算法的准确性取决于模型的初始值, 若初始模型选取不当, 则反演结果的准确性难以保证, 可能导致地质导向决策失误等问题[11]。相比于双界面模型, 多界面模型的代价函数极小值更多, 且待反演参数已知信息少, 梯度类算法不再适用。随机Bayesian反演算法具有全局收敛特性、无需计算Jacobian矩阵等优点, 近年来在电测井反演中的应用逐渐增多[12]。目前, 随钻方位电磁波测井资料处理一般采用Gauss- Newton算法, Bayesian反演算法的应用尚未见报道, 其收敛性与反演速度等也有待研究。

本文首先介绍了复杂地层模型的反演降维策略和正则化Gauss-Newton算法、Bayesian算法原理; 然后, 针对降维后的1D反演问题, 对比分析了不同反演模型、反演算法的适用性与反演效果; 最后, 分别利用优选后的反演模型与最优化算法对现场资料进行处理, 取得了良好的应用效果。

1 随钻方位电磁波测井原理与反演降维策略
1.1 随钻方位电磁波测井原理

本节以PeriScope仪器为例, 简要介绍随钻方位电磁波测井测量原理。PeriScope包括两类测量信息：视电阻率信息与地质信息（见图1a）, 前者通过采集轴向发射线圈T在轴向接收线圈R1和R2处激发的电势 $V_{R_{1}}$和$V_{R_{2}}$ , 并利用（1）式求出两者的幅度比（Att）和相位差（PS）。最后, 利用转换图版, 将相位差和幅度比转换为视相位差电阻率和视幅度比电阻率。

$Att=20\text{lg}\frac{\sqrt{{{\text{ }\!\![\!\!\text{ }\operatorname{Re}al\text{(}{{V}_{{{\text{R}}_{1}}}}\text{) }\!\!]\!\!\text{ }}^{2}}+{{\text{ }\!\![\!\!\text{ }\operatorname{Im}ag\text{(}{{V}_{{{\text{R}}_{1}}}}\text{) }\!\!]\!\!\text{ }}^{2}}}}{\sqrt{{{\text{ }\!\![\!\!\text{ }\operatorname{Re}al\text{(}{{V}_{{{\text{R}}_{2}}}}\text{) }\!\!]\!\!\text{ }}^{2}}+{{\text{ }\!\![\!\!\text{ }\operatorname{Im}ag\text{(}{{V}_{{{\text{R}}_{2}}}}\text{) }\!\!]\!\!\text{ }}^{2}}}}$

$PS=\text{ta}{{\text{n}}^{-1}}\frac{\operatorname{Im}ag\text{(}{{V}_{{{\text{R}}_{2}}}}\text{)}}{\operatorname{Re}al\text{(}{{V}_{{{\text{R}}_{2}}}}\text{)}}-\text{ta}{{\text{n}}^{-1}}\frac{\operatorname{Im}ag\text{(}{{V}_{{{\text{R}}_{1}}}}\text{)}}{\operatorname{Re}al\text{(}{{V}_{{{\text{R}}_{1}}}}\text{)}}$ （1）

图1 PeriScope仪器基础线圈设计

图1b引入倾斜线圈, 线圈法向方向与仪器轴夹角为45° 。利用钻铤旋转, 测量不同方位角情况下接收线圈处的电势, 根据（2）式将信息转化为地下实际地质条件下的相位差和幅度比。该相位差和幅度比包含地层电阻率及界面位置等信息, 因此通常被称为地质信息。

$At{t}'=20\lg \frac{\sqrt{{{[Real({{V}_{{{\beta }_{1}}}})]}^{2}}+{{[Imag({{V}_{{{\beta }_{1}}}})]}^{2}}}}{\sqrt{{{\left[ Real({{V}_{{{\beta }_{2}}}}) \right]}^{2}}+{{[Imag({{V}_{{{\beta }_{2}}}})]}^{2}}}}$

$P{S}'={{\tan }^{-1}}\frac{Imag({{V}_{{{\beta }_{1}}}})}{Real({{V}_{{{\beta }_{1}}}})}-{{\tan }^{-1}}\frac{Imag({{V}_{{{\beta }_{2}}}})}{Real({{V}_{{{\beta }_{2}}}})}$ （2）

1.2 降维处理与反演模型

在电测井仪器探测尺度下, 沿某一方向地层性质可认为保持不变, 这样褶皱等地质构造的电磁波测井资料反演可简化至2D, 地层模型如图2a所示; 即认为地层性质沿x和y方向变化, 沿z方向同一地层性质保持不变。对2D模型进行逐步开窗处理, 将复杂的反演问题转化为一系列的1D地层反演（地层性质仅沿y方向变化, 沿x和z方向同一地层性质保持不变）。经降维处理后, 反演速度的提升主要来自3个方面：①对每个窗口采用1D解析解可以满足实时正演的需要[13], 每个测量点所需计算时间为10-4~10-3 s; ②简化后的1D反演模型待反演参数极少; ③1D反演过程中, 内存等占用量可以忽略不计。对于随钻方位电磁波测井, 1D反演模型包括：多界面反演模型（见图2b, 地层层数多于3层）、双界面反演模型（见图2c, 地层层数为3层）与单界面反演模型（见图2d, 地层层数为2层）。

图2 地层反演模型降维处理（图中不同颜色代表不同地层）

2 梯度反演与随机反演方法

本文首先采用常见的Gauss-Newton算法, 在代价函数中引入自适应乘子正则化项, 以解决传统方法正则化选取困难的问题。同时, 为获取复杂模型的全局最优解, 引入了随机Bayesian算法。

2.1 乘子正则化Gauss-Newton算法

针对随钻方位电磁波测井资料的非线性反演问题, 基于梯度反演方法, 第k步迭代时代价函数${{C}_{k}}\left( m \right)$可表示为[14]：

${{C}_{k}}\left( m \right)=F\left( m \right)\cdot {{R}_{k}}\left( m \right)$ （3）

式中 $F\left( m \right)$为实测数据${{d}^{\text{obs}}}$与正演响应$S\left( m \right)$的L2范数, 即：$F\left( m \right)=\frac{1}{2}{{\left[ S\left( m \right)-{{d}^{\text{obs}}} \right]}^{2}}$。${{R}_{k}}\left( m \right)$表示自适应正则化项：

$\begin{cases} {{R}_{k}}\left( m \right)={{\eta }_{k}}\left( m-{{m}_{k}}^{2}+\delta \right) \\{{\eta }_{k}}={{\left( {{m}_{k}}-{{m}_{p}}^{2}+\delta \right)}^{-1}}\end{cases}$(4)

上式中δ 为一个定值, 其可通过数值实验获取。${{m}_{p}}$为参考模型矢量, 反演过程中可将其设置为上一次迭代时的模型状态, 即：${{m}_{p}}={{m}_{k-1}}$。

采用Gauss-Newton最优化方法对（3）式求解, 令其对m导数为零, 第k次迭代${\partial \text{C}}/{\partial {{m}^{k}}}\; =\text{0}$, 得：

${{m}^{k+1}}={{m}^{k}}-\frac{{{J}^{\text{T}}}\left( {{m}^{k}} \right)\cdot \left[ S\left( {{m}^{k}} \right)-{{\mathbf{d}}^{\text{obs}}} \right]+{{\eta }_{k}}F\left( {{m}^{k}} \right)\cdot \left( {{m}^{k}}-{{m}^{k-1}} \right)}{{{J}^{\text{T}}}\left( {{m}^{k}} \right)\cdot J\left( {{m}^{k}} \right)+{{\eta }_{k}}F\left( {{m}^{k}} \right)+{{\eta }_{k}}\left( {{m}^{k}}-{{m}^{k-1}} \right)\cdot {{J}^{\text{T}}}\left( {{m}^{k}} \right)\cdot \left[ S\left( {{m}^{k}} \right)-{{\mathbf{d}}^{\text{obs}}} \right]}$ （5）

2.2 Bayesian反演算法

根据Bayesian理论, 模型解可以看作一随机变量, 且满足某一概率密度函数分布, 即后验分布。假定先验概率密度$p\left( \mathbf{m} \right)$满足均匀分布, 且似然概率密度函数$p\left( {{\mathbf{d}}^{\text{obs}}}|\mathbf{m} \right)$满足正态分布, 则后验分布$p\left( \mathbf{m}|{{\mathbf{d}}^{\text{obs}}} \right)$可由下式表示[12]：

$p\left( \mathbf{m}|{{\mathbf{d}}^{\text{obs}}} \right)\propto p\left( {{\mathbf{d}}^{\text{obs}}}|\mathbf{m} \right)=$

$\exp \left\{ -\frac{1}{2}{{\left[ {{\mathbf{d}}^{\text{obs}}}-\mathbf{S}\left( \mathbf{m} \right) \right]}^{\text{T}}}\left[ {{\mathbf{d}}^{\text{obs}}}-\mathbf{S}\left( \mathbf{m} \right) \right] \right\}$ （6）

对于Bayesian反演, 问题进一步转化为如何高效采样以获取似然函数分布, 本文采用马科夫链-蒙特卡洛算法（MCMC）实现上述采样。

MCMC的关键在于构建一条马科夫链, 使其满足如下条件：

$p\left( {{\mathbf{m}}_{1}}|{{\mathbf{d}}^{\text{obs}}} \right)q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)=p\left( {{\mathbf{m}}_{2}}|{{\mathbf{d}}^{\text{obs}}} \right)q\left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)$ （7）

为满足上式, 采用“ Metropolis-Hasting” （M-H）算法, 引入接受概率函数$\alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$, 以使：

$p\left( {{\mathbf{m}}_{1}}|{{\mathbf{d}}^{\text{obs}}} \right)q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)\alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)=$

$p\left( {{\mathbf{m}}_{2}}|{{\mathbf{d}}^{\text{obs}}} \right)q\left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)\alpha \left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)$ （8）

式中 $\alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)=\min \left\{ 1, \frac{p\left( {{\mathbf{d}}^{\text{obs}}}|{{\mathbf{m}}_{2}} \right)q\left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)}{p\left( {{\mathbf{d}}^{\text{obs}}}|{{\mathbf{m}}_{1}} \right)q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)} \right\}$

$\alpha \left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)=\min \left\{ 1, \frac{p\left( {{\mathbf{d}}^{\text{obs}}}|{{\mathbf{m}}_{1}} \right)q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)}{p\left( {{\mathbf{d}}^{\text{obs}}}|{{\mathbf{m}}_{2}} \right)q\left( {{\mathbf{m}}_{2}}, {{\mathbf{m}}_{1}} \right)} \right\}$ （9）

MCMC具体实现步骤如下：根据当前迭代模型${{\mathbf{m}}_{1}}$和转移函数$q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$, 生成下一候选模型${{\mathbf{m}}_{2}}$; 然后计算接收概率$\alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$, 并与随机数β 相比较, 若$\beta \ge \alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$则返回模型${{\mathbf{m}}_{2}}$, 反之返回模型${{\mathbf{m}}_{1}}$, 重复上述步骤直至满足终止条件。

3 反演模型与反演算法对比

首先从最常用的双界面模型入手, 探讨其适用范围, 同时探讨Gauss-Newton和Bayesian反演效果的差异; 然后考察单界面、多界面反演的可行性, 并与双界面模型的反演效果进行对比; 最后, 综合反演速度与反演精度, 给出反演模型的选取方法。

3.1 双界面反演

为验证双界面反演模型与反演算法的正确性, 建立如图3a所示的砂泥岩互层模型, 其中砂岩和泥岩电阻率分别为10.0 Ω · m和1.0 Ω · m。模型最上层和最下层层厚为2.5 m, 中间地层层厚均为2 m。仪器自上而下穿过地层时, 视电阻率与幅度比、相位差地质信息的响应见图3b— 3d, 可见在测井仪器穿过地层界面附近时（如图3中虚线所示）, 其响应值均发生了突变。图4a为采用Gauss-Newton算法反演得到的二维窗帘图, 其中每点的颜色表示地层电阻率值, 可见在测井仪器穿过地层界面时（如图中虚线所示）, 其电阻率值发生了明显突变。对比图3a、图4a可知, 图4a中电阻率突变点（如图中虚线所示）所对应的界面位置与图3a地层真实模型基本一致, 从而验证了Gauss-Newton算法的准确性与双界面反演模型的可行性。值得注意的是, 当地层界面距井筒较远时, 难以准确重构远处的地层界面, 但这并不影响地质导向的进行。二维窗帘图不仅给出了仪器附近地层的上/下界面位置, 也直观地初步显示了地层结构分布。图4b为基于Bayesian算法的反演结果, 与图4a中Gauss-Newton反演结果相比, Bayesian算法也可准确确定地层界面位置信息（如图4b中虚线所示）, 但其地层界面处的电阻率有“ 毛边” 现象, 稳定性略差。

图3 层厚为2 m砂泥岩互层模型的PeriScope响应

为研究双界面模型在不同层厚地层的适用性, 建立层厚1.0 m砂泥岩薄互层（最上层和最下层层厚为2.0 m）, 其中砂岩和泥岩电阻率分别为10.0 Ω · m和1.0 Ω · m（见图5a）。分别利用Gauss-Newton和Bayesian算法进行处理, 其反演结果如图5b— 5d, 可见两种算法均可确定地层界面位置, 但在地层界面处电阻率变化的稳定性变差, 尤其后者Bayesian算法反演结果“ 毛边” 现象严重（见图中圆圈处）。对比图4和图5可以得到以下结论：①双界面反演适用于层厚大于1.0 m的地层, 层厚小于1.0 m时, 反演结果的稳定性变差; ②Gauss- Newton反演结果略优于Bayesian结果, 对于Bayesian算法而言, 不断增加采样次数可以提高反演精度。

图4 中间层厚为2.0 m的砂泥岩互层双界面反演模型反演结果

图5 中间层厚为1.0 m的砂泥岩薄互层双界面反演模型反演结果

3.2 单界面反演

为验证单界面反演模型的可行性与适用性, 建立5层砂泥互层模型。基于Gauss-Newton算法与基于单、双界面模型, 反演后的地层界面位置如图6所示, 图中红、蓝色点分别为双、单界面模型反演地层界面位置。对层厚为3 m的储集层（见图6a）, 仪器距邻近地层界面距离小于0.9 m时, 单界面反演可获得准确结果。当仪器距邻近地层界面距离大于0.9 m时, 方位地质信息受多个界面的影响, 导致单界面反演得到的仪器距邻近地层界面距离偏小。随着地层厚度增大, 单界面反演结果准确性和稳定性不断提高, 地层厚度大于4.0 m时（见图6b— 6d）, 基于单界面和双界面模型的反演结果基本一致（反演得到的地层界面位置基本一致）。一般而言, 当储集层厚度大于4.0 m时, 基于单界面模型能够准确反演距离仪器1.8 m范围内地层上、下界面, 可以满足实时地质导向的需要。同样也可基于Bayesian算法和单界面反演模型获取与图6一致的反演效果, 此处不再赘述。

图6 单界面模型与双界面模型反演得到的地层界面位置

3.3 四界面反演

在层厚小于1.0 m时, 仪器探测范围内可能存在多个界面, 过度简化的单界面和双界面模型不再适用。此时, 需采用更为复杂的多界面模型方可精细刻画仪器附近的地层界面分布。而模型层数的增加一方面导致反演参数变多, 另一方面导致代价函数局部极小值数量迅速增加, 梯度类算法不再适用, 故在本节中笔者仅考虑Bayesian反演方法。为说明多界面反演模型的优越性, 建立图7a所示模型, 其中中间层为0.6 m的高电阻率薄层（在测井轨迹上对应的横向距离分别约为135 m和165 m）。分别采用单界面、双界面和四界面模型, 反演得到的二维窗帘图见图7b— 7d。由图7b可以看出, 单界面反演仅能在上、下无限厚地层处获取邻近地层界面位置, 难以确定中间薄层的存在（在横向距离135 m和165 m处电阻率几乎无突变）; 使用双界面模型时（见图7c）, 仪器由上而下接近或远离中间薄层时, 可得到靠近仪器的一侧界面位置, 但仪器在中间薄层时, 无法准确确定薄层上下界面位置（在横向距离135 m、165 m处电阻率值突变明显, 但135~165 m处电阻率值不稳定）; 相比之下, 采用四界面反演模型（见图7d）, 不仅可以给出仪器所在层的上下界面, 还能得到邻层界面（在横向距离135 m、165 m处电阻率值突变明显, 135~165 m处电阻率值稳定）。以横向距离135 m处为例, 经四界面反演后可见在垂直深度1.0、2.3和2.9 m处存在3个电阻率突变界面（见图7d横向虚线）, 即存在3个地层界面, 与图7a中地层界面位置一致。四界面模型精确重构了薄层处地层结构, 其可用于储集层结构的精细评价。

图7 四界面模型及Bayesian反演后的二维窗帘图

3.4 反演模型与反演算法优选

反演速度是决定反演模型选取的另一重要因素。Gauss-Newton反演速度主要取决于初始值的数量及迭代次数, 而Bayesian反演则是取决于MCMC采样量。对单界面和双界面反演模型, Gauss-Newton算法每个测量点用时分别为0.3 s和3.0 s, 而Bayesian需要采样点数分别为5 000和20 000次, 每个测量点耗时15 s和60 s。对四界面反演模型, Bayesian需要进行100 000次采样耗时500 s。综合考虑反演速度、反演精度与反演模型最简化的原则, 反演模型与反演算法选取可参见表1。

表1 反演模型与反演算法选取表
4 实际资料处理

本节通过图8所示的某水平井实例说明随钻方位电磁波测井仪器在地质导向中的应用。该水平井可分为3部分：A段, 水平井着陆; B段, 水平井钻穿靶层下界面后, 调整钻进方向重新钻入砂岩储集层; C段, 控制井眼轨迹, 使其一直处于目的层内。

图8 基于随钻方位电磁波测井资料的水平井实时解释结果（图8e中不同颜色代表不同地层）

若无相位差和幅度比地质信息曲线时, 一般利用实测电阻率、岩性曲线与邻近导眼井提取的地层序列进行相关对比, 建立地层模型。在A段和B段, 基于实时对比结果可以较好地确定地层界面与井眼轨迹的相对位置关系。而在C段横向距离2 455~2 535 m处, 随钻电磁波测井曲线出现异常下降（见图8d中两条蓝色虚线间电阻率）, 仅通过相关对比难以判断低电阻率围岩上下位置, 测井解释存在不确定性且钻井时很容易穿出目的层。相比而言, 利用随钻方位电磁波测井（见图8a— 8c蓝色虚线间测井曲线）可提前预知地层上界面位置（随钻方位电磁波测井值稍变大, 远小于钻穿目的层时异常值, 说明接近目的层上界面）, 从而调整井眼轨迹以避免钻出目的层, 同时也给出了明确的地层结构分布。

5 结论

通过降维可以把随钻方位电磁波测井资料3D反演问题转化为更为可行的1D反演处理, 进而实现对地层界面位置信息的实时确定。当靶层厚度大于4.0 m时, 采用单界面反演模型可获取距离仪器最近的地层界面; 靶层厚度为1.0~4.0 m时, 采用双界面反演模型可获取准确的地层上、下界面位置与储集层电阻率信息; 当靶层厚度小于1.0 m时, 单/双界面反演模型不再适用, 此时需采用多界面反演模型。对单界面和双界面反演模型, 采用正则化Gauss-Newton方法即可实现对随钻方位电磁波测井资料的地层界面实时提取, 每个采样点反演耗时分别为0.3 s和3.0 s; 对多界面反演模型, 基于MCMC采样的Bayesian最优化算法可实现地层界面的精细解释, 每个采样点反演耗时为500 s。实际资料处理结果表明, 基于优选后的反演模型和反演算法可以实现邻近地层界面的精确提取; 随钻方位电磁波测井降低了传统电磁波测井水平井解释结果的不确定性, 大大提高了随钻地质导向能力。

致谢：感谢美国德克萨斯大学奥斯汀分校Carlos Torres-Verdin教授在反演算法上的指导, 感谢Chevron公司汪涵明博士在方位电磁波测井数据处理过程中提供的帮助。

符号注释：

Att— — 幅度比, dB; Att° — — 实际条件下幅度比, dB; Att° SAD1, Att° SAD4— — 频率100 kHz和400 kHz、源距243.84 cm（96 in）的幅度比, dB; Att° SAS4— — 频率400 kHz、源距86.36 cm（34 in）的幅度比, dB; ${{C}_{k}}\left( m \right), {\partial \text{C}}/{\partial m}\; $— — 代价函数及其导数; ${{d}^{\text{obs}}}$— — 实测数据; $F\left( m \right)$— — 实测数据与正演响应的L2范数; H— — 靶层厚度, m; Hup, Hdown— — 仪器记录点与地层上、下界面的距离, m; Imag— — 取虚部函数; J— — 雅克比矩阵; k— — 迭代次数; m— — 待反演的参数矢量; m1、m2— — 任意模型; N— — 地层界面数; $p\left( \mathbf{m} \right)$— — 先验概率密度函数; $p\left( {{\mathbf{d}}^{\text{obs}}}|\mathbf{m} \right)$— — 似然概率密度函数; $p\left( \mathbf{m}|{{\mathbf{d}}^{\text{obs}}} \right)$— — 后验概率密度函数; PS— — 相位差, （° ）; PS° — — 实际条件下相位差, （° ）; PS° SPD1, PS° SPD4— — 频率100 kHz和400 kHz、源距243.84 cm（96 in）的相位差, （° ）; PS° SPS4— — 频率400 kHz、源距86.36 cm（34 in）的相位差, （° ）; $q\left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$— — 转移核函数; Real— — 取实部函数; RA28H, RP28H— — 频率2 MHz、源距71.12 cm（28 in）的幅度比和相位差电阻率, Ω · m; RA40H, RP40H— — 频率2 MHz、源距101.60 cm（40 in）的幅度比和相位差电阻率, Ω · m; ${{R}_{k}}\left( m \right)$— — 自适应正则化项; Rs1, Rs2— — 单界面地层模型中围岩电阻率, Ω · m; Rt— — 仪器所在层电阻率, Ω · m; Rup, Rdown— — 双界面地层模型上、下围岩等效电阻率, Ω · m; $S\left( m \right)$— — 正演响应; T— — 矩阵转置; ${{V}_{{{\text{R}}_{1}}}}$, ${{V}_{{{\text{R}}_{2}}}}$— — 同轴发射线圈在两个同轴接收线圈处产生的电动势, V; ${{V}_{{{\beta }_{1}}}}$, ${{V}_{{{\beta }_{2}}}}$— — 方位角为β 1和β 2时倾斜接收线圈处的电动势, V; x, y, z— — 直角坐标系; $\alpha \left( {{\mathbf{m}}_{1}}, {{\mathbf{m}}_{2}} \right)$— — 接受概率函数; β — — 随机数; β 1, β 2— — 钻铤的旋转方位角, （° ）; θ — — 仪器与地层夹角, （° ）; δ — — 实验获取初始值。

The authors have declared that no competing interests exist.

参考文献
[1] BITTAR M, AKI A. Advancement and economic benefit of geosteering and well-placement technology[J]. The Leading Edge, 2015, 34(5): 524-528. [本文引用:1]
[2] LI Q, OMERAGIC D, CHOU L, et al. New directional electromagnetic tool for proactive geosteering and accurate formation evaluation while drilling[R]. New Orleans, USA: SPWLA 46th Annual Logging Symposium, 2005. [本文引用:1]
[3] LI H, WANG H. Investigation of eccentricity effects and depth of investigation of azimuthal resistivity LWD tools using 3D finite difference method[J]. Journal of Petroleum Science & Engineering, 2016, 143: 211-225. [本文引用:1]
[4] 范宜仁, 吴俊晨, 吴飞, 等. 地层模块尺度下的新型钻井液侵入物理模拟系统[J]. 石油勘探与开发, 2017, 44(1): 125-129.
FAN Yiren, WU Junchen, WU Fei, et al. A new physical simulation system of drilling mud invasion in formation module[J]. Petroleum Exploration and Development, 2017, 44(1): 125-129. [本文引用:1]
[5] 王磊, 范宜仁, 黄瑞, 等. 各向异性介质多分量感应测井三维Born几何因子理论研究[J]. 物理学报, 2015, 64(23): 430-440.
WANG Lei, FAN Yiren, HUANG Rui, et al. Three dimensional Born geometrical factor of multi-component induction logging in anisotropic media[J]. Acta Physica Sinica, 2015, 64(23): 430-440. [本文引用:1]
[6] 谭茂金, 高杰, 邹友龙, 等. 盐水泥浆条件下定向井双侧向测井环境校正方法研究[J]. 地球物理学报, 2012, 55(4): 1422-1432.
TAN Maojin, GAO Jie, ZOU Youlong, et al. Environment correction method of dual laterolog in directional well[J]. Chinese Journal of Geophysics, 2012, 55(4): 1422-1432. [本文引用:1]
[7] ABUBAKAR A, HABASHY T M, DRUSKIN V L, et al. 2. 5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements[J]. Geophysics, 2008, 73(4): 165-177. [本文引用:1]
[8] THIEL M, OMERAGIC D. High fidelity real-time imaging with electromagnetic logging-while-drilling measurements[J]. IEEE Transactions on Computational Imaging, 2017, 3(2): 369-378. [本文引用:1]
[9] YANG J, OMERAGIC D, LIU C B, et al. Bed-boundary effect removal to aid formation resistivity interpretation from LWD Propagation measurements at all relative dip angles[R]. New Orleans, USA: SPWLA 46th Annual Logging Symposium, 2005. [本文引用:1]
[10] LI H, ZHOU J. Distance of detection for LWD deep and ultra-deep azimuthal resistivity tools[R]. Oklahoma, USA: SPWLA 58th Annual Logging Symposium, 2017. [本文引用:1]
[11] OMERAGIC D, HABASHY T, ESMERSOY C, et al. Real-time interpretation of formation structure from azimuthal EM measurements[R]. Veracruz, Mexico: SPWLA 47th Annual Logging Symposium, 2006. [本文引用:1]
[12] YANG Q, TORRES-VERDÍN C. Efficient 2D Bayesian inversion of borehole resistivity measurements[J]. SEG Technical Program Expand ed Abstracts, 2011, 30(1): 4424. [本文引用:2]
[13] HONG D, XIAO J, ZHANG G, et al. Characteristics of the sum of cross-components of triaxial induction logging tool in layered anisotropic formation[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(6): 3107-3115. [本文引用:1]
[14] HABASHY T M, ABUBAKAR A. A general framework for constraint minimization for the inversion of electromagnetic measurements[J]. Progress in Electromagnetics Research, 2004, 46: 265-312. [本文引用:1]