一种基于稳健奇异谱分析的地震不规则噪声去除方法与流程

本发明是关于一种基于稳健奇异谱分析的地震不规则噪声去除方法，属于地震数据处理领域。

背景技术：

地震数据的噪声压制是地震数据处理中的一个重要流程.针对不同的地震干扰类型有不同的去噪方法，但最大的挑战仍然是如何做到保幅去噪，能否做到保幅去噪主要取决于对混杂在一起的信号和干扰进行有效分离的程度，调和分析方法通过诸如旋转、伸缩、平移等变换处理，将地震数据从原始的记录空间投影到不同的正交子空间，为研究人员提供了另一种分析和处理地震数据的视角，因此，如何针对地震数据的特点设计有针对性的特征函数或特征记录成为一个重要的研究方向。

不规则/或相干噪声衰减主要通过f-x和t-x域中的预测误差滤波器来进行，当然不规则噪声衰减也可以通过降秩方法实现，降秩方法认为地震波在地下介质中传播遵循波动规律，某一频率切片上的地震数据应呈现线性相关性，其对应的轨迹矩阵为低秩矩阵，而地震数据的高斯噪声和不规则噪声破坏了地震数据的相关性，导致地震轨迹矩阵的秩增加，因此，在高保真的条件下，地震数据降噪的过程就是对轨迹矩阵降秩的过程。降秩方法可以分为不同的类别，例如通过karhunen-loève变换过滤的特征图像滤波可以直接在t-x或f-x-y域的地震数据上进行降秩。

最近，对地震数据的噪声压制处理引入了奇异频谱分析方法(ssa)，也称为cadzow滤波，以减弱非相干噪声和重建地震数据。奇异频谱分析方法首先将地震数据嵌入hankel矩阵(汉克尔矩阵)，导致hankel矩阵的秩增加；然后通过降秩算法进行降噪，秩降低通常通过奇异值分解(svd)来实现，svd是一种非鲁棒矩阵分解方法，当存在不规则噪声时，常常导致次优结果。上述不规则噪声是指非高斯噪声，由已知或未知分布的大型孤立事件引入。采用稳健范数之后，基于奇异频谱分析的降秩滤波也用于抑制相干噪声，归纳起来，现有的基于降秩方法的地震不规则噪声去除方法主要为通过对hankel矩阵进行基于l2范数的截断svd方法，或采用基于稳健范数(例如turkey范数)的截断svd方法，其核心问题是基于l2范数的截断svd方法无法真正有效地去除不规则噪声，因为l2范数的核心本质就是求平均。基于turkey范数的截断svd方法中，turkey范数是凹函数，其优化问题本质上就是非凸优化，无法找到全局最优解以及智能获得局部最优解，且该方法依赖初始值，从而会影响降噪的性能。

技术实现要素：

针对上述问题，本发明的目的是提供一种能够智能获取全局最优解且不会影响降噪性能的基于稳健奇异谱分析的地震不规则噪声去除方法。

为实现上述目的，本发明采取以下技术方案：一种基于稳健奇异谱分析的地震不规则噪声去除方法，其特征在于，包括以下步骤：将单频地震数据嵌入hankel矩阵，构建单频地震数据的hankel矩阵；基于l1和l2混合范数，得到单频地震数据的hankel矩阵的目标函数；根据预先设定的所述目标函数中的正则化参数，对所述目标函数进行求解，得到基于l1和l2混合范数的hankel矩阵的奇异谱分析；对所述奇异谱分析依次进行交替最小化和加权最小二乘法最小化，求解得到单频地震数据的hankel矩阵的因子矩阵，进而得到单频地震数据的hankel矩阵的估计，完成单频地震数据的降噪。

进一步，所述目标函数为：

其中，λ1和λ2均表示正则化参数，λ1用于调节l1范数的权值大小，λ2用于调节l2范数的权值大小；m表示hankel矩阵；表示hankel矩阵的估计。

进一步，当正则化参数λ1＝0时，所述目标函数退化为：

所述目标函数即为l1范数下的拟合解，用于压制非高斯噪声；当正则化参数λ2＝0时，所述目标函数退化为：

所述目标函数即为最小二乘解，用于压制地震数据的高斯背景噪声。

进一步，所述奇异谱分析为：

其中，argminj表示退化后的目标函数；m表示hankel矩阵的行向量；u和v分别表示奇异谱分析的因子矩阵。

进一步，对所述奇异谱分析依次进行交替最小化和加权最小二乘法最小化，求解得到单频地震数据的hankel矩阵的因子矩阵，进而得到单频地震数据的hankel矩阵的估计，完成单频地震数据的降噪，具体为：对基于l1和l2混合范数下的hankel矩阵的奇异谱分析进行交替最小化：

