一种基于地震资料共振成像中的时频分析方法

摘要

本发明涉及一种基于地震资料共振成像中的时频分析方法。该方法基于稀疏优化算法和Tikhonov正则化理论，通过引入拉格朗日乘子和控制步长，采用Newton迭代计算的方式获取频率系数的最优解。该方法提高了成像分析的速度、以及成像的分辨率和准确性，该方法能有效的应用于地下空间精细成像中。

说明书

技术领域

本发明涉及勘探技术领域，具体地说，涉及一种基于地震资料共振成像中的时频分析方法。

背景技术

自然界中，任何物质都有其自身的固有频率，包括地下各地质体。当有一个宽频的震动传播到该地质体，特征固有频率能量将被放大，通过观测被放大的特征频率信号，对特征频率信号成像获得地下精细成像效果。目前，共振成像法应用广泛，例如大型土木工程方面的查明地下介质结构并进行分层，评价地基的加固效果，建筑工程的抗震设计，公路、机场跑道的质量无损检测；环境地质方面的勘查岩溶洞穴的分布范围，滑坡活动面等地质灾害和环境地质的调查；地下空间精细探测；以及考古发现等。

基于地震频谱的动力学特征，1989年，Nakamura首次提出了基于地震资料的单点谱比法。该方法其中的关键技术之一就是地震数据的频谱分析。常用的地震数据频谱分析方法主要包括傅里叶变换(FFT)、短时傅里叶变换(STFT)、连续小波分析(CWT)、匹配追踪算法(MP)、以及魏格纳分布(WVD)等。其中，FFT是最早提出并用于信号频谱分析的技术，但由于其是对数据进行整体分析，缺乏时间分辨率，不适用于非平稳的地震信号分析；Potter(1947)对FFT加上窗函数，首次提出STFT的概念，进而能获得时间-频域信息，但由于每次计算的窗函数固定，不能兼顾时间分辨率和频率分辨率；Morlet(1982)提出的CWT是一种多分辨率的时频分析方法，但是计算量大，耗时长；Mallet等(1999)提出的MP方法是借助子波字典获得由匹配参数构成的视频谱，该方法十分依赖于所选择的子波字典，并且计算量大；WVD方法弱化窗函数的影响，提高了时间分辨率，但由于其是非线性变换会造成交叉项干扰是的频谱失真。

发明内容

在地震数据时频分析过程中，噪音干扰，资料品质差、先验信息缺乏、数据量不足等条件都会导致资料时频分析反演时出现解不存在、解不唯一、或者解不稳定等问题。为了克服地震面波数据时频分析反问题中可能出现的不适当性、计算过程的病态性以及大规模计算等问题。提出了一种基于稀疏优化算法和正则化理论的时频分析中地震频率反问题求解算法。

基于地震资料共振成像中的时频分析方法的具体步骤如下：

步骤1：设置初始频率系数f₁，f₁>0；迭代容忍度ε_f>0，ε_λ>0，ε_μ>0，初始控制步长τ₁＝1，迭代步长τ的修正因子θ∈(0,1)，衰减参数β∈[0,1]，迭代次数阈值i_max＝30；

步骤2：初始化(f₁,λ₁,s₁)^T，其中，拉格朗日乘子λ₁>0，中间参数s₁>0，迭代次数i＝1；

步骤3：计算余量：其中M为实数或者复数正弦基函数的核矩阵，d为带噪音的地震面波数据，N是数据d的长度，e为一个所有分量都为1的向量，μ为障碍函数参数，若或i＝i_max，结束迭代，否则，进入步骤4；

步骤4：取衰减参数β∈(0,1]，计算参数矫正量(Δf_i,Δλ_i,Δs_i)；

步骤5：计算取θ∈(0,1)，令τ_i+1＝min{θτ^max,1}，更新迭代，i＝i+1，返回步骤3，继续判断。

其中，步骤2中所述的拉格朗日乘子满足以下条件，其中，μ为障碍函数参数。

其中，障碍函数参数μ满足以下条件，

其中，步骤2中所述的中间参数s＝μD―1e，其中μ为障碍函数参数，d为带噪音的地震面波数据，N是数据d的长度，e为一个所有分量都为1的向量。

其中，步骤计4中所述的参数矫正量(Δf_i,Δλ_i,Δs_i)由以下公式获得其中

该方法基于稀疏优化算法和Tikhonov正则化理论，通过引入拉格朗日乘子和控制步长，采用Newton迭代计算的方式获取频率系数的最优解。该方法提高了成像分析的速度、以及成像的分辨率和准确性，该方法能有效的应用于地下空间精细成像中。

附图说明

图1：模拟数据分析图

图2：预处理后数据整体分布图

图3：桩号34点预处理后数据分布图

图4：桩号34点资料的频谱分析结果

图5：共振成像剖面与钻井资料的对比分析图

图6：共振成像剖面

图7：共振成像剖面的局部分析图(管道)

图8：共振成像剖面的局部分析图(塌陷)

具体实施方式

我们可以把地震面波时频分析问题归结为如下方程形式：

Mf＝d(1)

