本发明属于地震数据处理领域，特别涉及一种地震去噪技术。
背景技术：
随着油气勘探开发的不断深入，勘探地区的地质构造和地质环境愈发复杂多变。在陆上和海上中深层勘探目标区采集的地震数据中，受复杂地质条件和其他因素等影响，原始资料的品质通常比较差，主要表现为：剖面中存在着较严重的随机噪声和相干噪声干扰，有效反射信号较弱且被强噪声覆盖，资料信噪比和保真度较低。由于地震数据的信噪比会直接影响到地质解释的精度，因此地震数据去噪处理就成了十分重要且基础的环节，也有必要对地震去噪方法展开深入研究。
基于矩阵降秩理论的地震数据去噪方法是近年来国内外去噪研究领域的热点，它的研究热点集中在降秩约束算法、自适应秩选择方法、混合约束模型等方面。矩阵降秩类方法广泛应用于含噪和缺失矩阵恢复、地震数据去噪与插值等任务，低秩逼近理论将由线性同相轴构成的纯净地震数据投影到hankel矩阵或toeplitz矩阵后，此时矩阵具有低秩性质，矩阵的秩大小与数据中指数相关信号数目有关，而加性随机噪声的存在会使该矩阵出现满秩现象，通过对矩阵实施降秩处理就可以去除随机噪声。
最近，引入了多道奇异谱分析(mssa)方法，该方法是一种有效的地震随机噪声衰减算法，它通过截断奇异值分解(tsvd)将含噪信号的hankel矩阵向量空间分解为信号子空间和噪声子空间，比传统的基于奇异值分解(svd)的奇异谱分析(ssa)方法具有更好的去噪效果。然而，tsvd只是通过截断部分较小的奇异值来达到降秩去噪的目的，无法将含噪地震数据很好地分解成噪声子空间和信号子空间，在低信噪比的情况下效果不是很理想。
技术实现要素：
为解决上述技术问题，本发明提出一种奇异值衰减的降秩去噪方法，在低信噪比的情况下也能取得较好的效果，在地震随机噪声压制应用中有优异表现。
本发明采用的技术方案为：一种奇异值衰减的降秩去噪方法，包括：
s1、对原始地震数据进行截断奇异值分解，得到含噪观测矩阵；
s2、计算出含噪观测矩阵的左、右奇异向量的优化权值系数；
s3、根据优化权值系数仅依赖于纯噪声矩阵的极限奇异值分布这一特性，将含噪信号分解为噪声子空间和信号子空间；
s4、添加正则化算子来约束奇异值，得到低秩估计；
s5、利用反对角线平均处理，从而获得频率域的地震数据的稳健估计。
进一步地，步骤s1所述的含噪观测矩阵表达式为：
其中，
进一步地，步骤s2所述优化权值系数βopt的表达式为：
其中，
更进一步地，优化权值系数βopt的求解过程为：
a1、令ur＝[u1…ur]，vr＝[v1…vr]，θr＝[θ1…θr]，
a2、根据βi*的表达式将优化权值系数βopt的表达式转换为：
其中，
进一步地，步骤s4所述正则化算子表达式为：
其中，i是单位矩阵，p是控制正则化算子zi的正则化因子。
更进一步地，p取值为2～5。
本发明的有益效果：本发明的一种基于正则化奇异值衰减的降秩去噪方法，首先计算出含噪观测矩阵的左、右奇异向量的优化权值系数，利用优化权值系数仅依赖于纯噪声矩阵的极限奇异值分布这一特性，能够较好地将含噪信号分解为噪声子空间和信号子空间，然后在其基础上添加正则化算子来约束奇异值，以获得更稳健的低秩估计，最后利用反对角线平均处理，从而获得频率域的地震数据的稳健估计，在压制地震随机噪声的应用中具有突出的去噪能力。本发明通过实际数据测试了我们基于正则化奇异值衰减的降秩去噪方法，实验结果表明该方法在去噪性能上比传统mssa降秩约束方法有较大提升，在低信噪比的情况下也能取得较好的效果，在地震随机噪声压制应用中有优异表现。
附图说明
图1为本发明的方法流程图；
图2为本发明实施例提供的hankel矩阵构造示意图；
其中，图2(a)为频率分量数据，图2(b)为hankel矩阵，图2(c)为块hankel矩阵；
图3为本发明实施例提供的svd算法矩阵分解示意图；
图4为本发明实施例提供的mssa地震数据去噪示意图；
其中，图4(a)为原始地震数据(snr＝-10.859db)，图4(b)为mssa地震去噪结果(snr＝0.664db)，图4(c)为mssa去除的噪声；
图5为本发明实施例提供的奇异值衰减地震数据去噪示意图；
其中，图5(a)为原始地震数据(snr＝-10.859db)，图5(b)为奇异值衰减地震去噪结果(snr＝3.736db)，图5(c)为奇异值衰减去除的噪声；
图6为本发明实施例提供的原始叠前含噪地震剖面；
图7为本发明实施例提供的降噪后的叠前地震剖面；
图8为本发明实施例提供的去除的噪声；
图9为本发明实施例提供的原始叠前含噪地震剖面；
图10为本发明实施例提供的降噪后的叠前地震剖面；
图11为本发明实施例提供的去除的噪声。
具体实施方式
为便于本领域技术人员理解本发明的技术内容，下面结合附图1-11对本发明内容进一步阐释。
如图1所示，本发明的一种奇异值衰减的降秩去噪方法，包括：
s1、对原始地震数据进行截断奇异值分解，得到含噪观测矩阵；
s2、计算出含噪观测矩阵的左、右奇异向量的优化权值系数；
s3、根据优化权值系数仅依赖于纯噪声矩阵的极限奇异值分布这一特性，将含噪信号分解为噪声子空间和信号子空间；
s4、添加正则化算子来约束奇异值，得到低秩估计；
s5、利用反对角线平均处理，从而获得频率域的地震数据的稳健估计。
步骤s1具体实现过程为：基于tsvd(truncatedsingularvaluedecomposition，截断奇异值分解)的多道奇异谱分析方法，得到含噪观测矩阵：
如图2所示为hankel矩阵构造示意图，其中图2(a)为频率分量数据；图2(b)为hankel矩阵；图2(c)为块hankel矩阵。图2中frequencyslice表示频率分量，hankelmatrix表示hankel矩阵，blockhankelmatrix表示hankel矩阵块。
对于大小为t＝1…nt，x＝1…nx，y＝1…ny的三维地震数据体d(x,y,t)，首先需要将原始地震数据通过傅里叶变换由时间域转到频率域：
其中，nt表示地震数据体在时间轴的分量，nx表示地震数据体在x轴的分量，ny表示地震数据体在y轴的分量，d(x,y,t)表示时域的原始数据体，d(x,y,ω)表示经傅里叶变换由时间域转换到频率域的地震数据体，ω＝1…nω，nω表示频率域的地震数据体的频率分量。每个特定频率分量ω0可以表示为矩阵：
d(ny,nx)表示在(x,y)的傅里叶变换。
对该频率切片分量d(ω0)中每行构造hankel矩阵，例如d(ω0)第i行构造的hankel矩阵ri可以表示为：
通常需要将ri构造为方阵或近似方阵形式，即
矩阵h同样也要构造为方阵或近似方阵形式，因此
h＝s n(5)
其中，s表示h中的纯净信号分量，n表示h中的随机噪声分量。若h和n都满秩，即rank(h)＝rank(n)＝q，而s的秩为r，即rank(s)＝r＜q，如图3所示，对h作奇异值分解(singularvaluedecomposition，svd)可以表示为：
其中，
mssa采用基于eckart-young-mirsky(eym)原理的截断奇异值分解tsvd方法，通过将
从图4中实际地震数据测试可以看出，mssa得到的结果仍然有比较多的噪音，这是因为mssa无法将地震数据很好地分解成噪声子空间和信号子空间。因此我们需要引入奇异值衰减的降秩方法来解决这一问题。图4中inline表示主测线地震解释剖面，xline表示联络线地震解释剖面，tline表示时间。
步骤s2具体为：
含噪地震数据通过mssa方法，无法完全使信号和噪声分离，因此我们将采用奇异值衰减的降秩方法来解决这一问题。我们将式子(6)转化为最优化问题求解形式，通过求解βi来逼近纯净信号s：
其中，
通过式子(9)中优化权值参数βi*会在大维矩阵极限下收敛到极限噪声分布μm的特定积分变换(d-变换)，可以得到奇异值衰减求解优化权值系数βopt的解析表达式如下所示：
其中，
利用式子(9)求得优化权值参数βopt与其对应的左右向量可以得到降秩去噪后的hankel矩阵
如图5所示为奇异值衰减地震数据去噪示意图，图5(a)为原始地震数据(snr＝-10.859db)；图5(b)为奇异值衰减地震去噪结果(snr＝3.736db)；图5(c)为奇异值衰减去除的噪声；从图5中实际地震数据测试可以看出，由于奇异值衰减法通过计算左右奇异向量的优化权值来低秩逼近信号分量的奇异值，在数值上表现为对原始奇异值进行了衰减处理，从优化奇异向量权值的角度对奇异值进行了优化约束处理，因此该方法在简单的倾斜构造处损失信号相对较少，它的去噪结果相比于传统的mssa结果在残余噪声控制方面有不错的提升。但是在断层和不整合面等复杂构造处损失信号相对较多。因此我们需要引入正则化算子来解决这一问题。图5中inline表示主测线地震解释剖面，xline表示联络线地震解释剖面，tline表示时间。
步骤s4具体为：
地震数据中不仅仅有简单的倾斜构造这些简单构造，还有断层和不整合面等复杂构造，而由于奇异值衰减方法无法对复杂构造进行很好的去噪，因此本专利还将在奇异值衰减方法的基础上，引入正则化算子z，它可以在计算优化权值系数γi的过程中对奇异值再进行软阈值处理，能够对复杂构造处的噪声进行稳健的去噪处理，从而取得更好更稳健的去噪效果。
我们将正则化算子z引入式(7)，得到：
我们经过大量数值实验，发现当正则化算子取到如下情况时，能得到一个较好的去噪结果，其表达式如下：
其中，i是单位矩阵，p是控制正则化算子zi的正则化因子，根据本发明的研究发现，一般取2～5之间。
由式(10)、(12)和(13)得到最终基于正则化奇异值衰减的降秩去噪后的hankel矩阵
步骤s5具体为：
本发明核心是设计一种基于正则化奇异值衰减的降秩去噪方法来压制随机噪声，同时稳健的反问题求解算法也是非常重要的。本发明使用反对角线平均方法来将降秩后的hankel矩阵恢复为频域地震数据。假设s为恢复的频域数据，若降秩处理后矩阵
如图6所示为原始叠前含噪地震剖面，如图7所示为降噪后的叠前地震剖面，如图8所示为去除的噪声，如图9所示为原始叠前含噪地震剖面，如图10所示为降噪后的叠前地震剖面，如图11所示为去除的噪声，从图6-图11结果可以看出，实际地震剖面去噪结果符合理论假设，并且实际效果非常好，降噪效果很明显。图6-图8中trace表示记录道，time表示时间，单位为s；图9-图11中tracenumber表示记录道编号，time表示时间，单位为s。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。
1.一种奇异值衰减的降秩去噪方法，其特征在于，包括：
s1、对原始地震数据进行截断奇异值分解，得到含噪观测矩阵；
s2、计算出含噪观测矩阵的左、右奇异向量的优化权值系数；
s3、根据优化权值系数仅依赖于纯噪声矩阵的极限奇异值分布这一特性，将含噪信号分解为噪声子空间和信号子空间；
s4、添加正则化算子来约束奇异值，得到低秩估计；
s5、利用反对角线平均处理，从而获得频率域的地震数据的稳健估计。
2.根据权利要求1所述的一种奇异值衰减的降秩去噪方法，其特征在于，步骤s1所述的含噪观测矩阵表达式为：
其中，
3.根据权利要求2所述的一种奇异值衰减的降秩去噪方法，其特征在于，步骤s2所述优化权值系数βopt的表达式为：
其中，
4.根据权利要求5所述的一种奇异值衰减的降秩去噪方法，其特征在于，优化权值系数βopt的求解过程为：
a1、令ur＝[u1…ur]，vr＝[v1…vr]，θr＝[θ1…θr]，
a2、根据βi*的表达式将优化权值系数βopt的表达式转换为：
其中，
5.根据权利要求3所述的一种奇异值衰减的降秩去噪方法，其特征在于，步骤s4所述正则化算子表达式为：
其中，i是单位矩阵，p是控制正则化算子zi的正则化因子。
6.根据权利要求5所述的一种奇异值衰减的降秩去噪方法，其特征在于，p取值为2～5。
技术总结