其中，ew(v)表示自变量为因子矩阵v的目标函数，w表示权值矩阵，m^j表示hankel矩阵的第j列，且m＝(m¹，m²，…，mⁿ)＝(m1，m2，…，mn)^h，mn表示hankel矩阵的第n行的共轭转置；vj表示因子矩阵v的列向量；ew(u)表示自变量为因子矩阵u的目标函数；uj表示因子矩阵u的列向量；因子矩阵u和v为：

u＝(u¹，u²，…，u^k)＝(u1，u2，…，um)^h

v＝(v¹，v²，…，v^k)＝(v1，v2，…，vn)^h

其中，u^k表示因子矩阵u的列向量，um表示因子矩阵u的行向量，v^k表示因子矩阵v的列向量，vn表示因子矩阵v的行向量；对进行交替最小化后的奇异谱分析进行加权最小二乘法最小化：

其中，mi表示hankel矩阵的第i行的共轭转置；采用qr分解，迭代求解上述公式中的加权最小二乘法最小化问题，得到单频地震数据的hankel矩阵的因子矩阵u和v，进而得到单频地震数据的hankel矩阵的估计完成单频地震数据的降噪。

本发明由于采取以上技术方案，其具有以下优点：1、本发明基于l1和l2混合范数，代替截断svd采用的二次误差准则函数，对单频地震数据进行降噪，其中，l2范数能够对高斯噪声进行降噪，l1范数能够对非高斯规则噪声进行降噪，从而能够有效地同时去除高斯噪声和非高斯噪声。2、由于l1范数和l2范数都是凸函数，采用l1和l2混合范数得到的目标函数也是凸函数，而最优化问题本质上是凸优化问题，容易获得全局最优解，因此，本发明的方法能够获得更好的降噪性能，可以广泛应用于地震资料处理领域中。

附图说明

图1是采用本发明方法对地震数据的hankel矩阵进行降噪的对比示意图，其中，图1(a)表示干净地震数据的hankel矩阵，图1(b)表示含噪地震数据的hankel矩阵，图1(c)表示降噪后地震数据的hankel矩阵；

图2是原始的叠后地震剖面示意图；

图3是采用本发明方法对图2中叠后地震剖面进行降噪后的示意图；

图4是采用本发明方法对图2中叠后地震剖面去除的干扰的示意图；

图5是原始的叠后地震剖面示意图；

图6是采用本发明方法对图5中叠后地震剖面进行降噪后的示意图；

图7是采用本发明方法对图5中叠后地震剖面去除的干扰的示意图。

具体实施方式

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

小窗口内的地震数据dj(ω)通过频率-空间域内的平面波叠加得到：

其中，i和j表示空间轴上的地震道集索引值，且n表示地震道数目；ω表示时间频率，在上述公式中，假设地震数据dj(ω)由具有不同射线参数pk的k个线性事件(即平面波)组成，ak(ω)表示第k个平面波的复振幅，δx表示地震剖面之间的空间间隔，k表示线性事件的序数。

奇异谱分析方法是将地震数据dj(ω)变换后取单独频率的数据得到单频地震数据d(ω)＝(d1(ω)，d2(ω），…，dn(ω))^t，并将单频地震数据嵌入信号，构建轨迹矩阵，即hankel矩阵m(ω)：

其中，表示hankel算子，为方便起见，选择hankel矩阵的行向量数目使hankel矩阵逼近为方阵表示复数域。省略符号ω，对所有频率进行分析，并对k个平面波进行叠加得到hankel矩阵的秩rank(m)＝k。外加噪声d将增加矩阵m的秩k，需要通过降秩来减小外加噪声。

因此，奇异谱分析通过以下表达式表示：

其中，表示反对角平均算子；表示降秩操作符，通过秩k逼近hankel矩阵m；表示hankel算子；运算符通过平均对角线将hankel形式转换成向量形式，对于多维信号则采用块hankel矩阵和阻尼对角平均算子。如图1所示，降秩操作可以通过截断svd、随机svd或采用lanczos双对角化方法，以及高效矩阵向量乘法的快速算法来实现，其中，快速算法采用快速傅里叶变换。

hankel矩阵m的秩k可以通过求解下式近似得到：

其中，表示frobenius范数；为矩阵m和n表示矩阵x的大小，xij表示矩阵x的元素；公式(4)具有全局唯一的局部最小值，这个解具有解析表达式，由截断奇异值分解(tsvd)给出：

其中，uk和vk表示包含奇异值向量的矩阵与秩k矩阵最大奇异值σq相关的向量，且q表示特征值或奇异值序号，q＝1，…，k；∑k表示特征值矩阵(即奇异值)，这些奇异值也是矩阵(表示实数域)的对角元素，称为eckart-young定理。即使公式(4)的解决方案是简单的，二次匹配函数使得截断奇异值分解对非高斯噪声相当敏感，这个缺点限制了奇异谱分析方法在其中包含异常值的应用。

基于上述现有奇异谱分析存在的问题，本发明提供的基于稳健奇异谱分析的地震不规则噪声去除方法，包括以下步骤：

1)将单频地震数据嵌入hankel矩阵，构建单频地震数据的hankel矩阵m。

