何起祥先生生平及主要学术成就

天然气水合物是一种重要的战略能源和潜在未来能源，准确地探明中国的天然气水合物资源储量，有效地进行水合物资源的勘探开发，具有非常重要的战略意义、环境意义和经济意义。地震勘探是天然气水合物识别和预测的重要手段之一。地震资料上似海底反射层（BSR）、振幅空白带以及极性反转等地震反射特征^[1-2]，通常作为海域水合物储层的重要识别标志。地震反射主要是由储层波阻抗差异引起的，但在特定条件下岩性变化或者不整合面等原因也可能造成波阻抗差异，引起类似的地震响应。地震剖面上BSR常见于海域斜坡或者浅部海底沉积地层上，但在一些发现水合物的区域，并未表现出明显的BSR特征^[3]，如墨西哥湾、中国南海北部陆坡西部区域等，而且不同类型水合物储层的BSR特征也有较大差异^[4-5]。BSR等地震反射特征与水合物地层并非一一对应关系，仅依靠地震反射特征难以准确识别天然气水合物^[6]，这对海域水合物储层的地震识别和预测提出了更高的要求。

天然气水合物的微观赋存形态是影响其储层特性的主要因素^[7-8]。为提高水合物识别和预测的可靠性，一个重要的途径是从地震数据中准确提取天然气水合物地层的物性信息，为此需要解决两个关键问题：第一个问题是如何建立储层宏观弹性特征（包括体积模量、剪切模量、纵波速度、横波速度、密度等）与天然气水合物微观结构之间的定量关系。岩石物理建模是建立该定量关系的重要手段，目前应用较多的岩石物理模型包括适用于固结海洋沉积物的胶结模型^[9-10]和等效介质模型^[11-12]以及适用于未固结海洋沉积物的接触理论^[13-15]等。这些岩石物理模型各自具有不同的假设条件和适用条件。水合物的微观赋存形态会因沉积环境的不同而显著变化，实验研究表明，水合物的存在会增强地层的剪切强度^[16]。考虑水合物微观赋存形态及剪切性质的岩石物理型有待进一步研究。

第二个问题是如何建立水合物储层地震信息（包括振幅、相位、频率等动力学特征以及BSR、振幅空白带等地球物理响应）与弹性参数之间定量关系。地震反演方法是天然气水合物探测的有效技术手段之一。弹性阻抗反演^[17-23]、AVO（Amplitude variation with offset，振幅随偏移距的变化）及AVA（Amplitude variation with angle，振幅随入射角度的变化）反演方法^[24-26]以及全波形反演^[27-28]等地震反演方法已在天然气水合物探测中取得了成功应用。然而，天然气水合物地层的宏观和微观结构高度复杂、非均质性强，成藏模式随不同区域而变化，基于地震数据进行水合物地层物性信息的定量反演仍存在一定的多解性^[29-30]，这对水合物的定量解释和勘探提出了挑战。目前流体识别的叠前地震反演方法日趋完善，并在多个地区取得了较好的应用效果^[31]，而水合物地层的微观赋存形态以及宏观沉积环境区别于常规含流体储层，因此此类反演方法在水合物地层的应用还有待进一步研究。

针对第一个问题，本文根据海域天然气水合物的微观结构特征，在岩石物理建模的框架下，考虑水合物非零剪切模量的特征，使用基于非耦合微分有效介质(DEM)理论和斑块饱和理论进行岩石物理建模，计算水合物地层物性参数与弹性参数的定量关系；针对第二个问题，在岩石物理分析的基础上，使用叠前贝叶斯地震反演方法，从叠前地震数据中提取水合物地层的剪切模量参数。实际应用结果表明，剪切模量对水合物含量有较好的指示性。本文的方法对水合物储层的识别提供了方法和技术支撑，对后续水合物物性参数定量反演提供新的思路。

1. 工区概况

神狐海域位于南海北部陆坡神狐暗沙东南海域附近（图1），水合物钻探区水深为1 000～1 500 m，新生代沉积厚度达1 000～7 000 m，具有扩散型天然气水合物的成藏系统。适宜的温压条件、丰富的气源供给，加之底辟构造、高角度断裂和垂向裂隙系统等有利的运移通道，为水合物发育创造了良好的条件^[32]。2006年，广州海洋地质调查局在此区域实施了三维地震采集，2007年完成了中国首次天然气水合物钻探，成功获取了天然气水合物的实物样品。

图 1 研究区位置图

a：神狐海域研究区构造(红色方块) ，b：研究区内水合物钻探站位分布

