技术领域
[0001] 本
发明涉及勘探地球物理领域，具体涉及一种基于地震海洋学的海洋
水体物性参数测量方法，是将地震反射方法和海洋水体参数研究的重要方法流程。
背景技术
[0002] 本发明所涉及的方法是基于地震海洋学，结合反演方法预测海洋水体物性参数的方法。2003年前，各国科学家对
海水地震反射资料的处理中，发现弱反射结构。但直到Holbrook重新处理北大西洋中纽芬兰大滩附近的海水地震资料，利用海洋反射
地震勘探方法得到海水水体结构的
叠加剖面，才正式提出地震海洋方法。Holbrook等人发表在Science期刊的文章，引起众多国家的海洋地球物理学家和物理海洋学家的广泛关注。该篇文章分析发现北大西洋洋流(NAC)和拉布拉多寒流(LC)海洋峰面处的海水混合与交换过程的精细结构中，有至少4种反射层类型，这说明利用地震海洋方法可以发现并分析物理海洋学的相关海水细结构。国内地震海洋学的研究一直关注地震剖面的结构性，通过结构以及反射层的对比进一步研究海洋水体的特征。本方法从地震反射原始数据出发，结合地震反射精细处理方法，联合拟
测井CTD数据与地震合成记录，通过反演方法以及时频分析方法得到相应的三瞬、速度、
密度、
温度以及
盐度等海洋水体物性参数，完善了地震海洋学的方法体系，提高了海洋水体的地震海洋学研究手段精细程度，并且该方法适用于任何深度的海洋水体研究，不再局限于深水海洋水体范围。
发明内容
[0003] 本次发明目的是提供一种基于地震海洋学的海水物性测量方法，旨在提高完善地震海洋学处理与分析海洋水体的地震反射资料，通过联合拟测井CTD数据和残差计算控制处理叠加剖面的效果，最终反演预测海洋水体物性参数，应用到海洋学的水体研究。
[0004] 地震海洋学是利用海洋水体的弱反射数据进行物理海洋方面研究的交叉性学科方法。海水的温盐结构随深度的变化而呈现不同，从而影响到在不同温盐层上表现出
声波波阻抗的差别，弹性波在不同深度的水层上会表现出
波动反射现象。因此，水体温盐结构的变化在反射地震剖面上清晰体现出来。本方法首先对原始数据进行常规处理，利用真振幅恢复、顶切底切、
频率扫描、滤波，预测反褶积，速度谱扫描以及叠加等处理方法，得到叠加剖面。可以从叠加剖面上看出海洋水体弱反射的非连续同相轴或者海洋水体中的异常反射，从而确定精细处理部分的目标层位。后面可以有针对性的进行处理。通过变换到共检波点域，抓住噪声在该域中的特征，利用曲波变换进行干扰压制，达到原始资料精细处理与去噪效果的最优化。
[0005] 设原始地震数据为D，有效
信号部分为S，同理噪声部分为N，之间关系满足关系：
[0006] D＝S+N
[0007] 共检波点域曲波变换去噪方法步骤：
[0008] (1)根据不同数据解编方式，将原始数据共炮点域Dshot转换为共检波点域Drec；
[0009] (2)在各个数域中，根据数据大小确定尺度数j和方向数l，对数据进行曲波变换，将数据稀疏化，得到对应曲波系数Cshot(j,l,k)与Crec(j,l,k)；
[0010] (3)对稀疏系数进行
阈值筛选，得到新的系数矩阵C'shot(j,l,k)与C'rec(j,l,k)，对其进行曲波反变换，得到重构地震数据D′shot与D′rec；
[0011] (4)通过差值原始数据与地震重构数据，得到各域一次变换后的去除的噪音成分N'shot与N'rec，将两者转换至同一数域后，对比分析两噪音成分，取两者交集得到噪声估计，认为这就是原始数据中的噪音N1。
[0012] (5)根据去噪后数据
质量的反馈，确定方法返回
迭代次数。
[0013] (6)叠加剖面继续去噪，分离噪音为N2。则认为原始噪音与分离噪音的关系为：N＝N1+N2
[0014] 海洋水体资料经过处理分析后，利用基于小波域多输出系统空间子波估计方法，提取
地震道子波，子波估计的精确度直接会提高地震剖面
分辨率、反褶积效果以及波阻抗反演效果等。将提取出的子波与CTD资料合成地震记录。
[0015] 记一道地震记录为x(n),n＝1,2,…J，构造方式混合
相位子波估计步骤如下：
[0016] (1)构造多输出系统，将地震记录进行小波分解，获得高低频的多输出地震道；
[0017] (2)计算伪子道的二阶累积量Rx；
[0018] (3)子空间分解，得到S＝[s1,…,sN+L+q]，G＝[g1,…,gN+L+q]；
[0019] (4)求解构造矩阵Q，估计出地震子波。
[0020] 将合成地震记录与叠加剖面相对应的CDP道匹配做差，计算残差与阈值之间的大小关系。若残差大于阈值，证明两者存在较大差异，应当返回重新调整
数据处理参数，数据进行循环处理；若残差小于阈值，证明两者匹配较好，直接输出叠后剖面。
[0021] 得到叠加剖面后，利用自适应核时频分析方法，对叠加剖面提取瞬时频率，瞬时振幅以及瞬时相位属性。针对各个剖面进行频率、反射强度以及相位变化的分析，完成海洋水体的三瞬属性预测。
[0022] 连续
小波变换具备时域和频域双重良好的局部性和随尺度变化的自动调焦功能，而自动调焦可以通过母小波的伸缩和平移来实现。首先给出母小波函数ψ(t)的定义，下式小波称为母小波函数。
[0023]
[0024] 其中，Ψ(ω)是ψ(t)的Fourier变换。
[0025] 对ψ(t)做平移和伸缩，得到一簇小波函数ψa,b(t)。
[0026]
[0027] 其中，a为尺度参数，b为平移参数。
[0028] 对于给定的信号f(t)∈L2(R)，其连续小波变换的定义为
[0029]
[0030] 其中，Wf(a,b)是信号f(t)的小波系数， 是ψa,b(t)的共轭函数。内积运算反映了信号的局部特征和小波函数的相似性，内积越大，相似性越高。
[0031] 在计算过程中，规定连续小波变换尺度变量的表达式为
[0032]
[0033] 其中，fw是母小波的中心频率如下式所示，scal是尺度的个数，s＝scal,scal-1,...1。
[0034]
[0035] 其中，ξ(f)是母小波的Fourier变换，f1和f2分别为母小波频带的上下限。
[0036] 每个尺度下的小波函数都有一个中心频率，其关系式为：
[0037]
[0038] 其中，fa是尺度a对应的频率，dt是待分析信号的
采样间隔，在上式的约束下，尺度变量对应的频率fa的采样是等间隔的。某个尺度下小波变换的结果即是待分析信号与具有某个确定中心频率的小波函数的相似性，反映出待分析信号在该频率下的
能量分布特征。
[0039] 利用小波系数，可以重构原信号，重构式为：
[0040]
[0041] 叠后约束反演是一种常见的适用于叠后数据反演的方法。稀疏脉冲反演方法，首先假设
地层反射系数序列满足稀疏特性。通过计算出反射系数后，采用最大似然反演方法转换为波阻抗。其最小目标函数为：
[0042]
[0043] 式中，R2表示反射系数的均方值；N2表示噪音的均方值；r(k)表示第K个采样点的反射系数；n(k)表示第K个采样点的噪音；M表示反射层数；L表示采样总数；λ表示给定反射系数的似然值。然后通过多次迭代，求取反射系数。
[0044] 若利用最大似然反褶积获得的反射系数是ri，我们利用最大似然反演的
算法转换反射系数序列得到宽带波阻抗。
[0045] 约束稀疏脉冲反演的关键是使下面这个目标函数得到优化：
[0046] E＝∑(ri)p+λq∑(di-si)q+α2∑(ti-zi)2
[0047] ri为反射系数；λ为控制稀疏脉冲因子，为实际地震记录与合成地震记录残差权重因子，与地震资料的
信噪比的大小有关；di为地震数据；α为数据匹配的权因子；si为合成地震道数据；ti为井阻抗曲线的趋势；zi为实际的阻抗；p,q为L模因子，一般p＝1,q＝2；i地震道采样序号。
[0048] 本方法突出精细处理地震多道数据，利用原始CTD数据作为测井数据，改进提取子波方法，结合常规处理，完成对海洋水体地震资料的精细处理，提高数据真实性与可靠性。整合现有反演方法与地震属性提取方法，获得一套基于采集，处理，属性预测的完整方法，使得海洋水体成像横向分辨率大幅提高。该方法完整，科学，可行以及降低测量成本与提高横向分辨率。是一种优良的海水物性参数测量方法，弥补传统物理海洋测量方法的不足，完善地震海洋学方法体系，并且该方法适用于各种水深的水体资料，具有普适性。
附图说明
[0049] 图1基于地震海洋学的海水物性预测方法
流程图[0050] 图2地震多道精细处理叠加剖面
[0051] 图3叠后剖面的三瞬属性剖面，
[0052] 瞬时频率剖面(图3a)，瞬时振幅剖面(图3b)以及瞬时相位剖面(图3c)
[0053] 图4反演得到的速度、密度、温度、盐度剖面图
[0054] 速度剖面(图4a)，温度剖面(图4b)，盐度剖面(图4c)。
具体实施方式
[0055] 如图1，一种基于地震海洋学的海水物性测量方法，包括以下步骤：
[0056] 步骤100，通过常规方式获得海洋原始资料，包括24道地震原始数据，抛弃式CTD原始资料，以及对应于24道地震原始数据和抛弃式CTD的导航数据；
[0057] 步骤101，对24道地震原始数据进行常规处理，包括真振幅恢复、顶切、速度扫描、滤波，预测反褶积，完成初步地震资料的处理工作；
[0058] 步骤102，进行CRP曲波变换去噪：即基于曲波变换理论
基础，将24道地震原始资料变换为共检波点域(CRP)，对于24道地震原始资料中的环境噪声，利用各类噪声在该共检波点域中的特点，利用曲波变换进行去噪，以突出海水层的弱反射同相轴部分，进而叠加24道地震原始资料，形成地震数据的叠加剖面，达到对地震数据的精细处理效果；
[0059] 其特征在于该方法还包括以下步骤：
[0060] 步骤103，利用上述导航数据，提取出对应于CTD投放点的叠后剖面CDP道；然后利用改进的单输入多输出统计方法进行地震子波估计，提取子波与对应CTD合成地震记录，[0061] 所述的改进的单输入多输出统计方法，即利用基于小波域多输出系统空间子波估计方法，提取地震道子波，获得高精确度的子波估计从而提高地震剖面分辨率、反褶积效果以及波阻抗反演效果，具体如下：
[0062] 首先，将上述与CTD投放点对应的叠后剖面CDP道地震数据记为x(n),n＝1,2,…P，P表示地震数据的样点数，
[0063] (1)将上述地震数据进行小波分解，获得高低频的多输出伪子道X(k)
[0064]
[0065] k为1或2，当k为1时表示高频伪子道，k为2时表示低频伪子道，
[0066] (2)计算上述多输出地震道的伪子道X(k)的二阶累积量Rx
[0067]
[0068] 其中，E表示协方差矩阵，N表示时窗样点数， 为上述多输出伪子道X(k)的转置向量；
[0069] (2)将上述二阶累积量Rx进行特征值分解，并按特征值从大到小排列对应的
特征向量，得到信号空间S和噪音空间G，即S＝[λ1,…,λr]，G＝[λr+1,…,λNP]。其中λi为对应从大到小排序的特征值，r为上述二阶累积量Rx的秩，NP为上述二阶累积量Rx的维数；
[0070] (3)利用上述噪音空间G构造矩阵Q，
[0071]
[0072] 求矩阵Q最小特征值对应的特征向量，即为对应地震伪子波wi，求和
[0073]
[0074] 则估计出地震子波w；
[0075] 将估计出的地震子波w同与其对应的抛弃式CTD原始资料合成地震记录；
[0076] 步骤104，将上一步骤提取出的叠后剖面CDP道与合成地震记录匹配做差，形成残差，然后计算叠后剖面CDP道平均能量，再将叠后剖面CDP道平均能量设定为阈值，比较残差与阈值的大小，若残差小于阈值，则证明叠后剖面CDP道与合成的地震记录匹配好，则将步骤102得到的叠加剖面输出，作为最终叠加剖面(如图2)；若残差大于阈值，返回步骤101，对24道地震原始数据重新进行速度扫描，得到修正后的速度谱，并进行后续操作，步骤开始循环，直至残差小于阈值；
[0077] 步骤105，对由步骤104确定的叠加剖面进行三瞬(瞬时频率、瞬时振幅、瞬时相位)属性提取，最终得到瞬时频率剖面(图3a)，瞬时振幅剖面(图3b)以及瞬时相位剖面(图3c)：
[0078] 利用连续小波变换时频分析方法，对叠加剖面提取瞬时频率，瞬时振幅以及瞬时相位属性；
[0079] 具体步骤如下：
[0080] 连续小波变换具备时域和频域双重良好的局部性和随尺度变化的自动调焦功能，而自动调焦可以通过母小波的伸缩和平移来实现；
[0081] 首先设母小波函数ψ(t)(可根据效果选取具体小波)，Ψ(ω)是母小波ψ(t)的Fourier变换，且满足关系式：
[0082]
[0083] 对母小波ψ(t)做平移和伸缩，得到一簇小波函数ψa,b(t)；
[0084]
[0085] 其中，a为尺度参数，b为平移参数；
[0086] 设叠加剖面CDP道为f(t)，其连续小波变换的定义为
[0087]
[0088] 其中，Wf(a,b)是上述叠后剖面CDP道f(t)的小波系数， 是ψa,b(t)的共轭函数，内积运算反映了信号的局部特征和小波函数的相似性，内积越大，相似性越高；
[0089] 在计算过程中，规定连续小波变换尺度变量的表达式为
[0090]
[0091] 其中，fw是母小波的中心频率如下式所示，scal是尺度的个数，s＝scal,scal-1,...1；
[0092]
[0093] 其中，ξ(f)是母小波ψ(t)的Fourier变换，f1和f2分别为母小波频带的上下限频率；
[0094] 通过上述小波变换的时频分析方法，计算得到a时刻最大的时
频谱能量|Wf(a,bmax)|，对应的尺度为bmax，形成瞬时振幅谱；
[0095] 对于每一组a时刻，通过上述a时刻的表达式
变形式，
[0096]
[0097] 计算得到fw,则fw为对应于该时刻的瞬时频率值，形成瞬时频率谱；
[0098] 给定一个噪音判别经验系数a0，当某一时刻时频谱能量的最大值小于分析信号的时频谱的平均值的a0倍时，判定该时刻的信号是随机噪音，规定其瞬时相位为零；反之，如果大于分析信号的时频谱的平均值的a0倍时，计算该时刻最大时频谱对应的瞬时相位，噪音判别条件和瞬时相位计算公式如下：
[0099] |Wf(a,bmax)|≤a0*mean(|Wf(a,b)|)
[0100]
[0101] 其中|Wf(a,bmax)|是a时刻最大的时频谱能量，对应的尺度是bmax，mean(|Wf(a,b)|)表示信号二维时频谱绝对值的平均值，Iang(a)是a时刻的瞬时相位，imag(Wf(a,bmax))是时刻a，尺度bmax对应的时频谱的
虚部，real(Wf(a,bmax))是时刻a，尺度bmax对应的时频谱的
实部；
[0102] 最终通过上述公式(1)获得瞬时相位谱；
[0103] 步骤106，将步骤100得到的抛弃式CTD原始资料作为拟测井数据，利用瞬时相位谱的同相轴连续性作为约束层位依据，利用叠后约束稀疏反演方法，最终得到海洋水体物性参数剖面，包括速度剖面(图4a)，温度剖面(图4b)，盐度剖面(图4c)。
[0104] 如以上所述的基于地震海洋学的海水物性测量方法，其特征在于在步骤100中，通过地震
数据采集仪器、抛弃式CTD，卫星导航以及
船舶的工作，得到24道地震原始数据，抛弃式CTD原始资料，以及导航数据。
[0105] 如以上所述的基于地震海洋学的海水物性测量方法，其特征在于在步骤104中，计算叠后剖面CDP道平均能量的具体步骤如下：
[0106] 基于Stein无偏似然估计的Sure阈值，对叠后剖面CDP道x(n)中的每个元素按照平方升序排列，得到新的信号序列，
[0107] W(n)＝sequence(x2),n＝1,…P
[0108]
选定W(n)中第i个元素的平方根为阈值λ，得到其
风险估计函数
[0109]
[0110] 再将风险估计函数求一次导数得到其最小风险点，记为imin，则其阈值为：
[0111] λ＝σW(imiin)0.5
[0112] 其中σ为噪声标准差。
[0113] 以上所述的基于地震海洋学的海水物性测量方法，其特征在于在步骤102中，曲波变换是一种保幅、多尺度精细图像数据处理方法，应用到地震海洋中，通过数据变换到共检波点域(CRP)，抓住各类噪声在不同数域中的特征，利用曲波变换进行干扰压制，达到原始资料的去噪效果最优化；
[0114] 设原始地震数据为D，有效信号部分为S，同理噪声部分为N，之间关系满足关系：
[0115] D＝S+N
[0116] 共检波点域曲波变换去噪方法步骤：
[0117] (1)根据不同数据解编方式，将原始数据共炮点域Dshot转换为共检波点域Drec；
[0118] (2)在各个数域中，根据数据大小确定尺度数j和方向数l，对数据进行曲波变换，将数据稀疏化，得到对应曲波系数Cshot(j,l,k)与Crec(j,l,k)；
[0119] (3)对稀疏系数进行阈值筛选，得到新的系数矩阵C'shot(j,l,k)与C'rec(j,l,k)，对其进行曲波反变换，得到重构地震数据D′shot与D′rec；
[0120] (4)通过差值原始数据与地震重构数据，得到各域一次变换后的去除的噪音成分N'shot与N'rec，将两者转换至同一数域后，对比分析两噪音成分，取两者交集得到噪声估计，认为这就是原始数据中的噪音N1；
[0121] (5)根据去噪后数据质量的反馈，确定方法返回迭代次数；
[0122] (6)叠加剖面继续去噪，分离噪音为N2。则认为原始噪音与分离噪音的关系为：N＝N1+N2
[0123] 以上所述的基于地震海洋学的海水物性测量方法，其特征在于在步骤106中，叠后约束反演是一种常见的适用于叠后数据反演的方法。
[0124] 稀疏脉冲反演方法，首先假设地层反射系数序列满足稀疏特性。通过计算出反射系数后，采用最大似然反演方法转换为波阻抗。其最小目标函数为：
[0125]
[0126] 式中，R2表示反射系数的均方值；N2表示噪音的均方值；r(k)表示第K个采样点的反射系数；n(k)表示第K个采样点的噪音；M表示反射层数；L表示采样总数；λ表示给定反射系数的似然值。然后通过多次迭代，求取反射系数。
[0127] 若利用最大似然反褶积获得的反射系数是ri，我们利用最大似然反演的算法转换反射系数序列得到宽带波阻抗。
[0128] 约束稀疏脉冲反演的关键是使下面这个目标函数得到优化：
[0129] E＝∑(ri)p+λq∑(di-si)q+α2∑(ti-zi)2
[0130] ri为反射系数；λ为控制稀疏脉冲因子，为实际地震记录与合成地震记录残差权重因子，与地震资料的信噪比的大小有关；di为地震数据；α为数据匹配的权因子；si为合成地震道数据；ti为井阻抗曲线的趋势；zi为实际的阻抗；p,q为L模因子，一般p＝1,q＝2；i地震道采样序号。