2)基于l1和l2混合范数，代替截断svd采用的二次误差准则函数，得到单频地震数据的hankel矩阵m的目标函数，具体为：

如图1所示，传统的基于l2范数的奇异谱分析只能压制噪声中的高斯部分，而基于l1范数的奇异谱分析对呈现出厚拖尾的非高斯部分具有更好的鲁棒性，因此，结合l1范数与l2范数在不同噪声环境下的优势，得到基于l1和l2混合范数的hankel矩阵m的目标函数j：

其中，λ1和λ2均表示正则化参数，λ1用于调节l1范数的权值大小，λ2用于调节l2范数的权值大小，表示hankel矩阵m的估计。

3)设定单频地震数据的hankel矩阵m的目标函数中的正则化参数λ1和λ2，其中，上述正则化参数λ1和λ2可以根据实际情况进行设定。

当λ1＝0时，上述目标函数(6)退化为：

公式(7)即为l1范数下的拟合解，用于压制非高斯噪声，噪声的概率密度函数可用拉普拉斯分布描述。

当λ2＝0时，上述目标函数(6)退化为：

公式(8)即为最小二乘解，用于压制地震数据的高斯背景噪声。

4)根据设定的正则化参数λ1和λ2，对单频地震数据的hankel矩阵m的目标函数进行优化求解，得到基于l1和l2混合范数的hankel矩阵m的奇异谱分析，即特征向量矩阵(u，v)：

其中，argminj表示退化后的目标函数；m表示hankel矩阵m的行向量；u和v分别表示特征向量矩阵的因子矩阵，且

5)对基l1和l2混合范数下的hankel矩阵m的奇异谱分析依次进行交替最小化(admm)和加权最小二乘法最小化(irls)，求解得到单频地震数据的hankel矩阵的因子矩阵u和v，进而得到单频地震数据的hankel矩阵的稳健估计，完成单频地震数据的降噪，具体为：

5.1)对基于l1和l2混合范数下的hankel矩阵m的奇异谱分析进行交替最小化：

其中，ew(v)表示自变量为因子矩阵v的目标函数，w表示权值矩阵，m^j表示hankel矩阵m的第j列，且m＝(m¹，m²，…，mⁿ)＝(m1，m2，…，mn)^h，mn表示hankel矩阵m的第n行的共轭转置；vj表示因子矩阵v的列向量；ew(u)表示自变量为因子矩阵u的目标函数；uj表示因子矩阵u的列向量；因子矩阵u和v可以表示为：

u＝(u¹，u²，…，u^k)＝(u1，u2，…，um)^h(12)

v＝(v¹，v²，…，v^k)＝(v1，v2，…，vn)^h(13)

其中，u^k表示因子矩阵u的列向量，um表示因子矩阵u的行向量，v^k表示因子矩阵v的列向量，vn表示因子矩阵v的行向量。

5.2)对进行交替最小化后的奇异谱分析进行加权最小二乘法最小化，采用qr分解，迭代求解直到收敛或到达最大迭代步数，得到单频地震数据的hankel矩阵m的因子矩阵u和v，进而得到单频地震数据的hankel矩阵的稳健估计，完成单频地震数据的降噪。

公式(10)和(11)可以分解成更小的加权最小二乘法最小化问题，即将公式(10)和(11)进行加权最小二乘法最小化：

其中，mi表示hankel矩阵m的第i行的共轭转置。

采用qr分解，迭代求解上述公式(14)和(15)中的加权最小二乘法最小化问题，得到单频地震数据的hankel矩阵m的因子矩阵u和v，进而得到单频地震数据的hankel矩阵的稳健估计完成单频地震数据的降噪。

如图2～7所示，通过对比两组原始的叠后地震剖面，采用本发明的地震不规则噪声去除方法进行降噪后的叠后地震剖面，以及采用本发明的地震不规则噪声去除方法对叠后地震剖面去除的干扰，可以看出本发明的地震不规则噪声去除方法能够有效去除地震数据中的高斯和不规则噪声，且地震数据本身没有受到损害。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。