本发明涉及一种地震资料处理技术领域，特别是关于一种能够适应不同噪声环境的基于反比例模型的地震数据同时插值与去噪方法。
背景技术：
地震数据的插值重建与去噪是地震资料叠前处理流程中非常重要的一个环节。由于障碍物，风化带，经济成本等因素，野外采集的地震数据通常是不规则的。不规则地震数据会导致假频的产生，从而对后续处理流程比如srme多次波去除，3d偏移与成像等造成不良影响。
为了使地震数据规则化，我们通常应用一些数学算法来对缺失的地震道进行插值重建。基于稀疏变换的插值方法是应用较为广泛的一类算法，其中pocs算法以其简便的特点而备受欢迎。pocs算法起源于图像处理领域，于2006年被首次应用到地震数据重建领域(abmaandkabir,2006)。pocs算法的缺点之一就是抗噪性差。为了减小噪声对重建的影响，oropeza与sacchi(2011)引入了权重因子并且利用一个线性模型来求取权重因子。gaoetal(2013)进一步验证了加权pocs算法在重建含噪声数据时相对传统pocs算法的优势。葛子建等(2015)提出了利用数据驱动模型来计算权重因子，但是随着噪声等级的改变，该方法不能保证得到高信噪比的重建数据。
由于现有的几种加权方法不具备很好的噪声适应性，我们提出了一种基于反比例模型的计算加权因子的方法，将其应用到加权pocs算法中对缺失含噪声地震数据进行重建。通过控制加权因子随迭代次数的取值变化，我们的方法在不同的噪声环境下都取得了良好的重建效果。
技术实现要素：
针对上述问题，本发明的目的是提供一种基于反比例模型的地震数据同时插值与去噪方法，其能含噪声缺失地震数据进行重建，同时压制噪声，更好的适应不同的噪声等级，在环境噪声变化时，通过调整加权因子仍然能取得良好的重建效果。
为实现上述目的，本发明采取以下技术方案：一种基于反比例模型的地震数据同时插值与去噪方法，其包括以下步骤：1)输入初始不规则含噪声地震数据，对输入的数据进行正seislet变换，将输入数据变换到seislet域后得到了稀疏的表示；2)通过阈值算子保留大于阈值的seislet系数，将小于阈值的seislet系数置零；3)对阈值处理后的seislet系数进行反seislet变换得到时间域的数据；4)根据反比例模型确定加权因子；5)根据加权pocs算法的迭代公式和加权因子的取值，将步骤3)中得到的时间域数据根据确定的加权因子按照(i-αr)的比例保留，原观测数据按照加权因子α的权重的比例回加到反seislet变换后的数据中，得到本次迭代的重建结果；并判断k是否小于预先设定的最大迭代次数n，小于则返回步骤1)，再将重建结果作为下次迭代的输入，继续迭代，反之则输出最终重建结果，使缺失的数据被重建，噪声也得到压制；其中，i是单位矩阵；r是采样矩阵。
进一步，所述步骤1)中，如果是第1次迭代，则输入原始观测数据dobs。
进一步，所述步骤2)中，采用指数硬阈值模型来确定阈值的大小，阈值大小随迭代次数而变换，第k次迭代所选的阈值τk为：
τk＝τiec(k-1)/(n-1),c＝ln(τf/τi)，k＝1,2,...,n.
其中，n为预先设定的最大迭代次数；k为迭代次数；τi为初始阈值；τf为终止阈值。
进一步，所述步骤4)中，加权因子α的取值由如下反比例模型确定：
其中，k为迭代次数，q为控制因子，控制因子q用来决定加权因子的下降速率。
进一步，所述控制因子q根据地震数据的噪声强度来确定。
进一步，所述步骤5)中，加权pocs算法的迭代公式为：
式中，mn是第n次迭代后的重建结果；dobs是观测数据；s为seislet正变换；s-1为seislet反变换；
进一步，所述硬阈值算子
其中，xk为seislet系数。
本发明由于采取以上技术方案，其具有以下优点：1、本发明采用pocs算法作为重建算法，稀疏变换选择seislet变换，阈值模型采用指数阈值模型。为了在重建的同时压制噪声，引入了加权因子，并且通过反比例模型控制加权因子随迭代次数的下降速率，从而实现了能够适应环境噪声变化的地震数据同时插值与去噪。2、本发明根据反比例模型确定加权因子，可以根据环境噪声的大致强度来控制加权因子的取值变化，相对于现有的几种加权方法，能够更好的适应不同的噪声等级，在环境噪声变化时，通过调整加权因子仍然能取得良好的重建效果。3、本发明根据反比例模型来确定加权因子，弥补了现有的线性模型、数据驱动模型在噪声适应性方面的不足，能够在不同的噪声强度下达到较好的重建效果并保证计算效率。4、本发明为了压制噪声引入加权因子，加权因子的值由反比例模型来确定。
附图说明
图1是本发明的整体流程示意图；
图2是完整无噪声的模拟地震数据；
图3是50％不规则缺失并且添加随机噪声的模拟地震数据；
图4是对图3数据利用本发明方法迭代30次之后的重建结果；
图5是噪声等级较高时，各加权方法重建的信噪比恢复曲线对比图；
图6是噪声等级中等时，各加权方法重建的信噪比恢复曲线对比图；
图7是噪声等级较低时，各加权方法重建的信噪比恢复曲线对比图；
图8是完整不含噪声的实际地震数据；
图9是30％不规则缺失并添加随机噪声的实际地震数据；
图10是利用本发明方法重建之后的结果；
图11是不同加权方法对图9实际数据分别进行重建后的信噪比恢复曲线对比图。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
如图1所示，本发明提供一种基于反比例模型的地震数据同时插值与去噪方法，该方法采用凸集投影重建算法，在每次迭代中，首先通过seislet变换将输入信号在seislet域中稀疏的表示出来，通过阈值算子保留大于阈值的系数，再对保留的系数进行反seislet变换回到时间域，根据环境噪声的强度，确定一个合适的控制因子q的值，根据反比例模型，确定每次迭代的加权因子α的取值。再根据加权pocs算法的迭代公式，进行相应的数据保留和原观测数据的部分回加，得到该次迭代的重建结果，将其作为输入进行下一次迭代。经过足够的迭代，可以实现缺失地震数据的重建以及噪声的压制。
本发明具体包括以下步骤：
1)输入初始不规则含噪声地震数据，对输入的数据进行正seislet变换，将输入数据变换到seislet域后得到了稀疏的表示；其中，如果是第1次迭代，则输入原始观测数据dobs；
2)阈值处理：通过阈值算子保留大于阈值的seislet系数，将小于阈值的seislet系数置零；
在本实施例中，采用指数硬阈值模型来确定阈值的大小，阈值大小随迭代次数而变换，第k次迭代所选的阈值τk的计算公式为：
τk＝τiec(k-1)/(n-1),c＝ln(τf/τi)，k＝1,2,...,n.
其中，n为预先设定的最大迭代次数；k为迭代次数；τi为初始阈值，一般选最大的系数作为初始阈值；τf为终止阈值，通常选为一个接近或等于零的小值。随着迭代的进行，阈值大小从初始阈值按照指数规律衰减到终止阈值。
3)对阈值处理后的seislet系数进行反seislet变换得到时间域的数据
4)为了得到良好的重建结果，需要让加权因子的值随着迭代由大到小递减，则根据反比例模型确定加权因子α的取值；
加权因子α具体取值由如下反比例模型确定：
其中，k为迭代次数，q为控制因子，控制因子q用来决定加权因子的下降速率，一般从0和1之间选择q的值。
需要根据地震数据的噪声强度来确定q的取值：如果噪声强度较高，则需要选择一个较大的q值(比如0.8左右)，这样所求取的加权因子可以在前几次迭代中迅速下降；如果噪声强度中等，那么相应调整q的取值(大致在0.3-0.6之间选取)，使加权因子以稍微平缓的速率下降；如果噪声强度较弱，那么相应选择较小的q值(比如0.15左右)，使加权因子以更为平缓的速率下降。需要注意的是，该方法对环境噪声强度的判断只是定性的估计，因此对q值的选取可能需要若干次尝试和调整，q值选取范围并不局限于前述给出的选取范围。选定q值之后，每次迭代的加权因子α也就可以通过反比例模型所确定。
5)根据加权pocs算法的迭代公式和加权因子的取值，将步骤3)中得到的时间域数据
上述步骤5)中，加权pocs算法的迭代公式为：
式中，mn是第n次迭代后的重建结果；dobs是观测数据；i是单位矩阵；r是采样矩阵；s为seislet正变换；s-1为seislet反变换；
其中，硬阈值算子
其中，xk为seislet系数。
实施例：
下面通过数值模拟和实际数据的例子来进一步说明。如图2所示，是原始完整模拟地震道，对其进行50％随机缺失并添加随机噪音(如图3所示)，用基于反比例模型的加权pocs算法对其进行重建，迭代30次之后的重建结果如图4所示，可以看出缺失的地震道被重建，而且大部分噪声得到压制。继续比较几种不同加权模型在3种不同噪声环境下的重建结果。通过计算信噪比并绘制snr曲线来对比重建质量。从图5、图6及图7可以看出，通过调整控制因子q的值，本发明的方法在3种不同的噪声环境下都有很好的重建效果。可以观察到线性模型虽然最终都能得到高信噪比的重建结果，但是在噪声强度较高的环境下，其收敛速率比反比例模型要慢很多，数据驱动模型在噪声强度较高的情况下重建质量良好，但是在噪声强度较低时，重建结果却不理想。最后，通过实际数据来进行验证。如图8所示，是完整的无噪声实际地震数据，移除30％的地震道并添加随机噪声(如图9所示)，分别用四种加权方法对其进行重建；如图10所示，是本发明提出的方法的重建结果，可以看出重建效果良好。从图11中的信噪比恢复对比曲线可以看出，基于反比例模型的加权方法相比其他方法，能够达到更高的信噪比。
上述各实施例仅用于说明本发明，各个步骤都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别步骤进行的改进和等同变换，均不应排除在本发明的保护范围之外。