本发明属于地震勘探数据处理领域，涉及一种地震数据插值方法，特别是一种基于改进的变分正则化的地震数据规则化方法。
背景技术：
地震勘探是寻找石油和天然气能源的主要方法。地震勘探包括数据采集，数据处理和解释三部分。地震采集的数据直接影响着处理和解释的结果。由于经济条件的限制，地形的影响和坏道等因素，采集的数据在空间上往往不满足采样定理，而这样的数据会严重影响后续的处理质量，比如多次波消除和偏移成像。因此地震数据规则化是地震数据处理的重要部分。
在地震处理流程之初，可以通过在规则网格上恢复地震数据来降低空间假频和非规则采样带来的影响。恢复的方法有很多种，大致可以分为三类。波场算子类方法，该方法基于kirchhoff积分算子，通过一个连续算子的积分得到延拓的波场；滤波类方法，该方法通过数据和预测滤波算子做卷积重建地震数据；变换类方法，非规则采样的数据可以使用几种变换(fourier,radon,wavelet等)通过两步方法进行规则化。第一步估计出变换系数，当采样不规则时，直接的正变换给出的系数是失真的，通过反演可以给出一个更好的估计；第二步通过逆变换给出规则网格上的数据。
滤波类算法通过数据和预测滤波算子做卷积重建地震数据。spitz(1991)提出了一种t-x域地震道空间插值方法,它可以有效对付空间假频现象。在f-k域、τ-p域及t-x等内插方法中，t-x域的效果是最理想的。spitz揭示了这样一个事实：等间距地震剖面中的线性同相轴是可以通过地震数据本身来精确估计的，而与原空间采样间隔无关。spitz算法的数学基础是单位步长预测滤波。为了求得滤波因子，需要求解两组复杂的复线性系统方程，因而运算工作量很大，这在一定程度上影响了该方法的吸引力。porsani(1999)发现，这一问题可以通过针对利用偶数道数据分量来预测奇数道数据分量而设计的半步长预测滤波器来克服。新算法只需求解一组简单且只与地震道数据本身有关的线性系统方程,因而极大地提高了处理效率,并且使程序实现更容易。国九英等(1996)在spitz的基础上对该方法进行了改进，但该方法不但要两次用到最小平方算法,而且运算速度较慢。soubaras(1997)提出了f-x投影滤波插值，该方法可以增加规则假频数据的采样密度。gulunay等(1996，1997)提出f-k域抗假频的插值方法，分别用原始道集中奇数道所组成的f-k变换和由原始道集中偶数道组成的道集的f-k变换，进而再计算一个滤波因子，用于对原始道集的f-k变换进行滤波处理，从而可得内插道的f-k变换，用以完成含有空间假频数据的道内插。
技术实现要素：
1、针对上述现有技术，本发明所要解决的技术问题是，spitz方法是一种处理等间距道上地震数据空间假频的常用地震数据规则化(道插值)方法。由于一般的spitz插值方法在求解预测系数时需要求解线性方程组，会造成累积的误差比较大，插值后带来的噪声较大。因此我们将一种改进的变分正则化方法引入线性方程组的计算中，提出了一种新的插值方法，使得地震数据插值达到更好的效果。
2、发明采用的技术方案
为解决上述技术问题，本发明是按如下方式实现的：该方法包括以下步骤
步骤一，采集地震数据，做数据预处理将原始数据转化为2d规则道集记录。
假设u(t,x)为二维均匀采样地震数据，t和x分别为地震数据的时间和空间变量，数据含n地震道，有l个同相轴(l<n)。假设每个同相轴的子波为wk(t)，指第k个同相轴的子波函数，子波函数可从数据中提取。那么地震数据可表示为
步骤二，对地震数据每道做时间域fft，转化为频率-空间域数据。
公式(1)的fourier变换为
其中，bk为第k条同相轴在相邻道的时间延迟。令zk＝exp(i2πfbk)，对于某一个频率切片，上式可写为
u(f)＝w(f)z(f).(3)
步骤三，对于一个频率切片，基于改进的变分正则化方法估计预测系数。
将某个频率切片看作空间方向上的平稳随机过程，对该随机过程使用前向一步线性预测器可得如下结果
后向一步线性预测器可得
其中，pj(f)为预测系数，
假设插值后的数据为u′(f,x)，p′(f)为插值后的预测系数，b′k是插值后第k条同相轴相邻道之间的时间延迟，
由于
对比前面的方程式可得
这样，由推导出的这个性质可得，插值前低频的数据估计所得的预测系数可以用来预测插值后的数据。
由上所述，一步前向预测系数为(p1(f),p2(f),k,pl(f))，求解如下方程组可以求得这个预测系数：
其中，r(f)是地震道数据频域切片的自相关矩阵。
利用改进的变分正则化方法求解该线性方程组，其步骤如下：
步骤三一，将线性方程组(8)写作一般形式
kx＝yδ(9)
步骤三二，为了克服系统不适定性带来的问题，将问题变为求解极小化泛函
minmα[x,yδ]＝||kx-yδ||2+α||lx||2(10)
使得问题变为适定问题。这里α＞0被称作正则参数，l称为尺度算子。
步骤三三，记φ＝k*k，在本发明中，取
(φ2+αi)x＝φzδ(11)
这里zδ＝k*yδ。
步骤三四，使用预条件共轭梯度法求解方程(11)
步骤四，利用求出的预测系数，对该频率切片做频率域插值。
步骤五，遍历所有的频率切片，重复步骤一到四
步骤六，做时间域逆fft，将数据转到时域。
至此，完成了地震数据规则化，从而去除了地震数据中的假频成分。
3、发明起到的积极效果
在本发明中我们将改进的变分正则化方法引入线性方程组的计算中，提出了一种新的插值方法。由于改进的变分正则化方法有比tikhonov正则化方法更强的稳定性，因此带来了地震数据规则化(道插值)的更好的效果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明基于改进的变分正则化的地震数据规则化的流程图。
图2是实例说明中带有假频的原始地震道数据。
图3是实例说明中传统spitz方法地震数据规则化(道插值)后的结果。
图4是实例说明中改进的变分正则化方法数据规则化(道插值)后的结果
具体实施方式
下面结合附图对本发明做进一步详细描述
如图1所示，本发明实施的步骤包括
s101：采集地震数据，做数据预处理将原始数据转化为2d规则道集记录。
假设u(t,x)为二维均匀采样地震数据，数据含n地震道，有l个同相轴(l<n)。假设每个同相轴的子波为wk(t)，那么地震数据可表示为
s102：对地震数据每道做时间域fft，转化为频率-空间域数据。
公式(1)的fourier变换为
其中，bk为第k条同相轴在相邻道的时间延迟。令zk＝exp(i2πfbk)，对于某一个频率切片，上式可写为
u(f)＝w(f)z(f).(3)
s103：对于一个频率切片，基于改进的变分正则化方法估计预测系数。
将某个频率切片看作空间方向上的平稳随机过程，对该随机过程使用前向一步线性预测器可得如下结果
后向一步线性预测器可得
其中，pj(f)为预测系数，上式中利用了前向和后向预测系数之间的关系。
假设插值后的数据为u′(f,x)，p′(f)为插值后的预测系数，b′k是插值后第k条同相轴相邻道之间的时间延迟，
由于
对比前面的方程式可得
这样，由推导出的这个性质可得，插值前低频的数据估计所得的预测系数可以用来预测插值后的数据。
由上所述，一步前向预测系数为(p1(f),p2(f),k,pl(f))，求解如下方程组可以求得这个预测系数：
其中，r(f)是地震道数据频域切片的自相关。
具体的，利用改进的变分正则化方法求解该线性方程组，步骤如下：
第一步，将线性方程组写作一般形式
kx＝yδ(9)
第二步，为了克服系统不适定性带来的问题，将问题变为求解极小化泛函
minmα[x,yδ]＝||kx-yδ||2+α||lx||2(10)
使得问题变为适定问题。这里α＞0被称作正则参数，l称为尺度算子。
第三步，记φ＝k*k，在本发明中，取
(φ2+αi)x＝φzδ(11)
这里zδ＝k*yδ。
第四步，使用预条件共轭梯度法求解方程(11)
s104：利用求出的预测系数，对该频率切片做频率域插值。
s105：遍历所有的频率切片，重复s101到s104。
s106：做时间域逆fft，将数据转到时域。
至此，完成了地震数据规则化，从而去除了地震数据中的假频成分。
实现本发明基于改进的变分正则化的地震数据规则化方法的一个具体实例如下：
如图2所示，是一个简单的三同相轴地震数据，共有16个地震道，检波器间的距离为10m，时间采样率为2ms，地震子波为主频为30hz的雷克子波。图3是传统的spitz方法得到的结果，图4是改进的变分正则化方法得到的结果。由图中可以看出，插值后数据虽然较准确地内插了同相轴数据，但是带来了较多的噪声污染。采用改进的变分正则化方法保证了数据插值的稳定性，得到了更好的结果。