Figure 1. The location of the study area

a: The tectonic map of the Shenhu area (red square), b: the distribution of hydrate drilling stations in the red square (enlarged)

下载: 全尺寸图片幻灯片

神狐海域水合物稳定带受到周围温压条件限制，含水合物沉积层基本上位于海底以下135～230 m，分布于强振幅BSR的上覆沉积地层中^[33]，大部分BSR与地层斜交，波形极性与海底同相轴极性相反，且上部发育弱振幅或者振幅空白带^[34]（图2）。横向上BSR反射振幅及垂向上BSR厚度均存在一定变化，主要是由于地层的孔隙度、渗透率等物性差异导致水合物饱和度不同而造成的。

图 2 神狐海域典型地震反射剖面

Figure 2. A typical seismic reflection profile in the Shenhu area

下载: 全尺寸图片幻灯片

2. 弹性参数计算

针对前述的第一个问题，本文采用岩石物理建模的方法，建立地层宏观弹性参数与天然气水合物微观结构之间的定量关系。陆红锋等^[35]对岩芯进行了X射线衍射成分分析，指出神狐海域沉积物组分以碎屑矿物、黏土矿物和碳酸盐矿物为主，并给出了各组分的含量，据此假设岩石基质由方解石、石英、黏土等组成。水合物含量主要受地层岩性、孔隙的微观结构及表面性质的影响^[8]，神狐海域水合物主要分散状分布于矿物颗粒的孔隙中^[36-37]，由于南海沉积物中存在大量有孔虫壳体，水合物主要占据有孔虫壳体内部空间并堵塞了孔隙间的连通喉道，显著降低了沉积物的气、水渗透能力^[38]；研究还发现，黏土含量的增加会导致被黏土颗粒束缚的孔隙水含量也随之增加^[39]。据此考虑孔隙几何形状和连通性，假设地层中孔隙空间包含不连通的束缚水孔隙以及连通的硬币状孔隙和椭球形孔隙，其中硬币状孔隙纵横比较小，受到压力易趋向于闭合，连通性较差，而椭球形孔隙纵横比较大，连通性较好。孔隙空间中可能存在水合物、游离气和水共存的情况^[40]，该地区的水合物呈现出典型的I型结构，甲烷气体的含量超过99.3%^[41]，据此，结合前人对水合物微观结构的研究结果（图3a、b），将水合物地层等效为图3c所示的等效介质。

图 3 实验室获取的水合物CT图像^[42]及水合物等效介质模型

a：原始图像，b：对a图红色框线内区域进行灰度分析、着色后放大图，c：等效介质模型。

Figure 3. CT images of gas hydrate obtained in laboratory and equivalent medium model of gas hydrate-bearing sediments

a: Original image, b: enlarged image of the area in the red frame after grey analysis and coloring, c: equivalent medium model of gas hydrate-bearing sediments.

下载: 全尺寸图片幻灯片

在地震岩石物理建模的理论框架下建立天然气水合物地层的岩石物理模型（主要流程见图4）：

图 4 水合物地层岩石物理建模流程示意图

Figure 4. Flowchart of the rock physics modelling

下载: 全尺寸图片幻灯片

①利用Hashin-Shtrikman-Walpole界限^[43]将矿物组分混合形成颗粒基质。

②在基质中加入充填束缚水的不连通孔隙，作为新的基质。束缚水孔隙纵横比假设为0.1。根据Kuster-Toksöz模型计算基质中加入束缚水孔隙之后的弹性参数：

$$\begin{split}& \left( {{\boldsymbol{K}} - {{\boldsymbol{K}}_{gr}}} \right){\left( {{\boldsymbol{K}} + 4{{{{\boldsymbol{R}}_{gr}}} \mathord{\left/ {\vphantom {{{{\boldsymbol{R}}_{gr}}} 3}} \right. } 3}} \right)^{ - 1}}\left( {{{\boldsymbol{K}}_{gr}}+ 4{{{{\boldsymbol{R}}_{gr}}} \mathord{\left/ {\vphantom {{{{\boldsymbol{R}}_{gr}}} 3}} \right. } 3}} \right) = \\&\quad{\phi _{wc}}\left( {{{\boldsymbol{K}}_w} - {{\boldsymbol{K}}_{gr}}} \right){{\boldsymbol{P}}^{wc}} \end{split}$$

(1)

