Tectonic characteristics and extensional-compressional rates of the North-east Sag of South Yellow Sea Basin
-
摘要: 通过选取南黄海盆地东北凹典型地震剖面，开展精细的构造解释，系统梳理了东北凹构造样式特征。采用平衡剖面恢复技术和伸缩率计算方法，恢复了东北凹各时期的地质演化剖面，分析了东北凹不同构造演化阶段的伸缩率变化特征。研究表明，南黄海盆地东北凹主要发育伸展构造、走滑构造（负花状）和反转构造等多种构造组合样式，经历了晚侏罗世的仪征运动和渐新世末的三垛运动，相应地在中—上侏罗统和渐新统沉积时期，东北凹处于明显的收缩阶段，伴随发育TK40和T20不整合界面。同时，本文结合区域应力场特征，探讨了南黄海盆地东北凹的构造演化历程：以两次构造运动为界，划分为3个构造演化阶段（晚三叠世—侏罗纪的初始断陷阶段、白垩纪—渐新世的裂陷-反转阶段、新近纪—第四纪的区域沉降阶段）。南黄海盆地东北凹伸缩率的时空变化及构造演化过程，是对“晚中生代以来，古太平洋板块相对欧亚板块俯冲汇聚速率和方向的改变”的局部响应。Abstract: Typical seismic sections are selected for delicate structural interpretation upon the basis of structural features of the North-east Sag of the South Yellow Sea Basin. The balanced cross-section technique and the extensional - compressional rates calculation method are used to reconstruct the tectonic evolution cross-sections for different tectonic evolutionary stages. The characteristics of extensional - compressional rates are analyzed. The results suggest that the North-east Sag of the South Yellow Sea basin is characterized by extensional structures, strike slip structures (negative flower patterns), inversion structures etc. by the Yizheng movement in Late Jurassic and the Sanduo movement by the end of Oligocene. Accordingly, the sag was obviously in compressional status during the deposition of Late Jurassic and Oligocene when the Tk40 and T20 unconformity interfaces were formed. In order to quantify the characteristics of the regional stress field, the paper also reconstructed the tectonic evolution of the north-east sag. Taking the two tectonic movements as the boundaries, the sag is subdivided into three stages of tectonic evolution i.e. the initial faulting- depression stage from Late Triassic to Jurassic, the rifting inversion stage from Cretaceous to Oligocene, and the regional subsidence stage from Neogene to Quaternary. The temporal and spatial variations in the extensional - compressional rates and the tectonic evolution of north-east sag in South Yellow Sea basin represent the local response to the regional stress field, such as the change in convergence rate and motion azimuth between the Pacific plate and the Eurasian plate since Late Mesozoic.
-
随着海洋地震勘探的发展，勘探目标逐渐趋于复杂化，高精度成像方法已经成为海洋油气勘探的瓶颈技术。各向异性现象普遍存在于地下介质中，基于各向同性的传统偏移方法难以满足实际生产中高精度成像的需求，会严重影响最终偏移结果及后续的偏移速度分析精度[1]。由于海洋拖缆地震勘探观测系统的特殊性，为了确保偏移精度，对偏移方法需提出更高要求。因此，有必要发展一种适应于海洋观测系统的各向异性高精度偏移成像方法，进一步提高地震资料解释的精度，解决海洋复杂油气藏精细勘探的问题。
近年来在三维实际资料处理中，以射线理论为基础的高斯束偏移方法取得了较好的应用效果，已成为各大油田公司关注的焦点，同时也是众多地球物理学家研究的热点之一。一方面，高斯束方法提供了地震波在地下传播的物理解释，更易在复杂的地球物理模型中描述波动现象，且由于射线追踪角度无限制，适合陡倾角成像问题；另一方面，与波动方程偏移相比，高斯束方法较为高效，非常适合三维实际资料成像处理。目前，工业界常用的Kirchhoff偏移也具有很高的成像效率，但其用于复杂探区成像处理时，射线焦散现象[2-4]及多路径问题[5-8]难以避免，从而影响成像精度。而高斯束偏移采用复值初始束参数动力学射线追踪, 能避免焦散点附近振幅奇异问题，且不同方向的局部平面波分解使其能够独立地从不同方向传播到地下成像点，不会产生多路径问题。总之，高斯束偏移不仅具有接近波动方程偏移的精度，而且具有Kirchhoff偏移灵活高效的特点。
高斯束方法起源于电磁学，在20世纪80年代初期，Červený和Popov将其引入到地球物理学领域，并取得了飞速发展[9-13]。Hill[14]基于高斯束表征的格林函数推导了叠后深度偏移公式，奠定了高斯束偏移的理论框架。其成像思想为：基于局部倾斜叠加近似，将分解后的平面波分量沿着不同方向以半速度反传到地下成像点处进行成像。后来，考虑到地下介质的复杂性，各向异性介质中的高斯束偏移方法的相关研究不断涌现，Alkhalifah[15]最早通过引入各向异性射线追踪方程发展了各向异性高斯束叠后偏移方法。Zhu等[16]通过引入相速度和群速度简化的各向异性射线追踪体系，提出了共炮域各向异性高斯束叠前深度偏移方法，进一步降低了基于弹性参数算法的计算量。段鹏飞等[17]在局部角度域提出了一种相对经济的动力学射线追踪近似方案。Protasov[18]提出了各向异性介质中针对多分量数据的高斯束偏移方法。张凯等[19]利用各向异性介质运动学和动力学射线追踪体系，提出了角度域各向异性高斯束逆时偏移方法。李振春等[20]基于转换波射线追踪算法实现了各向异性介质角度域转换波高斯束偏移方法。随着勘探精度要求逐步提高，常规的高斯束偏移难以满足高精度成像的要求，因此许多新的方法及优化算法相继出现。一方面，基于波动方程偏移中逆时偏移波场延拓思想的高斯束逆时偏移[19,21-23]极大地提高了高斯束偏移在复杂构造中的成像精度，有利于地下复杂构造中层位及断层识别等解释工作；另一方面，基于振幅保真理论[24-26]及最小二乘理论[27-29]的高斯束偏移在保幅成像及提高成像分辨率方面取得重大突破，极大地提高了地震解释过程中的AVA分析及AVO反演精度。此外，黄建平等[30]通过修改高斯束局部倾斜叠加提出了双复杂条件下非倾斜叠加精确束偏移方法。张瑞等[31]结合高斯束传播过程中反射张角信息的计算，发展了高斯束高精度叠前深度偏移方法，并在复杂海洋模型中取得了较好的应用效果。Yang和Zhu[32]在传统高斯束偏移基础上基于菲涅尔带思想提出了一种优化的高斯束偏移方法，提高了稀疏数据的偏移成像精度。
为了弥补高斯束偏移在处理低信噪比资料上的不足，出现了一种改进的高斯束偏移方法—控制束偏移方法。Vinje等[33]提出了各向同性介质的控制束偏移方法，并将其应用到裂缝基底成像、盐丘底边界成像、低覆盖次数及宽方位角数据成像等方面。Zhou等[34]将各向同性介质控制束偏移推广到各向异性介质中，发展了一种各向异性高保真度控制束偏移方法。Casasanta等[35]实现了基于共偏移距的PS转换波高斯束偏移。黄建平等[36]认为可以在τ-p域内设定阈值，针对不同信噪比的实际数据进行适当的切除来实现控制束偏移。
本文在传统高斯束偏移的基础上，考虑到地下介质的各向异性和实际资料中含有较强噪声干扰，根据有效信号和干扰信号在τ-p域中的相干性差异，基于各向异性运动学和动力学射线追踪，发展了一种VTI介质共偏移距域数据驱动控制束偏移理论方法。在高斯束偏移成像过程中，计算τ-p道集的相似系数，通过设定相似系数阈值对干扰信号进行控制，从而降低偏移剖面中的随机噪声，提高同相轴的连续性。通过各向异性洼陷模型、修改的SEG/Hess VTI模型对新方法进行测试，并与传统高斯束偏移结果进行对比验证本文方法的正确性及有效性；最后通过对中国东部某海上探区的实际资料进行处理，探讨了本文方法的优势及其对海上数据的适应性。
1. 方法原理
1.1 各向异性射线追踪
各向异性高斯束偏移方法的核心在于利用各向异性运动学和动力学射线追踪求解高斯束的复值走时和振幅。Červený[37]由各向异性介质无源弹性波波动方程经过一系列数学推导，最后得到各向异性运动学射线追踪方程组
$$\left\{ {\begin{aligned} &{\frac{{{\rm{d}}{x_i}}}{{{\rm{d}}\tau }} = {a_{ijkl}}{p_l}{g_j}{g_k}}\\ &{\frac{{{\rm{d}}{p_i}}}{{{\rm{d}}\tau }} = - \frac{1}{2}\frac{{\partial {a_{njkl}}}}{{\partial {x_i}}}{p_n}{p_l}{g_j}{g_k}} \end{aligned}} \right.$$ (1) 其中，aijkl=cijkl/ρ表示密度归一化的弹性参数；τ表示沿中心射线路径的旅行时；pi=∂τ/∂xi为慢度矢量；gi为极化矢量。
在二维情况下，式（1）中特征向量可以表示为
$$\left\{ {\begin{aligned} & {{g_1}{g_1} = \frac{{{\Gamma _{33}} - G}}{{{\Gamma _{11}} + {\Gamma _{33}} - 2G}}} \\ &{{g_3}{g_3} = \frac{{{\Gamma _{11}} - G}}{{{\Gamma _{11}} + {\Gamma _{33}} - 2G}}} \\ & {{g_1}{g_3} = \frac{{ - {\Gamma _{13}}}}{{{\Gamma _{11}} + {\Gamma _{33}} - 2G}}} \end{aligned}} \right.$$ (2) 式中，Γik=aijklpjpl为Christoffel矩阵，G为程函方程，当其取不同特征值时代表不同的程函方程。
根据式（1）即可求得中心射线路径从而构建射线中心坐标系，在各向异性介质中射线中心坐标系不再正交，因此需要在求解过程中引入一个沿中心射线路径的权值变量以校正这种非正交性所带来的误差。基于各向同性介质动力学射线追踪方程组，Hanyga[38]推导出了各向异性介质中的动力学射线追踪方程组
$$\left\{ {\begin{aligned} &{\frac{{{\rm{d}}Q}}{{{\rm{d}}\tau }} = MP + VQ}\\ &{\frac{{{\rm{d}}P}}{{{\rm{d}}\tau }} = - VP - HQ} \end{aligned}} \right.$$ (3) 其中，P和Q是动力学射线追踪参量，M、V、H表示程函方程对n和pn的导数，满足下列关系式
$$\left\{ {\begin{aligned} &{H = \frac{1}{2}\frac{{{\partial ^2}{G_m}}}{{\partial p_n^2}} - \frac{1}{4}{{\left( {\frac{{\partial {G_m}}}{{\partial {p_n}}}} \right)}^2}} \\ & {M = \frac{1}{2}\frac{{{\partial ^2}{G_m}}}{{\partial {n^2}}} - \frac{1}{4}{{\left( {\frac{{\partial {G_m}}}{{\partial n}}} \right)}^2}} \\ & {V = \frac{1}{2}\frac{{{\partial ^2}{G_m}}}{{\partial {p_n}\partial n}} - \frac{1}{4}\frac{{\partial {G_m}}}{{\partial {p_n}}}\frac{{\partial {G_m}}}{{\partial n}}} \end{aligned}} \right.$$ (4) 上式中的Gm=aijklpiplgjgk为Christoffel方程中行列式det（Γjk−Gmδjk）=0的特征值，代表各向异性介质中存在的三种不同地震波的程函方程，当m取值1、2和3时分别表示qP波、qSV波和qSH波。式（1）和式（3）中的各向异性介质运动学和动力学射线追踪方程均以弹性参数表示，目前人们习惯使用物理意义更易理解的Thomsen参数来表征弱各向异性介质，同时可以通过Thomsen参数与弹性参数转化关系来消除这种不便。
1.2 数据驱动控制束偏移
叠前高斯束偏移是通过高斯束构建炮点和检波点在地下成像点处的格林函数，然后通过它们来构建上行波场和下行波场，最后利用上、下行波场的互相关得到偏移结果。本文从频率域的互相关成像条件出发来推导出数据驱动控制束偏移公式，根据Zhang等[39]的文章可知，二维频率域互相关成像条件可以表示为
$$I({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}}) = - i\int {{\rm{d}}\omega {\rm{sgn}}(\omega ){P_{{\rm{up}}}}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega )P_{{\rm{down}}}^*} ({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega )$$ (5) 式中，Pup和Pdown分别表示上、下行波场，*表示复共轭。Gray和Bleistein[24]给出了其表达式如下：
$$\begin{aligned} {P_{{\rm{up}}}}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega ) =& 2\int {{\rm{d}}{x_{\rm{r}}}{P_{{\rm{up}}}}({\bf{x}_{\bf{g}}},{\bf{x}_{\bf{s}}},\omega )\frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}}\times\\&{\sqrt { - i\omega } {A_{{\rm{ART,g}}}}{\rm{exp}}[ - i\omega {\tau _{\rm{g}}}]} \\=& - 2\int {{\rm{d}}{x_{\rm{g}}}{P_{{\rm{up}}}}({\bf{x}_{\bf{g}}},{\bf{x}_{\bf{s}}},\omega )\frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}}\times\\&i\omega {G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{g}}},\omega ) \end{aligned}$$ (6) $$\begin{aligned} P_{{\rm{down}}}^*({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega ) = & - 2\frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\sqrt { - i\omega } {A_{{\rm{ART,s}}}}{\rm{exp}}[ - i\omega {\tau _{\rm{s}}}]\\ = &2\frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}i\omega {G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega ) \end{aligned}$$ (7) 式中，Vs和Vg分别表示炮点xs和检波点xg处的速度，θs和θg分别表示炮点和检波点处的出射角，τs、τg和AART,s、AART,g分别表示渐进射线理论下的实值走时和振幅，G（x0,x,ω）是波动方程中的格林函数，将式（6）和式（7）带入到式（5）中可得如下共炮域偏移成像公式：
$$\begin{aligned} I({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}}) = & - 4\int {{\rm{d}}\omega i\omega \left| \omega \right|\frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}{G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega )}\times \\ & \int {{\rm{d}}{x_{\rm{g}}}{P_{{\rm{up}}}}({\bf{x}_{\bf{g}}},{\bf{x}_{\bf{s}}},\omega )\frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}{G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{g}}},\omega )} \end{aligned}$$ (8) 在高斯束方法中，格林函数是通过一系列由震源点x出射的且具有不同出射方向的高斯束的叠加积分表示的：
$$G({\bf{x}_{\bf{0}}},{\bf{x}},\omega ) = \frac{i}{{4\pi }}\int {\frac{{{\rm{d}}{p_x}}}{{{p_z}}}} A({\bf{x}_{\bf{0}}},{\bf{x}}){\rm{exp}}[i\omega T({\bf{x}_{\bf{0}}},\bf{x})]$$ (9) 其中，px和pz分别为射线参数的水平和垂直分量，T（x0,x）和A（x0,x）分别是复值走时和振幅，它们可以通过式（1）和式（3）表示的各向异性射线追踪方程得到。对于海上接收的共偏移距地震记录，可以通过如下变换：
$${\bf{m}} = \frac{1}{2}({\bf{x}_{\bf{g}}} + {\bf{x}_{\bf{s}}}),{\bf{h}} = \frac{1}{2}({\bf{x}_{\bf{g}}} - {\bf{x}_{\bf{s}}})$$ (10) 将式（8）中的积分变量由炮点-检波点域变换到中心点-半偏移距域，
$$\begin{aligned}I({\bf{x}_{\bf{0}}},{\bf{h}}) =& - 4\int {{\rm{d}}\omega i\omega \left| \omega \right|\frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}{G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega )}\times\\&\int {{\rm{d}}{m_x}P({\bf{h}},{\bf{m}},\omega )\frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}{G^*}({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{g}}},\omega )} \end{aligned}$$ (11) 其中，P（h,m,ω）表示中心点坐标为m，半偏移距为h的共偏移距道集，这种类型的道集可以理解为半偏移距h固定，沿中心点坐标m变化的地震记录。接着，需要沿中心点坐标来确定中心点位置Lm，在确定Lm后，结合式（9）表示的格林函数，根据如下关系：
$${\bf{x}_{\bf{s}}} = {\bf{m}} - {\bf{h}},{\bf{x}_{\bf{g}}} = {\bf{m}} + {\bf{h}}$$ (12) 将中心点Lm有效范围内的炮点、检波点格林函数分别用Lm−h、Lm+h处出射的高斯束来近似表示，
$$\begin{aligned} G({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{s}}},\omega ) =& \frac{i}{{4\pi }}\int {\frac{{{\rm{d}}{p_{sx}}}}{{{p_{sz}}}}} A({\bf{x}_{\bf{0}}},{\bf{L}_{\bf{m}}} - {\bf{h}})\times\\&{\rm{exp}}[i\omega T({\bf{x}_{\bf{0}}},{\bf{L}_{\bf{m}}} - {\bf{h}})]\\ G({\bf{x}_{\bf{0}}},{\bf{x}_{\bf{g}}},\omega ) =& \frac{i}{{4\pi }}\int {\frac{{{\rm{d}}{p_{gx}}}}{{{p_{gz}}}}} A({\bf{x}_{\bf{0}}},{\bf{L}_{\bf{m}}} + {\bf{h}})\times\\&{\rm{exp}}[i\omega T({\bf{x}_{\bf{0}}},{\bf{L}_{\bf{m}}} + {\bf{h}})] \end{aligned}$$ (13) 其中，psx、psz和pgx、pgz分别为炮点和检波点射线参数的水平和垂直分量。由于高斯束具有一定的宽度，所以检波点处会接收到不同波束的能量贡献。为了区分这些波束，可以结合高斯束在初始位置波前面为平面的特点，采用局部倾斜叠加方法来将这些波束的能量反传到地下成像点。根据Hill[14]文章中给出的一维高斯窗函数：
$$\frac{{\Delta {x_L}}}{{\sqrt {2\pi } {w_0}}}\sqrt {\left| {\frac{\omega }{{{\omega _r}}}} \right|} \sum\limits_{{x_L}} {{\rm{exp}}\left[ { - \left| {\frac{\omega }{{{\omega _r}}}} \right|\frac{{{{\left( {{x_g} - {x_L}} \right)}^2}}}{{2w_0^2}}} \right]} \approx 1$$ (14) 式中，xL为高斯窗中心位置，w0是高斯窗宽度，ωr为参考频率，将式（13）和式（14）代入到公式（11）中可得
$$\begin{aligned} &I({\bf{x}_{\bf{0}}},{\bf{h}}) = \frac{{\Delta L}}{{4{\pi ^2}\sqrt {2\pi } {w_0}}}\sum\limits_{{x_L}} {\int {{\rm{d}}\omega i\omega \frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\int {\frac{{{\rm{d}}{p_{sx}}}}{{{p_{sz}}}}} A_s^*} }\times\\ & {\rm{exp}}\left[ { - i\omega T_{\rm{s}}^*} \right]\times \frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}\int {\frac{{{\rm{d}}{p_{gx}}}}{{{p_{gz}}}}} A_{\rm{g}}^*{\rm{exp}}\left[ { - i\omega T_{\rm{g}}^*} \right]{D_{\rm{S}}}({x_L},{p_m},\omega ) \end{aligned}$$ (15) 其中，pm=ps+pg表示中心点射线参数，DS（xL,pm,ω）表示共偏移距道集加窗局部倾斜叠加，
$$\begin{aligned}{D_{\rm{S}}}({x_L},{p_m},\omega ) =& {\left| {\frac{\omega }{{{\omega _r}}}} \right|^{3/2}}\int {{\rm{d}}{m_x}} P({\bf{h}};{\bf{m}};\omega )\times\\&{\rm{exp}}\left[ { - i\omega {p_m}({m_x} - {x_L}) - \left| {\frac{\omega }{{{\omega _r}}}} \right|\frac{{{{({m_x} - {x_L})}^2}}}{{2w_0^2}}} \right]\end{aligned}$$ (16) 式（15）表示的是先对频率、炮点出射的高斯束以及检波点处出射的高斯束进行三重积分，然后再对不同中心点所对应的三重积分求和的过程。考虑到高斯束中的振幅和走时均为复数，可以令
$$\begin{array}{l} {T_{{\rm{re}}}} = {\rm{Re}}(T_{\rm{s}}^* + T_{\rm{g}}^*),{A_{{\rm{re}}}} = {\rm{Re}}(A_{\rm{s}}^*A_{\rm{g}}^*)\\ {T_{{\rm{im}}}} = {\rm{Im}}(T_{\rm{s}}^* + T_{\rm{g}}^*),{A_{{\rm{im}}}} = {\rm{Im}}(A_{\rm{s}}^*A_{\rm{g}}^*) \end{array}$$ (17) 直接求解公式（15）需要面对一系列频率求和的巨大计算量，根据高成等[40]，可以通过改变积分次序，并且应用傅里叶反变换和希尔伯特变换，推导出式（15）对应的时间域表达式
$$\begin{aligned} & I({\bf{x}_{\bf{0}}},{\bf{h}}) =\frac{{\Delta L{\omega _r}}}{{4\pi \sqrt {2\pi } {w_0}}}\sum\limits_{{x_L}} {\int {\frac{{{\rm{d}}{p_{sx}}}}{{{p_{sz}}}}} \int {\frac{{{\rm{d}}{p_{gx}}}}{{{p_{gz}}}}} } \frac{{\cos {\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\frac{{\cos {\theta _{\rm{g}}}}}{{{V_{\rm{g}}}}}\times\\ &[{A_{{\rm{re}}}}{{\bar D}_{\rm{S}}}({x_L},{p_m},{T_{{\rm{re}}}},{T_{{\rm{im}}}}) - {A_{{\rm{im}}}}{{\bar D}_{{\rm{S}},{\rm{H}}}}({x_L},{p_m},{T_{{\rm{re}}}},{T_{{\rm{im}}}})] \end{aligned}$$ (18) 根据Hale[41]的研究可知，成像过程中可首先评估虚值走时Tim的范围，然后对Tre和Tim进行离散采样，并计算
${\bar D_{\rm{S}}}$ 和${\bar D_{{\rm{S,H}}}}$ 的离散样板表，最后根据地下成像点的旅行时信息与离散样板表进行匹配成像。该过程可以理解为一个数据驱动的过程，即从道集中提取每一样点采集“属性特征”和旅行时，然后找到这些特定数据点所指向的成像点。对含有噪声干扰的海上共偏移距道集而言，噪声会在一定程度上影响样点信息的正确提取，从而影响最终的偏移结果。通过对含噪声的共偏移距道集进行加窗局部倾斜叠加可以发现，τ-p坐标图版上的峰值显示在局部共偏移距道集中存在相干同相轴，并且这个同相轴来自地下的某一特定反射段，而一些相对较弱的能量团可能是一些背景和随机噪声等干扰信号。为了降低随机噪声对偏移结果产生的影响，本文采用相干体技术来提高偏移剖面的信噪比。相干体技术不仅能用于检测地震波同相轴的不连续性，也可以用来识别断层、特殊岩性体以及河道等，在地震资料解释中扮演着重要角色[42]。目前，相干体算法主要可以划分为三类：基于相关的算法、基于相似系数的算法以及基于本征结构分析的算法[43-45]。本文采用基于相似系数的算法，推导τ-p道集的相似系数计算公式。根据Neidell等[46]可知，相似系数可以定义为归一化的输出/输入能量比值，其中输出道能量表示为输入道的一个简单合成或者求和。在信号分析中相似系数的表达式
$${S_{\rm{C}}} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left(\sum\limits_{j = 1}^J {{f_{i,j}}} \right)}^2}} \bigg/\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^J {f_{i,j}^2} } $$ (19) 式中，N表示在一个短时窗内的样点数，J表示总道数，下标i和j分别表示样点序号和道号，fij为样点（i, j）处的振幅值。在此基础上，根据式（19）可以推导出τ-p道集的相似系数计算公式（20），同时根据有效信号和干扰信号在τ-p域中的相干性差异，可以应用相似系数阈值滤波公式（21）对离散样板表进行滤波处理，去除τ-p道集中的部分干扰能量，从而提高偏移剖面中同相轴的连续性。
$${S_{\rm{C}}}(\tau ,{x_L},{p_m}) = \frac{{\displaystyle\sum\nolimits_{i = 1}^N {\left[ {\displaystyle\sum\nolimits_{j = 1}^J {u(\tau + i\Delta t - {p_m}\Delta x,{x_j})} } \right]} }}{{N\displaystyle\sum\nolimits_{i = 1}^N {\displaystyle\sum\nolimits_{j = 1}^J {{u^2}(\tau + i\Delta t - {p_m}\Delta x,{x_j})} } }}$$ (20) $$W({\tau _i},{x_L},{p_j}) = \left\{ {\begin{array}{*{20}{c}} {1,{S_{\rm{C}}}{\text{≥}} {\rm{tp}}}\\ {0,\;{\rm{others}}\;} \end{array}} \right.$$ (21) 其中，u表示局部共偏移距道集，xj为高斯窗内第j道的笛卡尔坐标，xL为参考道的位置；τ表示走时，Δt和Δx分别表示时间采样间隔和空间采样间隔；tp为相似系数阈值参数，范围为0～1。
如图1所示，VTI介质共偏移距域数据驱动控制束偏移过程可以概括为如下4个步骤：（1）计算τ-p道集对应的相似系数并选取大于相似系数阈值参数对应的τ-p道集中的能量；（2）利用各向异性射线追踪分别计算来自炮点和检波点的中心射线以及旁轴射线；（3）组合两组旁轴射线来计算两个“胖射线”相交区域上所有点的走时和振幅；（4）根据旅行时及射线参数关系选取τ-p道集中的能量进行匹配得到偏移后的子波束，所有偏移后的子波束叠加得到最终的成像结果。单偏移距数据驱动控制束具体的实现流程如图2所示。
![]()
图 1 共偏移距域数据驱动控制束偏移原理图a. 偏移后的子波束，b. τ-p道集，c. 相似系数。Figure 1. Diagram of the data-driven CBM in common-offset domaina. migrated sub-beam，b. τ-p gather，c. semblance coefficient.![]()
图 2 单偏移距数据驱动控制束偏移程序流程图Figure 2. Flow chart of the single-offset data-driven controlled beam migration program2. 应用实例
从三个方面来进一步测试本文方法的正确性、有效性及实用性。首先，在编程实现算法的基础上，对简单的各向异性洼陷模型进行偏移处理，并与常规高斯束偏移结果对比，验证本文方法的正确性及其在处理存在各向异性地层中的优势；其次，对修改的SEG/Hess VIT复杂模型进行测试，验证本文方法的有效性及处理低信噪比数据的优势。最后，通过对中国东部某海上探区实际资料进行成像处理，验证本文方法对海上数据的适应性。
2.1 各向异性洼陷模型
为了验证本文方法的正确性，建立了如图3a所示的各向异性洼陷模型。正演模拟采用交错网格各向异性高阶有限差分正演方法：震源为主频20 Hz的Ricker子波；接收方式采用共偏移距接收方式，偏移距范围从−2 000 m到2 000 m，偏移距间隔为20 m，共201个共偏移距道集，每个共偏移距道集有401道，道间距为20 m，总记录时间为3 s，采样率为4 ms。图3b和3c分别显示了一个近偏移距（200 m）和远偏移距（1 500 m）道集。基于相同的成像参数，分别采用各向同性和各向异性共偏移距域高斯束偏移方法对其进行偏移成像，结果如图4所示。从图4a和4b对比可以发现，对于近偏移距数据，采用各向同性方法时，由于忽略了各向异性的影响，弯曲底界面成像位置不准确，曲率发生了明显的变化（黑色箭头所示）。对比图4c和4d的远偏移距结果可以看出，各向异性参数对共偏移距道集中的大偏移距信息成像影响较大，各向同性结果中弯曲底界面处存在明显的能量发散问题。图4e和4f为所有偏移距叠加后成像结果，采用各向异性偏移方法可以很好地聚焦成像（白色箭头所示），界面同相轴清晰且位置准确，界面弯曲曲率与实际模型更加吻合。
![]()
图 3 各向异性洼陷模型a. 各向异性参数场，b. 200 m偏移距道集，c. 1 500 m偏移距道集。Figure 3. Anisotropic sag modela. anisotropic parameters，b. common-offset gather （200 m），c. common-offset gather （1 500 m）.![]()
图 4 各向异性洼陷模型偏移结果a. 各向同性（偏移距为200 m），b. 各向异性（偏移距为200 m），c. 各向同性（偏移距为1 500 m），d. 各向异性（偏移距为1 500 m），e. 各向同性所有偏移距叠加，f. 各向异性所有偏移距叠加。Figure 4. Results of GBM on anisotropic sag modela. isotropic （offset=200 m），b. anisotropic （offset=200 m），c. isotropic （offset=1 500 m），d. anisotropic （offset=1 500 m）， e. all offset stack isotropic result，f. all offset stack anisotropic result.2.2 修改的SEG/Hess VTI模型
为了进一步验证本文方法的有效性，接着对不包含盐丘构造的部分SEG/Hess VTI模型进行测试，其各向异性参数如图5所示，可以看出，该模型具有典型的正断层、岩性低速体及各向异性透镜体，具有一定的构造代表性。共偏移距道集通过交错网格各向异性高阶有限差分正演模拟得到，震源为主频20 Hz的Ricker子波，为了得到完整的共偏移距道集，采取模型扩边处理，总共201个共偏移距道集，偏移距范围从−2 000 m到2 000 m，偏移距间隔为20 m，每个共偏移距道集有451道，道间距为20 m，总记录时间为3 s，采样率为4 ms。
基于相同的成像参数，分别采用各向同性和各向异性共偏移距域高斯束偏移对该模型进行处理，结果如图6a和6b所示。通过对比可以看出，各向异性相比各向同性而言，在断层附近成像位置更加准确且断点绕射能量收敛（图6红色箭头），地层的成像能量更加聚焦且同相轴更加清晰（图6蓝色箭头）。为了更加清楚地对比，将图6中绿色矩形框内成像结果进行放大显示（图6c和图6d）。从放大图中可以看出，由于各向异性高斯束偏移考虑地下各向异性的影响，可以大幅改善断点的成像效果（图6c和图6d蓝色箭头），并且图6c中透镜体未能很好成像，而各向异性的成像效果（图6d红色箭头）远好于各向同性偏移结果。
![]()
图 6 修改的SEG/Hess VIT模型共偏移距域高斯束偏移结果及局部放大显示对比a. 各向同性，b. 各向异性，c. 各向同性局部放大，d. 各向异性局部放大。Figure 6. Results of common-offset GBM and magnified views on modified SEG/Hess VTI modela. isotropic，b. anisotropic，c. isotropic magnified view，d. anisotropic magnified view.为了测试共偏移距域控制束偏移的成像效果，选取了偏移距为200 m的共偏移距道集（图7a）作为测试对象，并使用Seismic Unix软件包中的suaddnoise命令在原始道集中加入随机噪声得到低信噪比共偏移距道集（图7b）,加噪后的共偏移距道集中的同相轴变得模糊并混叠在背景噪声中。图8为图7中局部共偏移距道集（蓝色矩形框）相干性分析结果，通过对低信噪比局部共偏移距道集（图8a）进行倾斜叠加得到图8b所示的τ-p道集，与图8d所示的原始道集生成的标准τ-p道集对比可以发现，随机噪音在τ-p道集中表现为一些比较弱的能量，同时可以注意到即使是在标准τ-p道集能量也不能很好聚焦（图8d红色箭头），这主要是因为局部τ-p变换本身存在噪声泄露问题，然而这些能量却可以在控制束中得到压制（图8e、8f红色箭头）。随后通过计算图8b对应的相似系数，结果如图8c所示，可以发现不同时刻、不同中心点射线参数对应大小不同的相似系数，而不同中心点射线参数又分别对应着不同出射方向的一组炮检点高斯束，对于其中某一时刻，与最大相似系数对应的中心点射线参数控制着高斯束在炮检点主要的出射方向，即优势出射方向。而一些较小相似系数可能对应一些随机噪音的出射方向。因此，如果想要在偏移结果中压制由随机噪音引起的干扰，需要计算得到高斯束的优势出射方向，即找到最大相似系数对应的τ-p道集。为了解决上面的问题，我们通过设定相似系数阈值对τ-p道集进行滤波，图8e和图8f是分别采用阈值参数tp=0.4和0.6进行滤波得到的结果，对比可以发现，当取合适的阈值参数时随机噪声在τ-p道集中的能量可以得到一定压制（蓝色椭圆）；但是，过大的阈值参数也可能压制了部分有效能量（红色椭圆），使得有效能量得不到成像。
![]()
图 7 交错网格各向异性高阶有限差分正演模拟的共偏移距（200 m）道集a. 无噪共偏移距道集，b. 加噪共偏移距道集。Figure 7. Common-offset gathers （200 m） obtained by staggered grids anisotropic finite difference forwarda. no-noise common-offset gathers，b. common-offset gathers with noise.![]()
图 8 图7中局部共偏移距道集相干性分析a. 图7b蓝框中局部共偏移距道集，b. τ-p道集，c. 相似系数，d. 图7a蓝框中局部共偏移距道集生成的标准τ-p道集，e. 相似系数阈值tp=0.4滤波后的τ-p道集，f. 相似系数阈值tp=0.6滤波后的τ-p道集。Figure 8. Semblance analysis of local common-offset gathers in Fig.7a. local common-offset gathers marked by blue box in Fig. 7b，b. τ-p gather，c. semblance coefficient，d. standard τ-p gather，e. filtering by tp=0.4，f. filtering by tp=0.6.图9为偏移距为200 m共偏移距域各向异性高斯束及控制束偏移结果，整体上看，常规高斯束偏移结果（图9b），可以发现存在明显的偏移噪声，控制束偏移结果（图9c和9d）中噪声得到明显压制，成像剖面的信噪比显著提高，这对于复杂地区构造解释是非常有用和不可或缺的。同时可以注意到，tp=0.6偏移剖面（图9d）与tp=0.4偏移剖面（图9c）相比而言，阈值较大时信噪比更高，特别是浅层同相轴更加清晰（蓝色矩形框），更接近于图9a原始道集偏移结果，但同时也会损失部分有效地层能量，如模型深部蓝色箭头所示。
![]()
2.3 海上数据处理结果
最后测试新方法对海洋共偏移距实际资料的适应性。图10a、10b和10c分别为中国东部某海上探区的二维各向异性模型的VP0、ε和δ参数场，该地区包含一个大断层及洼陷构造。此各向异性模型参数场是通过传统射线层析方法构建的，海上接收的共偏移距道集已经过带通滤波等预处理，偏移距范围从150 m到4 650 m，偏移距间隔为100 m，共46个共偏移距道集，记录总时间为3.6 s，采样间隔为2 ms，其中偏移距为150 m的共偏移距道集如图10d所示。由图10d可知，该共偏移距道集中含有较为明显的背景及随机噪声，资料信噪比较低，适合进行控制束偏移处理。
![]()
图 10 中国东部某海上探区实际数据a. VP0，b. ε，c. δ，d. 150 m偏移距道集。Figure 10. The field data of marine survey in Eastern Chinaa. VP0，b. ε，c. δ，d. common-offset gather （150 m）.基于相同的成像参数，分别采用各向同性和各向异性共偏移距域高斯束偏移对该实际资料进行处理，结果如图11a和11b所示。为了更加清楚地对比，将图11中白色矩形框内成像结果进行放大显示（图11c和图11d）。从放大图中可以看出，由于考虑地下各向异性参数（更接近于真实地下介质情况）对偏移结果的影响，可以大幅改善高角度构造的成像质量，其结果更能反映地下实际地质构造，陡倾断面更加清晰、连续。为了测试不同相似系数阈值参数对偏移结果的影响，分别采用阈值参数tp=0、0.4、0.6和0.8进行测试得到如图12a-d所示的偏移结果。同样通过对其蓝色矩形框内的偏移结果局部放大对比（图12e-h）可以看出，当阈值参数较小时，成像噪音较严重，当取合适的阈值参数时可大幅提高成像结果的信噪比。但是，过大的阈值参数也可能压制了部分有效能量，使得有效能量得不到成像。所以针对不同的数据资料，需要通过局部共偏移距道集相干性分析来选取合适的阈值参数，使噪音得到压制，有效能量得到保持，最终得到一个较为理想的高精度偏移剖面。
![]()
图 11 海上数据共偏移距域高斯束偏移结果及局部放大对比a. 各向同性，b. 各向异性，c. 各向同性局部放大，d. 各向异性局部放大。Figure 11. Results of GBM and magnified views of the marine surveya. isotropic，b. anisotropic，c. isotropic magnified view，d. anisotropic magnified view.![]()
图 12 海上数据共偏移距域控制束偏移结果及局部放大对比a. tp=0，b. tp=0.4，c. tp=0.6，d. tp=0.8，e. tp=0局部放大，f. tp=0.4局部放大，g. tp=0.6局部放大，h. tp=0.8局部放大。Figure 12. Results of CBM and magnified views of the marine surveya. tp=0，b. tp=0.4，c. tp=0.6，d. tp=0.8，e. tp=0 magnified view，f. tp=0.4 magnified view，g. tp=0.6 magnified view，h. tp=0.8 magnified view.3. 结论
本文根据有效信号和干扰信号在τ-p域中的相干性差异，基于各向异性运动学和动力学射线追踪，发展了一种VTI介质共偏移距域数据驱动控制束偏移理论方法。通过各向异性洼陷模型、修改的SEG/Hess VTI模型及海上实际资料成像处理，得到如下几点认识：
（1）各向异性参数对共偏移距道集中的大偏移距信息成像影响较大。
（2）当地层各向异性不能忽略时，新方法能够更加准确地恢复地下的复杂构造。
（3）新方法能够在一定程度上提高低信噪比数据的偏移成像效果。同样也应该注意到，在实际资料处理中，相似系数阈值参数的选取也十分关键。阈值较小时偏移噪音较强，但过大阈值也可能压制了部分有效能量或者产生偏移假象，选取合适的阈值参数才能够得到较为理想的控制束偏移剖面。因此，下一步拟针对相似系数阈值参数的选取开展深入研究。
-
![]()
图 1 南黄海盆地及邻区构造区划简图
Figure 1. Tectonic map of the South Yellow Sea Basin and its adjacent areas
![]()
![]()
图 3 南黄海盆地东北凹构造演化
Figure 3. The tectonic evolution section of north-east sag in the South Yellow Sea Basin
![]()
![]()
表 2 南黄海盆地东北凹伸缩率计算结果
Table 2 The extensional and compressional rate of north-east sag，the South Yellow Sea Basin
时代 T3 J1 J2+3 K E1+2 E3 剖面A-A' 2.6% 3.5% −3.2% 3.2% 5.6% −4.5% 剖面B-B' 3.7% 4.5% −6.2% 1.5% 4.4% −5.2% 剖面C-C' 1.8% 2.6% −7.6% 2.3% 7.0% −3.0% 
下载: 导出CSV
-
[1] 高顺莉, 周祖翼. 南黄海盆地东北凹侏罗纪地层的发现及其分布特征[J]. 高校地质学报, 2014, 20(2):286-293. [GAO Shunli, ZHOU Zuyi. Discovery of the Jurassic strata in the north-east sag of South Yellow Sea [J]. Geological Journal of China Universities, 2014, 20(2): 286-293. doi: 10.3969/j.issn.1006-7493.2014.02.013 [2] 高顺莉, 徐曦, 周祖翼. 南黄海北部盆地晚白垩世以来构造变形与盆地成因[J]. 石油与天然气地质, 2015, 36(6):924-933. [GAO Shunli, XU Xi, ZHOU Zuyi. Structural deformation and genesis of the northern sub-basin in South Yellow Sea since Late Cretaceous [J]. Oil & Gas Geology, 2015, 36(6): 924-933. doi: 10.11743/ogg20150607 [3] 张海军, 曹强, 孙建峰. 地层剥蚀厚度恢复研究: 以南黄海北部盆地东北凹陷为例[J]. 中国石油勘探, 2008, 13(5):21-25. [ZHANG Haijun, CAO Qiang, SUN Jianfeng. Study on strata eroded thickness restoration: a case study of northeast depression in northern South Yellow sea basin [J]. China Petroleum Exploration, 2008, 13(5): 21-25. doi: 10.3969/j.issn.1672-7703.2008.05.004 [4] 郭飞飞, 曹强, 唐文旭. 南黄海北部盆地东北凹陷地层剥蚀厚度恢复[J]. 海洋地质与第四纪地质, 2009, 29(5):99-105. [GUO Feifei, CAO Qiang, TANG Wenxu. Study on strata erosion thickness restoration in Northeast Depression of South Yellow Sea basin [J]. Marine Geology & Quaternary Geology, 2009, 29(5): 99-105. [5] 李楠, 李巍然, 龙海燕. 南黄海盆地北部坳陷正反转构造[J]. 海洋地质与第四纪地质, 2013, 33(3):95-100. [LI Nan, LI Weiran, LONG Haiyan. Study of positive inversion structures in the North depression of the South Yellow Sea basin [J]. Marine Geology & Quaternary Geology, 2013, 33(3): 95-100. [6] 雷宝华, 陈建文, 梁杰, 等. 印支运动以来南黄海盆地的构造变形与演化[J]. 海洋地质与第四纪地质, 2018, 38(3):45-54. [LEI Baohua, CHEN Jianwen, LIANG Jie, et al. Tectonic deformation and evolution of the South Yellow Sea basin since Indosinian movement [J]. Marine Geology & Quaternary Geology, 2018, 38(3): 45-54. [7] 韩宗珠, 刘涵, 许红, 等. 南黄海盆地盆山耦合作用与地质演化[J]. 海洋地质前沿, 2015, 31(4):55-62. [HAN Zongzhu, LIU Han, XU Hong, et al. Geological evolution and basin-range coupling relationship of the south Yellow sea basin [J]. Marine Geology Frontiers, 2015, 31(4): 55-62. [8] Pang Y M, Zhang X H, Guo X W, et al. Basin modeling in the initial stage of exploration: A case study from the North Subbasin of the South Yellow Sea Basin [J]. Acta Oceanologica Sinica, 2017, 36(9): 65-78. doi: 10.1007/s13131-017-1112-1
[9] Zhang M H, Xu D S, Chen J W. Geological structure of the Yellow Sea area from regional gravity and magnetic interpretation [J]. Applied Geophysics, 2007, 4(2): 75-83. doi: 10.1007/s11770-007-0011-1
[10] 王巍, 陈高, 王家林, 等. 苏北南黄海盆地区域构造特征分析[J]. 地震学刊, 1999(1):47-55. [WANG Wei, CHEN Gao, WANG Jialin, et al. Analysis for regional structural characteristics of north Jiangsu-South Yellow Sea basin [J]. Journal of Seismology, 1999(1): 47-55. [11] 郝天珧, SUH M, 王谦身, 等. 根据重力数据研究黄海周边断裂带在海区的延伸[J]. 地球物理学报, 2002, 45(3):385-397. [HAO Tianyao, SUH M, WANG Qianshen, et al. A study on the extension of fault zones in Yellow Sea and its adjacent areas based on gravity data [J]. Chinese Journal of Geophysics, 2002, 45(3): 385-397. doi: 10.3321/j.issn:0001-5733.2002.03.010 [12] 胥颐, 郝天珧, 李志伟, 等. 中国边缘海域及其邻区的岩石层结构与构造分析[J]. 地学前缘, 2008, 15(3):55-63. [XU Yi, HAO Tianyao, LI Zhiwei, et al. Analysis of lithospheric structure and tectonics of the Chinese marginal seas and adjacent region [J]. Earth Science Frontiers, 2008, 15(3): 55-63. doi: 10.3321/j.issn:1005-2321.2008.03.003 [13] 李慧君. 南黄海盆地海相中—古生界构造区划[J]. 海洋地质前沿, 2014, 30(10):8-13. [LI Huijun. The tectonic division of the Mesozoic-Paleozoic in the south Yellow sea basin [J]. Marine Geology Frontiers, 2014, 30(10): 8-13. [14] 周建勋. 同沉积挤压盆地构造演化恢复的平衡剖面方法及其应用[J]. 地球学报, 2005, 26(2):151-156. [ZHOU Jianxun. The balanced cross-section method for restoration of structural evolution in compressional basins with synkinematic sedimentation and its application [J]. Acta Geoscientica Sinica, 2005, 26(2): 151-156. doi: 10.3321/j.issn:1006-3021.2005.02.008 [15] 李三忠, 索艳慧, 戴黎明, 等. 渤海湾盆地形成与华北克拉通破坏[J]. 地学前缘, 2010, 17(4):64-89. [LI Sanzhong, SUO Yanhui, DAI Liming, et al. Development of the Bohai Bay Basin and destruction of the North China craton [J]. Earth Science Frontiers, 2010, 17(4): 64-89. [16] 万天丰, 朱鸿. 中国大陆及邻区中生代—新生代大地构造与环境变迁[J]. 现代地质, 2002, 16(2):107-120. [WAN Tianfeng, ZHU Hong. Tectonics and environment change of Meso-Cenozoic in China continent and its adjacent areas [J]. Geoscience, 2002, 16(2): 107-120. doi: 10.3969/j.issn.1000-8527.2002.02.001 [17] Li S Z, Jahna B M, Zhao S J, et al. Triassic southeastward subduction of North China Block to South China Block: Insights from new geological, geophysical and geochemical data [J]. Earth-Science Reviews, 2017, 166: 270-285. doi: 10.1016/j.earscirev.2017.01.009
[18] 朱光, 刘国生, 李双应, 等. 下扬子地区盆地的“四层楼”结构及其动力学机制[J]. 合肥工业大学学报: 自然科学版, 2002, 23(1):47-52. [ZHU Guang, LIU Guosheng, LI Shuangying, et al. “Four-story” structure of the basins in the Lower Yangtze region and their geodynamic setting [J]. Journal of Hefei University of Technology: Natural Science, 2002, 23(1): 47-52. [19] 朱光, 刘程, 顾承串, 等. 郯庐断裂带晚中生代演化对西太平洋俯冲历史的指示[J]. 中国科学: 地球科学, 2018, 61(4):386-405. [ZHU Guang, LIU Cheng, GU Chengchuan, et al. Oceanic plate subduction history in the western Pacific Ocean: Constraint from late Mesozoic evolution of the Tan-Lu Fault Zone [J]. Science China Earth Sciences, 2018, 61(4): 386-405. doi: 10.1007/s11430-017-9136-4 [20] Maruyama S, Isozaki Y, Kimura G, et al. Paleo-geographic maps of the Japanese islands: plate tectonic synthesis from 750 Ma to the present [J]. The Island Arc, 1997, 6(1): 121-142. doi: 10.1111/j.1440-1738.1997.tb00043.x
-
期刊类型引用(3)
1. 王首良，李元昊，马婷钰，段祎乐. 滨里海盆地早二叠世孔谷期盐构造特征及其形成机制. 海洋地质前沿. 2024(01): 65-78 . 
百度学术
2. 韩续，索艳慧，李三忠，丁雪松，宋双双，田子晗，付新建. 新近纪以来华北东部古地貌演化数值模拟及陆架海沉降控制. 古地理学报. 2024(01): 192-207 . 百度学术
3. 陈念楠，李满根，关宝文，宋志杰，段建兵，李西得，刘武生，刘颖，范鹏飞. 二连盆地塔北凹陷西部早白垩世断—坳发育特征研究. 物探与化探. 2022(05): 1149-1156 . 百度学术
其他类型引用(0)
计量
- 文章访问数: 1922
- HTML全文浏览量: 504
- PDF下载量: 25
- 被引次数: 3
