本发明涉及一种地球物理勘探中的地震数据处理方法，尤其是基于pca-emd的并行震源地震数据随机噪声压制方法。
背景技术：
为了提高效率和降低成本，现在的地震勘探技术正在从单震源地震勘探方式向并行地震勘探方式发展。但并行震源地震勘探在采集过程中往往受到随机噪声的影响，采集到的并行地震数据中常伴有随机噪声，因此得到的地震数质量不高，进而影响了后期地震数据处理解释和偏移成像质量。并行震源地震数据涉及空间范围广，时间跨度大，造成不同地震道的随机噪声在能量强弱和稳定性方面差异性大，同时这种空间范围广的特点导致随机噪声的频率分布特性随着坏境的不同而发生变化，而常规单震源勘探方式受到的随机噪声在能量强弱、平稳性和噪声频率分布方面相对稳定，这两种勘探方式受到的随机噪声显著的差异使得常规单震源地震数据压噪方法无法适用于并行震源地震数据噪声压制。目前针对并行震源采集方式一般只能通过增加覆盖次数提高信噪比，但该方法会显著增加采集成本。而专门研究并行震源数据随机噪声压制方法的文献尚少见，在其它领域对随机噪声压制的方法主要有频域法和时域法两大类，频域法如小波滤波，维纳滤波等都在频率域进行噪声压制，虽然该类方法简单，但当有用信号频带与噪声频带混叠时，这类方法易损害有用信号细节。时域法如svd、k-l等方法虽然可有效保护有用信号细节，但该类方法只适应于随机噪声能量较弱条件。可见，上述方法都不适应于并行震源数据随机噪声的压制。
技术实现要素：
本发明的目的就在于针对上述现有技术的不足，提供一种基于pca-emd的并行震源地震数据随机噪声压制方法。
本发明的思想是：并行震源地震勘探技术提高了工作效率和降低了生产成本，但采集得到的地震数据往往受到随机噪声的影响，这样获得的地震数据经常影响了后期地震数据处理解释和偏移成像质量，本发明首先通过信噪间频谱特征确定经emd分解得到的有用信号占主导地位的模态分量并进行重构，将重构结果根据相空间理论构造hankel矩阵，并对其进行主成分分解与有效信号的恢复，从而实现了并行震源地震勘探数据随机噪声的压制。
本发明的目的是通过以下技术方案实现的：
一种基于pca-emd的并行震源地震数据随机噪声压制方法，包括以下步骤：
a、对并行震源数据单道信号x(l)进行频谱分析，估计有用信号频谱范围，其中l为采样序列，l＝1,2,…,n，n为最大采样点；
b、对x(l)进行emd分解，得到若干模态分量及余项，x(l)通过下式进行emd分解：
其中imfk为模态分量中第k个模态分量，k＝1,2,…,k，k为模态分量总数，r为余项；
c、按imf1～imfk的顺序分别进行频谱分析得到对应的频谱，若首次出现第s个模态分量(imfs)的频谱主要位于有用信号频谱范围内，则imfs～imfk为有用信号占主导地位的模态分量，其中s≤k；
d、将imfs～imfk及余项进行重构得到重构结果x'(l)，如公式
x'(l)＝imfs+imfs+1+…+imfk+r(2)
e、据相空间重构理论，对x'(l)构造hankel矩阵
该矩阵的行数记为m，列数记为n，记m＝n-n+1，若n为偶数，则令m＝n/2+1，n＝n/2，若n为奇数，则令m＝(n+1)/2，n＝(n+1)/2；
f、计算h的协方差矩阵γ，如公式
其中ht为h的转置矩阵，“·”表示矩阵乘法；
g、用奇异值分解法，计算协方差矩阵γ的特征值矩阵λ和特征向量矩阵r，则存在公式
γ＝r·λ·rt(5)
其中λ为特征值由大到小排列的对角矩阵，λ＝diag[λ1,λ2,…,λn]，λ1,λ2,…,λn为特征值，r为各个特征值对应的特征向量矩阵，rt为r的转置矩阵，且满足rt·r＝r·rt＝e，其中e为单位矩阵；
h、h经线性映射，得到主成分矩阵φ，如公式
φ＝rt·h(6)
i、计算前p个特征值累计贡献率：
其中λj为特征值，p为所取特征值个数，1≤p≤n，j＝1,2,…,n；
j、若满足
y＝r·φ'(8)
则y为压制随机噪声后hankel重构矩阵，其具体形式记为
h、令x”(l)＝[y(1),y(2),···,y(n)]，则x”(l)即为对应x(l)的压制随机噪声后的压噪信号。
有益效果：经试验，本发明公开的一种基于pca-emd的并行震源地震数据随机噪声压制方法能够实现在并行震源地震勘探数据中压制随机噪声，该算法处理数据快，相比于emd压制随机噪声方法，本方法能够在全频带范围压制随机噪声，并且在信噪频带混叠时，不仅能够压制噪声能量，还能有效保护信号细节。由于其良好的信噪比改善能力，使得处理后的目标数据定位误差更小，降低了数据处理成本，此外在强噪声条件下本方法更具优势。
附图说明：
图1单道信号及频谱，(a)为单道信号，(b)为单道信号频谱
图2部分模态分量及对应频谱，(a)为imf1，(b)为imf1频谱，(c)为imf2，(d)为imf2频谱
图3压噪处理前后结果对比，(a)压噪前，(b)为压噪后
具体实施方式：
下面结合附图和实施例对本发明做进一步的详细说明：
在本实施例中使用2个震源为一组的方法进行激发，记录时间为3s，采样率1000hz。
一种基于pca-emd的并行震源地震数据随机噪声压制方法，包括以下步骤：
a、对并行震源数据单道信号x(l)进行频谱分析，估计有用信号频谱范围，其中l为采样序列，l＝1,2,…,n，n为最大采样点，本例中l＝1,2,…,3001，n＝3001，有用信号频谱范围为0～100hz；
b、对x(l)进行emd分解，得到若干模态分量及余项，x(l)通过下式进行emd分解：
其中imfk为模态分量中第k个模态分量，k＝1,2,…,k，k为模态分量总数，r为余项，本例中k＝1,2,…,10，k＝10，r为余项；
c、按imf1～imfk的顺序分别进行频谱分析得到对应的频谱，若首次出现第s个模态分量(imfs)的频谱主要位于有用信号频谱范围内，则imfs～imfk为有用信号占主导地位的模态分量，其中s≤k，本例中有用信号占主导地位的模态分量为imf2～imf10；
d、将imfs～imfk及余项进行重构得到重构结果x'(l)，如公式
x'(l)＝imfs+imfs+1+…+imfk+r(2)
本例中x'(l)＝imf2+imf3+…+imf10+r；
e、据相空间重构理论，对x'(l)构造hankel矩阵
该矩阵的行数记为m，列数记为n，记m＝n-n+1，若n为偶数，则令m＝n/2+1，n＝n/2，若n为奇数，则令m＝(n+1)/2，n＝(n+1)/2，本例中m＝1501，n＝1501；
f、计算h的协方差矩阵γ，如公式
其中ht为h的转置矩阵，“·”表示矩阵乘法，本例中n＝1501；
g、用奇异值分解法，计算协方差矩阵γ的特征值矩阵λ和特征向量矩阵r，则存在公式
γ＝r·λ·rt(5)
其中λ为特征值由大到小排列的对角矩阵，λ＝diag[λ1,λ2,…,λn]，λ1,λ2,…,λn为特征值，r为各个特征值对应的特征向量矩阵，rt为r的转置矩阵，且满足rt·r＝r·rt＝e，其中e为单位矩阵；
h、h经线性映射，得到主成分矩阵φ，如公式
φ＝rt·h(6)
i、计算前p个特征值累计贡献率：
其中λj为特征值，p为所取特征值个数，1≤p≤n，j＝1,2,…,n；
j、若满足
y＝r·φ'(8)
则y为压制随机噪声后hankel重构矩阵，其具体形式可记为
本例中p＝300，保留φ的前300行主成分；
h、令x”(l)＝[y(1),y(2),···,y(n)]，则x”(l)即为对应x(l)的压制随机噪声后的压噪信号。本例中x”(l)＝[y(1),y(2),···,y(3001)]。