其中，${\boldsymbol{K}} = (K,\mu )$，$ K $、$ \mu $分别为基质中加入束缚水孔隙之后的体积模量和剪切模量；${{\boldsymbol{K}}_{{\text{gr}}}} = ({K_{{\text{gr}}}},{\mu _{{\text{gr}}}})$，${{\boldsymbol{K}}_{\text{w}}} = ({K_{\text{w}}},{\mu _{\text{w}}})$分别为岩石基质和水的弹性模量；$ {\phi _{{\text{wc}}}} $为束缚水孔隙的孔隙度，${\boldsymbol{R}} = \left( {\mu ,\dfrac{\mu }{8}\dfrac{{9K + 8\mu }}{{K + 2\mu }}} \right)$；系数${{\boldsymbol{P}}^{{\text{wc}}}} = ({P^{{\text{wc}}}},{Q^{{\text{wc}}}})$是几何因子，由孔隙纵横比确定。

③利用非耦合DEM模型^[44]向②得到的基质中加入干燥的硬币状孔隙和椭球形孔隙，形成干燥骨架：

$$ K_{\text{d}}^{}\left( \phi \right){\text{ = }}{K_1}{(1 - \phi )^{\left( {{P^0} + {P^1}} \right)}}{e^{{P^1}\phi }} $$

(2)

$$ \mu _{\text{d}}^{}\left( \phi \right){\text{ = }}{\mu _1}{(1 - \phi )^{\left( {{Q^0} + {Q^1}} \right)}}{e^{{Q^1}\phi }} $$

(3)

其中，$K_{\text{d}}^{}$和$\mu _{\text{d}}^{}$分别是干燥骨架的体积模量和剪切模量；${P^0} = \dfrac{{{A_1}}}{{{A_4}}}$，${P^1} = \dfrac{{{A_2}{A_4} - {A_1}{A_5}}}{{A_4^2}}$，${Q^0} = \dfrac{{{C_1}}}{{{C_4}}} + \dfrac{{{D_1}}}{{{D_4}}}$，${Q^1} = \dfrac{{{C_2}{C_4} - {C_1}{C_5}}}{{C_4^2}} + \dfrac{{{D_2}{D_4} - {D_1}{D_5}}}{{D_4^2}}$，$ A_1,A_2.\ldots A_5 $，$ C_1,C_2.\ldots C_5 $和$ D_1,D_2.\ldots D_5 $等是与固体基质和孔隙纵横比有关的参数，具体表达式见文献[44]。根据Xu和White（1995）的岩石物理建模思想^[45]，假设硬币状孔隙的孔隙度与黏土体积含量成正比，其余的孔隙为椭球形孔隙。

④利用广义Gassmann方程^[46]将水合物填充到硬币状孔隙中，形成水合物斑块：

$$ K_{\text{s}}^{ - 1} = K_{\text{d}}^{ - 1} - {\left( {K_{\text{d}}^{ - 1} - K_{\text{m}}^{ - 1}} \right)^2}{\left[ {\phi \left( {K_{if}^{ - 1} - K_\phi ^{ - 1}} \right) + \left( {K_{\text{d}}^{ - 1} - K_{\text{m}}^{ - 1}} \right)} \right]^{ - 1}} $$

(4)

$$ \mu _{\text{s}}^{ - 1} = \mu _{\text{d}}^{ - 1} - {\left( {\mu _{\text{d}}^{ - 1} - \mu _{\text{m}}^{ - 1}} \right)^2}{\left[ {\phi \left( {\mu _{if}^{ - 1} - \mu _\phi ^{ - 1}} \right) + \left( {\mu _{\text{d}}^{ - 1} - \mu _{\text{m}}^{ - 1}} \right)} \right]^{ - 1}} $$

(5)

其中，$ {\mu _{\text{s}}} $是水合物充填斑块的剪切模量，$ {K_\phi } $和$ {\mu _\phi } $分别是与孔隙相关的体积模量和剪切模量，$ {K_{if}} $和$ {\mu _{if}} $分别是第$i$种孔隙充填物的体积和剪切模量。利用广义Gassmann方程可以得到固体（如有机质、水合物等）充填于孔隙情况下的饱和岩石的弹性模量。

利用Gassmann方程将水/气填充到椭球形孔隙中，形成水/气斑块；最后，利用斑块饱和模型^[47]和Reuss平均^[48]将水合物斑块和水/气斑块混合，形成饱和岩石：

$$ {K_{\text{R}}} = {\left\langle {{{\left( {K + 4{\mu \mathord{\left/ {\vphantom {\mu 3}} \right. } 3}} \right)}^{ - 1}}} \right\rangle ^{ - 1}} - 4{\mu \mathord{\left/ {\vphantom {\mu 3}} \right. } 3} $$