其中，f表示频率系数，d表示带噪音的地震面波数据，M表示实数或者复数正弦基函数的核矩阵，

M(t,x)＝cos(2π·k·Δx·n·Δt)+isin(2π·k·Δx·n·Δt)(2)

其中，k表示频率域采样指数，Δx表示频率增量，n表示时间采样指数，Δt表示时间增量。

对于上式的求解一般可通过如下公式(3)、(4)求得：

M^*Mf＝M^*f(3)

f＝(M^*M)^―1M^*d(4)

如果M为一正交矩阵，则M^*M＝I，这时公式(1)的最小二乘解可由f＝M^*d获得，即离散傅里叶变换(DFT)。但是，在实际地震资料中，M的元素通常是相关的，M具有非正交性。这将导致公式(1)解的不唯一，以及直接求解公式(4)的病态性。

考虑地震观测系统的限制，地球物理反问题一般是离散的、稀疏的，以及频率系数f的非负性等特征。为此，本文利用Tikhonov正则化思想，对地震时频分析中求解频率f的反问题建立如下基于l1范数的稀疏优化模型：

公式(5)的原始对偶形式为

min e^Tf,s.t.Mf＝d,f≥0(6)

其中，f有N个分量，f＝[f₁,f₂,…f_N]；N是数据d的长度；e是一个所有分量都为1的向量。

利用障碍函数(μ>0)来求解公式(6)的原始对偶问题：

公式(7)中，当f<0时，表达式无意义；也就是借助障碍函数将不等式约束(公式6)转化为无约束的最小化问题(公式7)。当对偶间隙μ→0，公式(7)的解会不断逼近原问题的全局最优解。

引入拉格朗日乘子λ，利用拉格朗日乘子法对公式(7)的最小化问题进行求解，

由于

grad_fL(f,λ)＝e―μD^―1e-M^Tλ(9)

grad_λL(f,λ)＝d―Mf(10)

令s＝μD^―1e，可得到公式(7)的最优性条件：

s+M^Tλ＝e,Mf＝d,Ds＝μe,f>0(11)

其中，

该系统的一个最优解(f,λ,s)^T即是原始的可行解又是对偶可行的。在线性规划领域，向量(f,λ,s)^T称为中心路径。

考虑到对偶间隙

e^Tf―d^Tλ＝f^Ts＝μN(13)

因此，数值上，求解非线性方程组

其中，β是衰减参数，β∈[0,1]，控制下降速度使得迭代点的轨迹在中心路径附近。

借助Newton迭代法对上式进行求解，假设第i次迭代的(f_i,λ_i,s_i)^T给定，对公式(14)求导，其梯度如下：

其中，

Newton迭代方法中的参数校正量(Δf_i,Δλ_i,Δs_i)^T可由以下公式获得：

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

利用图1中的模拟信号进行测试分析，令ε_f＝0.0001，ε_λ＝0.0001，ε_μ＝0.00001，频率系数均为0.1，μ＝4，θ＝0.85，β＝0.9，s₁可由s＝μD^―1e获得，λ₁可由公式(9)获得。对比常规STFT，本算法在不降低计算效率的情况下，能清晰的区分10-15Hz段中的11.108Hz、12.573Hz、13.062Hz三个频率信息，改善了资料的频域分辨率，证实算法的有效性。进一步，根据浅层勘探(地下200米以上)中共振频率-深度的经验关系，，A为一常数，取值范围是50-300。以A＝180为例，三个频率分别对应地下8.499米、7.264米、6.921米，也因此表明该算法能清晰区分、并刻画此地下三个不同深度的异常体，证实了该算法在浅层勘探共振成像中的有效性。

研究区位于成都市空港区牧华路沿线，共振激发方式为天然源，采用的是高精度2.5HZ三分量低频检波器，采集时间长度为18分钟，采样率为4毫秒，图2预处理后数据整体分布图，图3是桩号34点预处理后数据分布图。

同样，在本研究区的数据测试中，令ε_f＝0.0001，ε_λ＝0.0001，ε_μ＝0.00001，频率系数均为0.1，μ＝4，θ＝0.85，β＝0.9，s₁可由s＝μD―1e获得，λ₁可由公式(9)获得。对预处理之后的数据，进行逐点单分量的频谱分析。图4是桩号34点处地震资料的频谱分析结果。由图5分析可见，本次研究的共振成像剖面与相应的钻井资料对比，成像深度误差小于0.5米，证实了本方法结果的可靠性与准确性。

从成都空港区的共振成像剖面(图6)可见清晰的鹅卵石层、粘土层、强风化基岩和弱风化基岩等地质现象，为该研究区地下地质空间建模提供高精度成像。同时，从共振成像局部剖面(图7、8)的精细对比分析可见沿测线的下水管道、雨水管道、塌陷等，成像刻画清晰，且成像分辨率较高，这可为城市地下管线探测，特别是深埋管线探测，以及道路灾害探测，提供高精度成像。此外，从分析结果可推断地面下陷原因与下覆雨水管道掏空造成地层疏松有关，这可以为城市建设规划提供重要的参考。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。