Recursive Algorithm for Electromagnetic Fields from Magnetic Dipole in Layered Triaxial Anisotropic Medium and Its Application
-
摘要:
为快速、精确模拟复杂三轴各向异性介质中的电磁波测井响应规律,研究了一种适用于水平层状地层的磁偶极子源电磁场伪解析递推算法。该算法通过双重傅里叶变换,将空间域的三维电磁场正演转换为一系列对一维谱域场的求解,将地层上、下界面位置引入谱域电磁场通解公式,克服了传统方法存在的数值溢出问题;进一步利用传播矩阵方法,递推获得各个地层界面处的幅度系数;针对积分核函数存在的强烈振荡问题,提出了一种双重正余弦数值滤波积分方法,实现了谱域到空间域电磁场的准确、快速转换。模拟结果表明:三轴各向异性导致电测井响应更为复杂,传统各向异性电测井解释模型不再适用。新递推算法是复杂各向异性储层测井响应分析的基础,也为电性参数准确提取和精准地质导向提供了正演手段。
Abstract:In order to quickly and accurately simulate the response of electromagnetic (EM) wave logging in complex triaxial anisotropic media, a pseudo-analytical recursive algorithm for the EM fields generated by magnetic dipole was developed, which was suitable for horizontal layered strata. By using the dual Fourier transform technique, the algorithm transformed the forward modeling of the three-dimensional EM field into a series of solutions of a one-dimensional spectral field. The upper and lower boundary positions were introduced into the general expressions of the spectral field, which avoided the numerical overflow problem in traditional algorithms. The propagation matrix method was then adopted to recursively derive the amplitude coefficients at each boundary. To handle the strong oscillation of the integral kernel function, a dual sine/cosine numerical filtering integral technique was developed to achieve the accurate and fast conversion between EMs from the spectral domain to the spatial domain. The simulation results show that the existence of triaxial anisotropy makes the electrical logging response more complex, and the traditional anisotropic electrical logging interpretation model is no longer applicable. The new recursive algorithm is the basis of complex anisotropic reservoir logging response analysis and provides a forward modeling method for accurate extraction of electrical parameters and accurate geosteering.
-
致密砂岩、页岩等地下岩石的发育具有广泛的各向异性特征[1–4]。沉积过程中,地层电性通常沿水平方向保持不变,而在纵向上存在非均质性,该类地层称为横向各向同性介质(VTI)[5–10]。若各向同性地层中存在相互平行的垂直裂缝,则会引起横向电导率的各向异性,该类地层常等效为水平向横向各向同性介质(HTI)[11]。页岩及致密岩性复杂地层既存在由于沉积或成岩作用引起的层理各向异性,又存在裂缝发育引起的诱导各向异性[12],不同方向上的电导率均不相同,表现为更复杂的三轴各向异性[13]。电测井是评价地层各向异性的重要手段,其本质是测量偶极子源激发的磁场[14–15]。传统感应测井仅能探测地层水平电阻率;多分量感应测井通过测量磁场张量,可在任意井斜获取地层的各向异性。随钻方位电磁波测井则借助于磁场交叉分量的方位敏感性,实现了对未钻遇地层实时“环视”和“前视”[16–17]。然而,目前电测井仪研究主要针对VTI介质,对HTI及三轴各向异性介质响应规律的研究极少。
准确求解层状介质中偶极子源电磁场,不仅能够指导仪器设计、储层定性解释,也是反演的基础。早期,研究人员主要针对VTI介质提出并发展了多种求解方法[18–20]。相比于VTI介质,三轴各向异性介质电磁场的横电波与横磁波相互耦合,电磁波表达形式更复杂,需求解双重无穷域振荡积分,积分难度大。L. Løseth等人[21] 给出了三轴各向异性地层广义反射矩阵递推公式。Hu Yunyun等人[22]提出了基于高斯离散点和一次场扣除技术的双重积分方法。总体而言,现有的三轴各向异性地层磁偶极子源电磁场求解方法较为单一,双重积分方法所需节点数仍较多,积分效率仍有待进一步提升。
笔者首先推导了磁偶极子源在一维地层的谱域电磁场的通解形式,通过将地层上/下界面引入至下/上行波的指数项,避免了数值溢出问题;然后将系数传播矩阵推广至三轴各向异性地层,获得了各层内的幅度系数递推表达式,并提出了双重无限域振荡核函数的快速积分方法,实现了任意位置空间域电磁场的准确计算;最后,将新的伪解析算法应用于电磁波测井,分析了不同各向异性地层模型对感应和随钻方位电磁波测井响应的影响。
1. 层状各向异性介质电磁场伪解析解
1.1 控制偏微分方程
假定空间中仅存在低频磁偶极子,且其按
e−iωt 形式随时间变化,相应的麦克斯韦方程为:\left\{\begin{array}{l} \nabla \times {\boldsymbol{E}} = {\rm{i}}\omega \mu \left( {{\boldsymbol{M}} + {\boldsymbol{H}}} \right)\\ \nabla \times {\boldsymbol{H}} = {\boldsymbol{\sigma E}} \end{array}\right. (1) 式中:H为磁场强度,A/m;E为电场强度,V/m;μ为介质的磁导率,H/m;ω为介质的角频率,rad/s;M为磁矩,A·m2;σ为三轴各向异性介质的电导率,S/m。
水平层状介质仅在z方向存在电性非均质性。对式(1)作关于x和y的傅里叶变换,并消除谱域场纵向分量,分别得到满足谱域电场
{\tilde E_x} 和谱域磁场{\tilde H_x} 的微分方程:\left\{ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^4}{{\tilde E}_x}}}{{\partial {z^4}}} - {a_{\rm{a}}}\dfrac{{{\partial ^2}{{\tilde E}_x}}}{{\partial {z^2}}} + {a_b}{{\tilde E}_x} = \dfrac{{{\partial ^2}{e_c}}}{{\partial {z^2}}} + {e_b}\dfrac{{\partial {b_c}}}{{\partial z}} - {e_c}{b_a}} \\ {\dfrac{{{\partial ^4}{{\tilde H}_x}}}{{\partial {z^4}}} - {a_a}\dfrac{{{\partial ^2}{{\tilde H}_x}}}{{\partial {z^2}}} + {a_b}{{\tilde H}_x} = \dfrac{{{\partial ^2}{b_c}}}{{\partial {z^2}}} + {b_b}\dfrac{{\partial {e_c}}}{{\partial z}} - {e_a}{b_c}} \end{array}} \right. (2) \, 其中\;\;{a_{\rm{a}}} = - \frac{{k_z^2\left( {k_x^2 + k_y^2} \right) - {\xi ^2}\left( {k_x^2 + k_z^2} \right) - {\eta ^2}\left( {k_y^2 + k_z^2} \right)}}{{k_z^2}} (3) {a_{\rm{b}}} = \dfrac{{\left( {{\xi ^2} + {\eta ^2} - k_z^2} \right)\left( {{\xi ^2}k_x^2 + {\eta ^2}k_y^2 - k_x^2k_y^2} \right)}}{{k_z^2}} (4) {e_{\rm{a}}}{\text{ = }}\dfrac{{\left( {{\xi ^2}k_x^2 + {\eta ^2}k_y^2 - k_x^2k_y^2} \right)\left( {k_z^2 - {\xi ^2}} \right)}}{{k_z^2\left( {k_y^2 - {\xi ^2}} \right)}} (5) {e_{\rm{b}}}{\text{ = }}\dfrac{{{\rm{i}}\omega \mu \xi \eta \left( {k_y^2 - k_z^2} \right)}}{{k_z^2\left( {k_y^2 - {\xi ^2}} \right)}} (6) {{e_{\rm{c}}} = {\rm{i}}\omega \mu {M_y}\delta \left( z \right) + \dfrac{{\omega \mu \eta k_y^2\left( {k_z^2 - {\xi ^2}} \right)}}{{k_z^2\left( {k_y^2 - {\xi ^2}} \right)}}{M_z}\delta \left( z \right)} (7) {b_{\rm{a}}}{\text{ = }}\dfrac{{\left( {k_y^2 - {\xi ^2}} \right)\left( {{\xi ^2} + {\eta ^2} - k_z^2} \right)}}{{k_z^2 - {\xi ^2}}} (8) {b_{\rm{b}}}{\text{ = }}\dfrac{{\xi \eta \left( {k_z^2 - k_y^2} \right)}}{{{\rm{i}}\omega \mu \left( {k_z^2 - {\xi ^2}} \right)}} (9) {{b_{\rm{c}}} = \left( {{\xi ^2} - k_y^2} \right){M_x}\delta \left( z \right) + \dfrac{{\xi \eta \left( {k_y^2 - {\xi ^2}} \right)}}{{k_z^2 - {\xi ^2}}}{M_y}\delta \left( z \right) - {\rm{i}}\xi {M_z}\delta \left( z \right)} (10) k_j^2 = {\rm{i}}\omega \mu {\sigma _j} (11) 式中:ξ和η为x和y方向的波数,cm−1;δ(z)为狄拉克函数。
1.2 谱域电磁场的通解形式
式(2)的解由非齐次方程的特解和齐次方程的通解叠加而成,前者是源在背景介质中激发的一次电磁场,后者则为地层非均质引起的散射场。为获取谱域一次场,对式(2)沿z向进行傅里叶变换和逆变换。首先,采用围线积分方法,获得谱域电磁场的表达式。以z>0为例,
\tilde E_x^{\rm{b}}\left( {\boldsymbol{k}} \right) 和\tilde H_x^{\rm{b}}\left( {\boldsymbol{k}} \right) 写为:\begin{split} \tilde E_x^{\rm{b}} =& \frac{{\text{π}} }{{k_z^2}}\left[ {{{\left. {\frac{{{{\rm{e}}^{ - \alpha z}}\left( {{a_1} + {a_2} + {a_3}} \right)}}{{\alpha \left( {{\beta ^2} - {\alpha ^2}} \right)}}} \right|}_{\varsigma = {\rm{i}}\alpha }} + {{\left. {\frac{{{{\rm{e}}^{ - \beta z}}\left( {{a_1} + {a_2} + {a_3}} \right)}}{{\beta \left( {{\alpha ^2} - {\beta ^2}} \right)}}} \right|}_{\varsigma = {\rm{i}}\beta }}} \right], \\ &\qquad \qquad\qquad\qquad\quad z > 0 \end{split} (12) \begin{split} \tilde H_x^{\rm{b}} =& \frac{{\text{π}} }{{k_z^2}}\left( {{{\left. {\frac{{{{\rm{e}}^{ - \alpha z}}\left( {{b_1} + {b_2} + {b_3}} \right)}}{{\alpha \left( {{\beta ^2} - {\alpha ^2}} \right)}}} \right|}_{\varsigma = {\rm{i}}\alpha }} + {{\left. {\frac{{{{\rm{e}}^{ - \beta z}}\left( {{b_1} + {b_2} + {b_3}} \right)}}{{\beta \left( {{\alpha ^2} - {\beta ^2}} \right)}}} \right|}_{\varsigma = {\rm{i}}\beta }}} \right),\\ &\qquad\qquad\qquad\qquad\quad z > 0 \end{split} (13) 式中:α和β分别为I型波和II型波的波数,cm−1。
\alpha ,\beta = \sqrt {{{\left( {{a_{{a}}} \mp \sqrt {a_{{a}}^2 - 4{a_{{b}}}} } \right)} \mathord{\left/ {\vphantom {{\left( {{a_{{a}}} \mp \sqrt {a_{{a}}^2 - 4{a_{{b}}}} } \right)} 2}} \right. } 2}} (14) 式(12)和式(13)中的a1,a2,a3,b1,b2和b3为中间变量,表达式如下:
\left( {\begin{array}{*{20}{c}} {{a_1}} \\ {{a_2}} \\ {{a_3}} \end{array}} \right) = \left( {\begin{array}{*{20}{l}} {\omega \mu \xi \eta \varsigma \left( {k_y^2 - k_z^2} \right){M_x}} \\ {\omega \mu \varsigma \left( {k_z^2\left( {{\xi ^2} + {\varsigma ^2} - k_y^2} \right) + k_y^2{\eta ^2}} \right){M_y}} \\ {\omega \mu \eta \left( {k_y^2k_z^2 - {\varsigma ^2}k_z^2 - k_y^2\left( {{\xi ^2} + {\eta ^2}} \right)} \right){M_z}} \end{array}} \right) (15) \left( {\begin{array}{*{20}{c}} {{b_1}} \\ {{b_2}} \\ {{b_3}} \end{array}} \right) = \left( {\begin{array}{*{20}{l}} {\left( {{\varsigma ^2}k_z^2\left( {k_y^2 - {\xi ^2}} \right) + \left( {{\xi ^2}k_x^2 + {\eta ^2}k_y^2 - k_x^2k_y^2} \right)\left( {k_z^2 - {\xi ^2}} \right)} \right){M_x}} \\ {\xi \eta \left( { - {\varsigma ^2}k_z^2 - \left( {{\xi ^2}k_x^2 + {\eta ^2}k_y^2 - k_x^2k_y^2} \right)} \right){M_y}} \\ {\xi \varsigma \left( { - {\varsigma ^2}k_z^2 + k_x^2k_z^2 - {\eta ^2}k_y^2 - {\xi ^2}k_x^2} \right){M_z}} \end{array}} \right) (16) 然后,求取谱域散射场的一般表达式,其通解满足指数形式
{{\rm{e}}^{\lambda z}} ,则式(12)变为:{\lambda ^4} - {a_{\rm{a}}}{\lambda ^2} + {a_{\rm{b}}} = 0 (17) 式(17)存在4个互不相同的解α,–α,β和–β。将地层界面位置引入指数项,第j层散射场通解的表达式为:
\begin{split} \left( {\begin{array}{*{20}{c}} {\tilde E_{xj}^{\rm{s}}} \\ {\tilde H_{xj}^{\rm{s}}} \end{array}} \right) = &\left( {\begin{array}{*{20}{c}} {{C_j}} \\ {{U_j}} \end{array}} \right){{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + \left( {\begin{array}{*{20}{c}} {{D_j}} \\ {{V_j}} \end{array}} \right){{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} +\\ &\left( {\begin{array}{*{20}{c}} {{G_j}} \\ {{W_j}} \end{array}} \right){{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} + \left( {\begin{array}{*{20}{c}} {{H_j}} \\ {{X_j}} \end{array}} \right){{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}}\\[-14pt] \end{split} (18) 式中:zj–1和zj分别为第j层地层下界面和上界面的位置,m;Cj,Dj,Gj和Hj为界面处的电场幅度,V/m;Uj,Vj,Wj和Xj为界面处的磁场幅度,A/m。
式(3)—式(10)中的指数项能保证波总是在衰减的,避免了数值溢出现象。
将一次场和二次场合并,得到磁偶极子源在三轴各向异性介质的通解:
\left\{ {\begin{array}{*{20}{l}} \begin{split} {{\tilde E}_{xj}} =& {C_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + {D_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} + \\ &{G_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} +{H_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + \gamma \tilde E_x^{\rm{b}} \end{split} \\ \begin{split} {{\tilde H}_{xj}} =& {U_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + {V_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} + \\ &{W_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} + {X_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + \gamma \tilde H_x^{\rm{b}} \end{split} \end{array}} \right. (19) 其中,
\gamma =\left\{\begin{array}{c}0\;\; 源不在第j层\\ 1\;\; 源在第j层时\end{array} \right. 。进一步,可得y方向的谱域分量:\begin{split} {\tilde E_y} = &a_{1_j}\left( {C_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + {D_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} + {G_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} +\right.\\ &\left.{H_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + \gamma \tilde E_x^{\rm{b}} \right) + \left( - a_{2_j}{U_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} - \right.\\ &\left. a_{3_j}{V_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} +a_{2_j}{W_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}}+\right. \\ & \left. a_{3_j}{X_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + \gamma a_{4_j} \right)\\[-12pt] \end{split} (20) \begin{split} {\tilde H_y} =& \left( b_{2_j}{C_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + b_{3_j}{D_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} -\right. \\ & \left.b_{2_j}{G_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} - b_{3_j}{H_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + b_{4_j}\gamma \right) +\\ &b_{1_j}\left( {U_j}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}} + {V_j}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}} +\right. \\ & \left.{W_j}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}} + {X_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}} + \gamma \tilde H_x^{\rm{b}} \right) \end{split} (21) \,其中\qquad\qquad\qquad\quad a_{1_j} = - \dfrac{{\xi \eta }}{{k_y^2 - {\xi ^2}}} \quad (22) a_{2_j} = \dfrac{{{\rm{i}}\omega \mu {\alpha _j}}}{{k_y^2 - {\xi ^2}}} (23) a_{3_j} = \dfrac{{{\rm{i}}\omega \mu {\beta _j}}}{{k_y^2 - {\xi ^2}}} (24) a_{4_j} = \dfrac{{{\rm{i}}\omega \mu }}{{k_y^2 - {\xi ^2}}} (25) b_{1_j} = - \dfrac{{\xi \eta }}{{k_z^2 - {\xi ^2}}} (26) b_{2_j} = - \dfrac{{k_z^2{\alpha _j}}}{{{\rm{i}}\omega \mu \left( {k_z^2 - {\xi ^2}} \right)}} (27) b_{3_j} = - \dfrac{{k_z^2{\beta _j}}}{{{\rm{i}}\omega \mu \left( {k_z^2 - {\xi ^2}} \right)}} (28) b_{4_j} = \dfrac{{k_z^2}}{{\rm{i}}\omega \mu \left( {k_z^2 - {\xi ^2}}\right)} (29) 式 (19)—式(21)中存在8个未知幅度,其中仅有4个系数独立,具体推导公式见附录A。第j层中,谱域水平电磁场分量可表达为:
{\left[ {{{\tilde E}_{xj}},{{\tilde H}_{xj}},{{\tilde E}_{yj}},{{\tilde H}_{yj}}} \right]^{\rm{T}}} = {{\boldsymbol{P}}_j}\left( z \right){{\boldsymbol{\varLambda }}_j} + \gamma {\boldsymbol{S}}\left( z \right) (30) \, 其中\qquad {\boldsymbol{S}}\left( z \right) = {\left[ {{k_s}\left( z \right),{l_{\rm{s}}}\left( z \right),{m_{\rm{s}}}\left( z \right),{n_{\rm{s}}}\left( z \right)} \right]^{\rm{T}}} (31) {{\boldsymbol{\varLambda }}_j} = {\left[ {{C_j},{G_j},{V_j},{X_j}} \right]^{\rm{T}}} (32) 式中:
{{\boldsymbol{\varLambda }}_j} 为每层待求幅度系数构成的矢量。第j层的幅度传播矩阵
{{\boldsymbol{P}}_j}\left( z \right) 为:{{\boldsymbol{P}}_j}\left( z \right) = \left( {\begin{array}{*{20}{c}} {{k_{{\rm{a}},j}}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}}}&{{k_{{\rm{b}},j}}{G_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}}}&{{k_{{\rm{c}},j}}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}}}&{{k_{{\rm{d}},j}}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}}} \\ {{l_{{\rm{a}},j}}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}}}&{{l_{{\rm{b}},j}}{G_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}}}&{{l_{{\rm{c}},j}}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}}}&{{l_{{\rm{d}},j}}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}}} \\ {{m_{{\rm{a}},j}}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}}}&{{m_{{\rm{b}},j}}{G_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}}}&{{m_{{\rm{c}},j}}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}}}&{{m_{{\rm{d}},j}}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}}} \\ {{n_{{\rm{a}},j}}{{\rm{e}}^{ - {\beta _j}\left( {z - {z_{j - 1}}} \right)}}}&{{n_{{\rm{b}},j}}{G_j}{{\rm{e}}^{{\beta _j}\left( {z - {z_j}} \right)}}}&{{n_{{\rm{c}},j}}{{\rm{e}}^{ - {\alpha _j}\left( {z - {z_{j - 1}}} \right)}}}&{{n_{{\rm{d}},j}}{{\rm{e}}^{{\alpha _j}\left( {z - {z_j}} \right)}}} \end{array}} \right) (33) 式中:k,l,m和n取决于地层的电性参数和角频率,其具体表达见式(A-29)—式(A-48)。
1.3 多层介质中电磁波的传播与幅度系数递推方法
假定三轴各向异性地层有N层,源位于第m层,如图1所示。根据边界条件,界面两侧电磁场的切向分量连续,z=zj处有:
{\left[ {{{\tilde E}_{xj}},{{\tilde H}_{xj}},{{\tilde E}_{yj}},{{\tilde H}_{yj}}} \right]^{\rm{T}}} = {\left[ {{{\tilde E}_{x\left( {j + 1} \right)}},{{\tilde H}_{x\left( {j + 1} \right)}},{{\tilde E}_{y\left( {j + 1} \right)}},{{\tilde H}_{y\left( {j + 1} \right)}}} \right]^{\rm{T}}} (34) 若第j和第j+1层均为无源层,2层之间的幅度系数递推关系如下:
{{\boldsymbol{P}}_j}\left( {{z_j}} \right){{\boldsymbol{\varLambda }}_j} = {{\boldsymbol{P}}_{j + 1}}\left( {{z_j}} \right){{\boldsymbol{\varLambda }}_{j + 1}} (35) 式中:
{{\boldsymbol{P}}_j}\left( {{z_j}} \right) 及{{\boldsymbol{P}}_{j + 1}}\left( {{z_j}} \right) 为j层上、下界面处的幅度传播矩阵。同理,在第m层上界面和下界面处,式(33)变为:
{{\boldsymbol{P}}_{p - 1}}\left( {{z_{p - 1}}} \right){{\boldsymbol{\varLambda }}_{p - 1}} = {{\boldsymbol{P}}_p}\left( {{z_{p - 1}}} \right){{\boldsymbol{\varLambda }}_p} + {\boldsymbol{S}}\left( {{z_{p - 1}}} \right) (36) {{\boldsymbol{P}}_p}\left( {{z_p}} \right){{\boldsymbol{\varLambda }}_p} + {\boldsymbol{S}}\left( {{z_p}} \right) = {{\boldsymbol{P}}_{p + 1}}\left( {{z_p}} \right){{\boldsymbol{\varLambda }}_{p + 1}} (37) 利用有源区和无源区的边界条件,经逐层递推可以得到第一层和最后一层幅度系数的矢量关系式:
{{\boldsymbol{P}}_1}\left( {{z_1}} \right){{\boldsymbol{\varLambda }}_1} = {\boldsymbol{A}}{{\boldsymbol{P}}_N}\left( {{z_{N - 1}}} \right){{\boldsymbol{\varLambda }}_N} - {\boldsymbol{BS}}\left( {{z_p}} \right) + {\boldsymbol{CS}}\left( {{z_{p + 1}}} \right) (38) 式中:A,B和C为4阶方阵。
{\boldsymbol{A}} = {{\boldsymbol{P}}_2}\left( {{z_1}} \right){\left( {{{\boldsymbol{P}}_2}\left( {{z_2}} \right)} \right)^{ - 1}} \cdots {{\boldsymbol{P}}_{N - 1}}\left( {{z_{N - 2}}} \right){\left( {{{\boldsymbol{P}}_{N - 1}}\left( {{z_{N - 1}}} \right)} \right)^{ - 1}} (39) {\boldsymbol{B}} = {{\boldsymbol{P}}_2}\left( {{z_1}} \right){\left( {{{\boldsymbol{P}}_2}\left( {{z_2}} \right)} \right)^{ - 1}} \cdots {{\boldsymbol{P}}_p}\left( {{z_{p - 1}}} \right){\left( {{{\boldsymbol{P}}_p}\left( {{z_p}} \right)} \right)^{ - 1}} (40) {\boldsymbol{C}} = {{\boldsymbol{P}}_2}\left( {{z_1}} \right){\left( {{{\boldsymbol{P}}_2}\left( {{z_2}} \right)} \right)^{ - 1}} \cdots {{\boldsymbol{P}}_{p - 1}}\left( {{z_{p - 2}}} \right){\left( {{{\boldsymbol{P}}_{p - 1}}\left( {{z_{p - 1}}} \right)} \right)^{ - 1}} (41) 由于最外层是半无限厚介质,
{{\boldsymbol{\varLambda }}_1} 和{{\boldsymbol{\varLambda }}_N} 中各有2个值非零。求解式(37),得到最外层的幅度系数,结合式(35)—式(37),即可逐层递推各层的幅度系数谱域电磁场的切向和垂向分量。最后,沿ξ和η方向进行双重傅里叶逆变换,得到空间域电磁场响应:\begin{split} &\qquad\quad \left[ {{\boldsymbol{E}}\left( {x,y,z} \right),{\boldsymbol{H}}\left( {x,y,z} \right)} \right] = \frac{1}{{{{\left( {2{\text{π}} } \right)}^2}}}\cdot \\ & \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\left[ {{{\tilde {\boldsymbol{E}}}}\left( {\xi ,\eta ,z} \right),{\boldsymbol{\tilde H}}\left( {\xi ,\eta ,z} \right)} \right]} {{\rm{e}}^{{\rm{i}}\xi x + {\rm{i}}\eta y}}{\rm{d}}\xi {\rm{d}}\eta } \end{split} (42) 2. 双重无限域振荡函数的快速积分方法
三轴各向异性介质中,电磁场求解的另一难点在于快速、精确地实现双重傅里叶逆变换。为此,给出了2种双重无穷域振荡函数的高效积分方法。
2.1 双重高斯积分方法
低角度井中,积分核函数快速收敛,故采用高斯–勒让德积分方法。首先,将两重无限积分转化为半无穷积分与定积分的组合。以电场为例,令
\xi = {k_\rho }\cos\phi ,\eta = {k_\rho }\sin\phi ,式(42)变为:{\boldsymbol{E}}\left( {\boldsymbol{r}} \right) = \frac{1}{{{{\left( {2{\text{π}} } \right)}^2}}}\int_0^{2{\text{π}} } {\int_0^\infty {\widetilde {\boldsymbol{E}}\left( {{k_\rho },\phi ,z} \right)} {{\rm{e}}^{{\rm{i}}{k_\rho }\left( {x\cos\phi + y\sin\phi } \right)}}{k_\rho }{\rm{d}}\phi {\rm{d}}} {k_\rho } (43) 假定定积分[0,2π]内设置P个高斯积分点,式(43)写为:
{\boldsymbol{E}}\left( {\boldsymbol{r}} \right) \approx \frac{1}{{4{{\text{π}} ^2}}}\sum\limits_{k = 1}^P {{A_k}\int_0^\infty {\widetilde {\boldsymbol{E}}\left( {{k_\rho },{\phi _k},z} \right)} {{\rm{e}}^{{\rm{i}}{k_\rho }\left( {x\cos{\phi _k} + y\sin{\phi _k}} \right)}}{k_\rho }{\rm{d}}{k_\rho }} (44) 式中:Ak和ϕk分别为求积系数和高斯节点。
然后,选择合适的积分上限,实现空间域电磁场的精确计算,式(44)变为:
{\boldsymbol{E}}\left( {\boldsymbol{r}} \right) \approx \frac{1}{{4{{\text{π}} ^2}}}\sum\limits_{k = 1}^P {\sum\limits_{l = 1}^P {\left( {{A_k}{B_l}\widetilde {\boldsymbol{E}}\left( {{k_{l,}},{\phi _k},z} \right){{\rm{e}}^{{\rm{i}}{k_l}\left( {x\cos{\phi _k} + y\sin{\phi _k}} \right)}}{k_l}} \right)} } (45) 式中:Bl和kl为kρ方向的求积系数和高斯节点。
2.2 双重正余弦滤波积分方法
大斜度井中,积分核函数振荡性强、衰减慢,不适合用直接积分方法。为此,采用正余弦滤波系数方法。以第二重积分为例,令
\widetilde {\boldsymbol{e}}\left( {\xi ,z} \right) = \displaystyle \int_{ - \infty }^\infty {\widetilde {\boldsymbol{E}}\left( {\xi ,\eta ,z} \right)} {{\rm{e}}^{{\rm{i}}\eta y}}\cdot {\rm{d}}\eta ,有:\begin{split} \widetilde {\boldsymbol{e}}\left( {\xi ,z} \right) =& \int_0^\infty \left[ \left( {\widetilde {\boldsymbol{E}}\left( {\xi ,\eta ,z} \right) + \widetilde {\boldsymbol{E}}\left( {\xi , - \eta ,z} \right)} \right)\cos\left( {\eta y} \right) +\right.\\ & \left.{\rm{i}}\left( {\widetilde {\boldsymbol{E}}\left( {\xi ,\eta ,z} \right) - \widetilde {\boldsymbol{E}}\left( {\xi , - \eta ,z} \right)} \right)\sin\left( {\eta y} \right) \right] {\rm{d}}\eta \end{split} (46) 基于数值滤波理论,结合离散采样方法[35],式(42)变为:
\begin{split} \widetilde {\boldsymbol{e}}\left( {\xi ,z} \right) =& \frac{1}{y}\sum\limits_{l = 1}^M \left[ \widetilde {\boldsymbol{E}}\left( {\xi ,\frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} + {\rm{i}}W_l^{\rm{s}}} \right) +\right.\\ &\left.\widetilde {\boldsymbol{E}}\left( {\xi , - \frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} - {\rm{i}}W_l^{\rm{s}}} \right) \right] \end{split} (47) 式中:
W_l^{\rm{c}} 和W_l^{\rm{s}} 分别为正弦和余弦滤波的权重;M为滤波点个数;{\lambda _l} 为采样点位置。利用同样的方法,对ξ方向的无穷积分进行处理,即可得到空间中电场的双重数值滤波积分表达式:
\begin{split} {\boldsymbol{E}}\left( {\boldsymbol{r}} \right) = &\frac{1}{{{{\left( {2{\text{π}} } \right)}^2}}}\frac{1}{{xy}}\sum\limits_{p = 1}^M {\sum\limits_{l = 1}^M {\left[ \widetilde {\boldsymbol{E}}\left( {\frac{{{\lambda _p}}}{x},\frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} + iW_l^{\rm{s}}} \right) + \widetilde {\boldsymbol{E}}\left( {\frac{{{\lambda _p}}}{x}, - \frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} - iW_l^{\rm{s}}} \right) \right]\left( {W_p^{\rm{c}} + iW_p^{\rm{s}}} \right)} } +\\& \frac{1}{{{{\left( {2{\text{π}} } \right)}^2}}}\frac{1}{{xy}}\sum\limits_{p = 1}^M {\sum\limits_{l = 1}^M {\left[ {\left( \widetilde {\boldsymbol{E}}\left( { - \frac{{{\lambda _p}}}{x},\frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} + iW_l^{\rm{s}}} \right) + \widetilde {\boldsymbol{E}}\left( { - \frac{{{\lambda _p}}}{x}, - \frac{{{\lambda _l}}}{y},z} \right)\left( {W_l^{\rm{c}} - iW_l^{\rm{s}}} \right) \right)} \right]} \left( {W_p^{\rm{c}} - iW_p^{\rm{s}}} \right)} \end{split} (48) 建立无限厚各向异性地层模型,其中沿x,y和z方向的电阻率分别为1,1和4 Ω·m。假定源工作频率为10 kHz,接收点位为(0.5, 0.05, 0.01)(见图2(a))。以测量的磁场z向分量为例,双重高斯积分核函数在kρ-ϕ平面的展布如图2(b)所示。
从图2(b)可以看出,核函数沿kρ和ϕ轴均具有强烈的振荡性,且kρ大于75时衰减仍未结束。当采用双重滤波系数方法且滤波点数量为241时,积分核函数沿λp和λl的分布如图2(c)所示。与图2(b)相比,经数值滤波处理后,可以忽略核函数振荡性,极大地降低了积分难度。进一步对比2种方法的采样点数量与积分精度,结果见表1。由表1可以看出,数值滤波方法仅需81个采样点即可保证计算精度,且计算速度在高斯积分方法的200倍以上。
表 1 2种积分方法精度与速度的对比Table 1. Accuracy and speed comparison of two integral methods方法 采样点数量 误差,% 相对耗时 高斯积分 1 200 0.004 0 219.50 2 000 0.000 3 610.00 数值滤波 81 0.006 0 1.00 241 0.002 5 8.85 601 0.001 6 55.10 2.3 算法准确性验证
为进一步验证递推算法的准确性,建立2层三轴各向异性地层模型,第一层的电导率分别为1.00,0.50和0.25 S/m,第二层的电导率分别为0.10,0.05和0.025 S/m。采用三轴发射和接收,源距和频率分别为1 m和10 kHz。仪器轴与地层法向夹角θ分别为0°,45°和80°时,接收位置处的部分磁场强度分量如图3所示。从图3可以看出,2种方法计算结果的一致性高,表明递推算法准确、可靠,且能在任意角度下精确模拟磁偶极子源激发的电磁场。
3. 实例分析
3.1 电缆测井中的应用
针对直井电缆测井,研究不同方向各向异性条件下多分量感应测井响应特征。图4(a)所示为典型的多分量感应测井仪器的基本结构。其中,每个发射/接收天线系由3个相互正交的磁偶极子线圈组成。考虑2种各向异性地层,即VTI各向异性和HTI各向异性,分别见图4(b)和图4(c)。VTI介质中水平方向的电导率相等为σ1,垂直方向电导率为σ2。假定HTI介质xoz平面为电导率为σ1的各向同性平面,y方向电导率为σ2。
直井中多分量感应测井仅同轴和共面分量不为零,故仅考虑测量σxx,σyy和σzz。仪器工作频率为10 kHz,源距为0.635 m,对于VTI介质,多分量感应测井在俄克拉荷马模型中的响应特征如图5所示。从图5可以看出:zz分量很好地反映了地层水平方向的电导率,即σ1;σxx和σyy的响应完全一致,其受地层垂向电导率的影响更大。
对于HTI介质,仪器工作频率为10 kHz,源距为0.635 m,多分量感应测井响应特征如图6所示,可以看出xx和yy分量不再重合,表明两者受y向各向异性的影响不同。与图5相比,同轴分量电导率变低,其受σ2的影响增大,此时若采用zz分量评价地层电导率σ1,则存在较大的误差,从而影响储层饱和度等参数计算的准确性。
3.2 随钻地质导向中的应用
随钻方位电磁波测井仪(ARM)具有地层界面实时方位探测和地质导向能力,其基础在于测量磁场交叉分量(Hzx),即地质信号。为对比各向异性对地质信号的影响,建立3层地层模型,其中中间层为电导率为1.0 S/m的各向同性地层,两侧围岩为电导率一致的高阻各向异性地层。考虑各向同性、VTI、HTI和三轴各向异性等4种围岩情况建立地层模型,围岩的电阻率见表2。测井仪器频率和源距分别为100 kHz、2 m时,4种模型地质信号响应如图7所示。从图7可以看出,各向同性介质中测井仪器靠近层界面时,ARM幅度单调增加。利用非零幅度特征,可定量计算仪器距地层界面的距离(DTB)。同时,仪器自高阻地层进入低阻地层和自低阻进入高阻地层时,ARM测量信号符号相反。基于此特性,可以准确判断地层界面的上、下方位信息。
表 2 不同地层模型围岩的电导率Table 2. Conductivities of surrounding rock in different formation models地层模型 σx/(S·m−1) σy/(S·m−1) σz/(S·m−1) 各向同性 0.25 0.25 0.25 VTI 0.25 0.25 0.10 HTI 0.25 0.10 0.25 各向异性 0.20 0.10 0.05 受垂向各向异性的影响,ARM即使在无限厚VTI和三轴各向异性介质中的响应仍不为零,此时非零响应将导致邻近存在地层界面的假象,影响DTB和界面方位评价的准确性。相比而言,无限厚HTI地层中ARM响应基本不受y方向各向异性的影响;然而,测井仪器靠近界面时,ARM在HTI地层和各向同性地层中的响应曲线分离严重,导致DTB计算结果出现误差。
4. 结 论
1)幅度传播矩阵递推算法能够准确、稳定计算层状三轴各向异性介质谱域电磁场分布,克服了传统系数传播矩阵存在的数值溢出问题。该递推算法可以模拟多分量感应测井和随钻方位电磁波测井在VTI、HTI和三轴各向异性介质中的响应规律。
2)双重高斯积分适用于井斜角小于45°的低角度井,双重数值滤波系数方法可以解决大斜度井积分的问题。采样点为81个时,数值滤波系数方法的计算速度比传统方法快2个数量级以上。
3)不同方向电导率的各向异性对多分量感应测井响应的影响严重,传统基于VTI模型的电导率解释模型不再适用。对随钻方位电磁波测井,横向和纵向各向异性则可能引起地层界面的假象,导致DTB和界面方位评价的准确性降低。对页岩、致密砂岩等复杂油气藏电磁波测井资料进行定性解释和定量处理时,必须构建合适的各向异性模型。
-
表 1 2种积分方法精度与速度的对比
Table 1 Accuracy and speed comparison of two integral methods
方法 采样点数量 误差,% 相对耗时 高斯积分 1 200 0.004 0 219.50 2 000 0.000 3 610.00 数值滤波 81 0.006 0 1.00 241 0.002 5 8.85 601 0.001 6 55.10 表 2 不同地层模型围岩的电导率
Table 2 Conductivities of surrounding rock in different formation models
地层模型 σx/(S·m−1) σy/(S·m−1) σz/(S·m−1) 各向同性 0.25 0.25 0.25 VTI 0.25 0.25 0.10 HTI 0.25 0.10 0.25 各向异性 0.20 0.10 0.05 -
[1] 李宁,冯周,武宏亮,等. 中国陆相页岩油测井评价技术方法新进展[J]. 石油学报,2023,44(1):28–44. LI Ning, FENG Zhou, WU Hongliang, et al. New advances in methods and technologies for well logging evaluation of continental shale oil in China[J]. Acta Petrolei Sinica, 2023, 44(1): 28–44.
[2] 李潮流,武宏亮,李霞,等. 倾斜电各向异性地层阵列感应测井数据处理新方法[J]. 石油学报,2022,43(10):1439–1449. LI Chaoliu, WU Hongliang, LI Xia, et al. A new data processing method of array induction logging for inclined electrical anisotropic formation[J]. Acta Petrolei Sinica, 2022, 43(10): 1439–1449.
[3] 蔺敬旗,孟鑫,李晴晴,等. 砾岩储层电成像测井表征方法及应用:以准噶尔盆地玛湖凹陷砾岩油藏为例[J]. 石油钻探技术,2022,50(2):126–131. LIN Jingqi, MENG Xin, LI Qingqing, et al. Characterization method and application of electrical imaging logging in conglomerate reservoir: A case study in Mahu Sag of Junggar Basin[J]. Petroleum Drilling Techniques, 2022, 50(2): 126–131.
[4] 杨小兵,姚梦麟,王思静,等. 页岩气测井地质工程技术新需求及解决方案[J]. 天然气工业,2022,42(2):20–27. doi: 10.3787/j.issn.1000-0976.2022.02.003 YANG Xiaobing, YAO Menglin, WANG Sijing, et al. Shale gas logging, geology and engineering technologies: New requirements and solutions[J]. Natural Gas Industry, 2022, 42(2): 20–27. doi: 10.3787/j.issn.1000-0976.2022.02.003
[5] 张国艳,肖加奇,郝永杰. 三维感应测井数值计算与理论分析[J]. 测井技术,2012,36(1):15–19. ZHANG Guoyan, XIAO Jiaqi, HAO Yongjie. Numerical computation and theoretical analysis of three-dimensional induction logging tool[J]. Well Logging Technology, 2012, 36(1): 15–19.
[6] 于蕾,汪宏年,王浩森,等. 层状各向异性地层中含环状天线槽随钻方位电磁波测井几何因子算法[J]. 地球物理学报,2023,66(3):1298–1314. YU Lei, WANG Hongnian, WANG Haosen, et al. Algorithm of geometrical factor of the azimuthal electromagnetic wave logging while drilling with annular antenna recesses in layered anisotropic formation[J]. Chinese Journal of Geophysics, 2023, 66(3): 1298–1314.
[7] WANG Lei, QIAO Ping, LI Zhiqiang, et al. A new semianalytical algorithm for rapid simulation of triaxial electromagnetic logging responses in multilayered biaxial anisotropic formations[J]. Geophysics, 2023, 88(2): D115–D129. doi: 10.1190/geo2021-0714.1
[8] 李星翰,张文秀,陈鹏,等. 随钻方位电磁波电阻率测井数据精度分析及实时反演成像[J]. 地球物理学报,2023,66(9):3990–3998. LI Xinghan, ZHANG Wenxiu, CHEN Peng, et al. Accuracy analysis and real-time inversion and imaging logging while drilling azimuth electromagnetic wave resistivity measurements[J]. Chinese Journal of Geophysics, 2023, 66(9): 3990–3998.
[9] 岳喜洲,马明学,李国玉,等. 随钻方位电磁波电阻率测井技术与地质导向应用[J]. 测井技术,2021,45(2):122–127. YUE Xizhou, MA Mingxue, LI Guoyu, et al. Azimuthal electromagnetic wave resistivity logging while drilling technology and its application in geo-steering[J]. Well Logging Technology, 2021, 45(2): 122–127.
[10] 王磊,范宜仁,操应长,等. 大斜度井/水平井随钻方位电磁波测井资料实时反演方法[J]. 地球物理学报,2020,63(4):1715–1724. WANG Lei, FAN Yiren, CAO Yingchang, et al. Real-time inversion of logging-while-drilling azimuthal electromagnetic measurements acquired in high-angle and horizontal wells[J]. Chinese Journal of Geophysics, 2020, 63(4): 1715–1724.
[11] ANDERSON B I, BARBER T D, HABASHY T M. Interpretation and inversion of fully triaxial induction data; a sensitivity study[R]. SPWLA-2002-O, 2002.
[12] 邓少贵,刘天淋,王磊,等. 双轴各向异性介质多分量感应测井响应快速计算[J]. 地球物理学报,2020,63(1):362–373. DENG Shaogui, LIU Tianlin, WANG Lei, et al. Analytical solution of multicomponent induction logging response in biaxial anisotropic medium[J]. Chinese Journal of Geophysics, 2020, 63(1): 362–373.
[13] DAVYDYCHEVA S, KAMINSKY A. Triaxial induction logging: New interpretation method for biaxial anisotropic formations: part 1[J]. Interpretation, 2016, 4(2): SF151–SF164. doi: 10.1190/INT-2015-0136.1
[14] 高永德,王世越,常波涛,等. 基于随钻前视探测技术的异常高压气层综合识别方法[J]. 天然气工业,2022,42(10):98–106. GAO Yongde, WANG Shiyue, CHANG Botao, et al. An integrated identification approach of abnormally high-pressure gas zones based on the look-ahead while drilling technology[J]. Natural Gas Industry, 2022, 42(10): 98–106.
[15] 王磊,刘英明,王才志,等. 水平井随钻电磁波测井实时正反演方法[J]. 石油勘探与开发,2021,48(1):139–147. doi: 10.11698/PED.2021.01.12 WANG Lei, LIU Yingming, WANG Caizhi, et al. Real-time forward modeling and inversion of logging-while-drilling electromagnetic measurements in horizontal wells[J]. Petroleum Exploration and Development, 2021, 48(1): 139–147. doi: 10.11698/PED.2021.01.12
[16] 杨震,肖红兵,张智勇. 随钻方位电磁波电阻率仪器性能指标检测方法[J]. 测井技术,2020,44(5):448–452. doi: 10.16489/j.issn.1004-1338.2020.05.005 YANG Zhen, XIAO Hongbing, ZHANG Zhiyong. Main specifications test method of azimuthal electromagnetic logging while drilling tool[J]. Well Logging Technology, 2020, 44(5): 448–452. doi: 10.16489/j.issn.1004-1338.2020.05.005
[17] 刘天淋,岳喜洲,李国玉,等. 超深探测随钻电磁波测井地质信号特性研究[J]. 石油钻探技术,2022,50(6):41–48. LIU Tianlin, YUE Xizhou, LI Guoyu, et al. Study over the geo-signal properties of ultra-deep electromagnetic wave logging while drilling[J]. Petroleum Drilling Techniques, 2022, 50(6): 41–48.
[18] CHEN Y P, CHEW W C, JIANG Lijun. A new green’s function formulation for modeling homogeneous objects in layered medium[J]. IEEE Transactions on Antennas and Propagation, 2012, 60(10): 4766–4776. doi: 10.1109/TAP.2012.2207332
[19] HONG Decheng, XIAO Jiaqi, ZHANG Guoyan, et al. Characteristics of the sum of cross-components of triaxial induction logging tool in layered anisotropic formation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3107–3115. doi: 10.1109/TGRS.2013.2269714
[20] ANDERSON W L. Computer program; numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering[J]. Geophysics, 1979, 44(7): 1287–1305. doi: 10.1190/1.1441007
[21] LØSETH L O, URSIN B. Electromagnetic fields in planarly layered anisotropic media[J]. Geophysical Journal International, 2007, 170(1): 44–80. doi: 10.1111/j.1365-246X.2007.03390.x
[22] HU Yunyun, SUN Qingtao. Modeling of triaxial induction logging responses in multilayered anisotropic formations[J]. Geophysics, 2021, 86(4): E305–E314. doi: 10.1190/geo2020-0475.1
-
期刊类型引用(5)
1. 李新勇,李骁,赵兵,王琨,苟波. 顺北油田S井超深超高温碳酸盐岩断溶体油藏大型酸压关键技术. 石油钻探技术. 2022(02): 92-98 . 本站查看
2. 邱春阳,张翔宇,赵红香,王雪晨,张海青,陈二丁. 顺北区块深层井壁稳定钻井液技术. 天然气勘探与开发. 2021(02): 81-86 . 百度学术
3. 汤明光,刘清华,薛国庆,方小宇,唐慧敏. 海上低渗油藏大井距精细注水技术应用实践. 承德石油高等专科学校学报. 2020(02): 11-15+37 . 百度学术
4. 刘阳. 高温深层碳酸盐岩裸眼酸压完井封隔器研制与现场试验. 石油钻探技术. 2020(06): 76-81 . 本站查看
5. 刘彪,潘丽娟,张俊,白彬珍,李双贵. 顺北区块超深小井眼水平井优快钻井技术. 石油钻探技术. 2016(06): 11-16 . 本站查看
其他类型引用(0)