(6)

$$ {\mu _{\text{R}}} = {\left( {\sum\limits_{i = 1}^N {{{{f_i}} \mathord{\left/ {\vphantom {{{f_i}} {{\mu _i}}}} \right. } {{\mu _i}}}} } \right)^{ - 1}} $$

(7)

其中，$ {K_{\text{R}}} $和${\mu _{\text{R}}}$分别表示饱和沉积物的等效体积和剪切模量；${f_i}$为其体积含量；$ \left\langle \cdot \right\rangle $表示对每个斑块按体积含量进行加权平均；$ {\mu _i} $为第$ i $个斑块的剪切模量。

根据上述岩石物理模型计算水合物地层的弹性参数。假设地层的矿物组分为石英和黏土，孔隙度为35%，硬币状和椭球形孔隙的纵横比分别假设为0.04和0.35，孔隙中充填水合物和水，选取的其余弹性常数如表1所示。剪切模量和纵波速度随水合物剪切模量的变化如图5所示，两者都随水合物剪切模量的增大而增大，水合物饱和度越高，地层剪切模量随水合物剪切模量的变化越大。

表 1 数值计算所选用的弹性常数

Table 1. The elastic constants used for numerical calculation

组分	体积模量/GPa	剪切模量/GPa	密度/(g/cm³)	来源
石英	37.0	44.0	2.65	Carmichale, 1989
黏土	21.0	7.0	2.60	Tosaya and Nur, 1982
水合物	7.7	3.2	0.91	Waite et al, 2000
水	2.25		1000

下载: 导出CSV

| 显示表格

图 5 水合物地层纵波速度(a)和剪切模量(b)随水合物剪切模量的变化

Figure 5. P-wave velocity (a) and shear modulus (b) of gas hydrate-bearing sediments versus shear modulus of gas hydrate component

下载: 全尺寸图片幻灯片

使用下式计算弹性参数对水合物含量的敏感性评价指标：

$$ {D_{{S_i}}} = \left| {\frac{{{P_{{S_i}}} - {P_{{S_0}}}}}{{{P_{{S_0}}}}}} \right| $$

(8)

其中，${D_{{S_i}}}$和${P_{{s_i}}}$分别表示水合物饱和度为${S_i}$时的敏感性评价指标和弹性参数；${P_{{s_0}}}$表示地层不含水合物时的弹性参数。计算不同弹性参数对水合物含量的敏感性评价指标（图6）发现，纵波速度、横波速度、泊松比、体积模量和剪切模量等弹性参数对水合物含量都有一定的敏感性，但是泊松比的敏感性较低，纵波速度和横波速度对水合物含量的敏感性相当，剪切模量对水合物含量的敏感性最高。这与数值计算的结论是一致的，也证明了使用剪切模量参数来识别水合物的可行性。

图 6 不同弹性参数的敏感性评价指标

Figure 6. Sensitivity indicator for hydrate content of different elastic parameters

下载: 全尺寸图片幻灯片

2007年在神狐海域共实施了8个天然气水合物钻探站位，其中SH2、SH3和SH7站位钻遇了水合物样品。其中SH2站位水合物厚度约25 m，位于海底以下195～220 m，平均饱和度25%^[49]，最高值达48%，部分测井曲线如图7所示。

图 7 SH2站位测井曲线

绿色层段表示水合物地层。

Figure 7. Well logging data of SH2 site

The green section represents hydrate-bearing sediments.

下载: 全尺寸图片幻灯片

针对本研究区内的SH2站位，使用上述建模方法计算站位处的弹性参数。硬币状孔隙和椭球形孔隙各自的纵横比根据实测纵波速度反演求得。计算得到的纵横波速度以及剪切模量如图8所示。计算的纵波速度和与实测纵波速度有较好的一致性，两者相关系数为0.824，方差为0.0182，绝对误差平均值为48.49 m/s。对计算得到的弹性参数进行交会分析（图9），实测纵波速度和自然伽马对水合物含量的区分能力有限，而剪切模量对于水合物含量有较好的指示作用，较高水合物饱和度的地层对应着较高的剪切模量。

图 8 SH2站位水合物地层的水合物饱和度及弹性参数计算结果

a：根据实测电阻率曲线计算的水合物饱和度，b：纵波速度(红色曲线表示实测值，黑色曲线表示计算值)，c：横波速度计算值；d：剪切模量计算值。

Figure 8. Calculation results of hydrate saturation and elastic parameters of hydrate strata at SH2 site

a: hydrate saturation calculated according to the measured resistivity, b: P-wave velocity (red curve represents the measured value, black curve represents the calculated value), c: calculated shear wave velocity; d: calculated shear modulus.

下载: 全尺寸图片幻灯片

图 9 SH2站位水合物地层弹性参数交会图

a：实测纵波速度和实测自然伽马交会图，b：计算的体积模量和剪切模量交会图。

Figure 9. Crossplots of elastic parameters of hydrate-bearing sediments at SH2 site

a: Crossplot of measured P-wave velocity and measured Gamma, b: crossplot of calculated bulk modulus and shear modulus.

下载: 全尺寸图片幻灯片

3. 叠前地震反演

关于纵波模量和剪切模量的弹性阻抗方程，其标准化形式为^[50]：

$$ {\text{EI}}(\theta ) = {({M_0}{\rho _0})^{\frac{1}{2}}}{\left( {\frac{M}{{{M_0}}}} \right)^{a(\theta )}}{\left( {\frac{\mu }{{{\mu _0}}}} \right)^{b(\theta )}}{\left( {\frac{\rho }{{{\rho _0}}}} \right)^{c(\theta )}} $$

(9)

其中，$ {\text{EI}} $为弹性阻抗，$ M $和$ \mu $分别为饱和岩石的纵波模量和剪切模量（横波模量），$a = \frac{1}{2}{\sec ^2}\theta $，$ b = - 4\gamma _{}^2{\sin ^2}\theta $，$c = 1 - \frac{1}{2}{\sec ^2}\theta $，$\theta $为入射角。

结合岩石物理框架下弹性参数的计算结果，从测井资料以及地震资料出发，利用弹性阻抗贝叶斯反演方法计算得到弹性阻抗，然后使用上式进行纵横波模量的直接提取。该方法不需要先提取纵横波速度及密度信息，因此避免了从弹性阻抗计算纵横波速度再计算剪切模量过程中的累积误差。

弹性阻抗贝叶斯反演得到的弹性阻抗体以及直接反演得到的水合物地层剪切模量如图10所示。反演结果整体横向连续性较好，水合物地层分属高剪切模量数值区间，与岩石物理交会分析的认识相吻合；测井纵波速度增大的位置与剪切模量较高的位置匹配较好，反演结果与SH2井的实际测井钻遇结果保持较高一致性，验证了本工区使用剪切模量直接反演方法进行水合物识别的可靠性。

图 10 SH2站位水合物地层弹性参数交会图

蓝色曲线为SH2井的实测纵波速度曲线。a：中角度叠加地震剖面，中心角度为20°；b：中角度弹性阻抗剖面，中心角度为20°；c：反演得到的剪切模量剖面。

Figure 10. Seismic, elastic impedance, and inverted shear modulus profiles of SH2 site

The blue curve is the measured P-wave velocity of SH2 site. a: Partial stack seismic profile at centric angle of 20°; b: elastic impedance profile at centric angle of 20°; c: inverted shear modulus profile.

下载: 全尺寸图片幻灯片

4. 结论

（1）考虑神狐海域天然气水合物的微观结构特征，基于测井数据和地震数据，使用岩石物理建模以及叠前贝叶斯地震反演方法，从叠前地震数据中提取水合物地层的剪切模量参数，与SH2站位测井钻遇结果保持较高一致性，验证了使用剪切模量反演方法进行水合物识别的可靠性和有效性。

（2）在岩石物理建模过程中，对于水合物地层的微观结构，本文主要是通过矿物组分、孔隙结构以及水合物和矿物颗粒的接触关系来描述。实际水合物地层的微观结构随沉积环境以及富集条件的变化而差别非常大，因此等效岩石物理模型所用的一些物性参数，比如孔隙形状以及各类孔隙的孔隙度，都是对实际地层微观结构的等效近似。针对水合物地层的矿物组分以及水合物在地层中的赋存状态，许多学者通过CT图像以及实验等途径，给出了具体的描述和大量的数据，这些都是本文进行岩石物理建模的依据。

（3）本文的方法适用于孔隙充填型水合物地层，为该类水合物地层的识别和预测提供了一种新的思路和可靠的方法。但是海域水合物的赋存状态多样，包括颗粒支撑型、裂缝充填型等，不同类型的水合物需要通过岩石物理建模和交会分析等方法查明弹性参数对水合物含量的敏感性，确定敏感弹性参数，再使用地震反演方法进行敏感弹性参数的提取。

下载: 全尺寸图片幻灯片