专利名称：基于可靠高阶统计量样本的地震子波计算机估计方法
技术领域：
本发明属于地震信号处理领域。
背景技术：
地震子波估计是地震信号处理中的关键问题之一，是其他一些地震处理技术(如地震正演，地震反演，反褶积等)的基础。利用地震信号实现快速准确的地震子波估计具有重要的应用价值。
高阶统计量由于其能解决信号处理中的相位估计，且对加性高斯噪声不敏感，从而在许多领域中得到广泛的应用。人们提出了许多基于高阶统计量的信号处理算法，来替代传统的基于二阶统计量的信号处理方法，在许多领域中取得了很好的应用效果。在地震信号处理中，人们提出了基于高阶统计量的地震地震子波估计算法，相对于基于二阶统计量的地震子波估计算法，可以解决地震子波相位估计的问题(Lu，W.，2005a，Non-minimum-phase waveletestimation using second and higher order momentGeophysical prospecting，53，149-158.Velis，D.R.，and T.J.Ulrych，1996，Simulated annealing wavelet estimation via fourth-order cumulantmatchingGeophysics，61，1939-1948.Lazear，G.L.，1993，Mixed-phase wavelet estimation usingfourth-order cumulantsGeophysics，58，1042-1051.)。
利用高阶统计量进行地震子波估计，一般先要利用地震信号估计出地震子波的高阶统计量，然后利用地震子波的高阶统计量进行地震子波的估计。地震子波估计的精度受地震子波的高阶统计量估计精度的影响很大。相对于二阶统计量，高阶统计量的估计存在更大的方差，为了得到准确的高阶统计量的估计，往往需要较长的地震信号(Lu，W.，2005a，Non-minimum-phase wavelet estimation using second and higher order momentGeophysicalprospecting，53，149-158.Velis，D.R.，and T.J.Ulrych，1996，Simulated annealing waveletestimation via fourth-order cumulant matchingGeophysics，61，1939-1948.Lazear，G.L.，1993，Mixed-phase wavelet estimation using fourth-order cumulantsGeophysics，58，1042-1051.)。
而考虑到实际地震资料中，地震子波是时变和空变的，所以不能用太长的地震信号进行子波估计。人们已经发现由观测信号估计得到的信道(在地震信号处理中，信道被称为地震子波)的高阶统计量样本的可靠性是不一样的，提出了利用可靠性高的高阶统计量切片来进行信道估计的方法并取得了好的估计结果(Lu，W.，2005b，Blind channel estimation usingzero-lag slice of third order momentIEEE signal processing letters，12，725-727.Petropulu，A.P.，and Pozidis，H.，1998，Phase reconstruction from bispectrum slicesIEEE Trans.Signal Process.，46，527-530.)。
地震信号x＝[x1，...，xi，...xL]T(其中下标表示信号的样本点号)可以表示为反射系数s＝[S1，...，Si，...SN]T和地震子波b＝[b1，...，bi，...，bQ]T的卷积，并且带有加性白噪声n＝[n1，...，ni，...，nL]T，即x＝b*s+n (1)根据公式(1)，假设反射系数是独立同分布的序列，则可以用地震信号x的高阶累计量来估计地震子波b的高阶矩(Lu，W.，2005a，Non-minimum-phase wavelet estimation usingsecond and higher order momentGeophysical prospecting，53，149-158.Velis，D.R.，and T.J.Ulrych，1996，Simulated annealing wavelet estimation via fourth-order cumulant matchingGeophysics，61，1939-1948.Lazear，G.L.，1993，Mixed-phase wavelet estimation usingfourth-order cumulantsGeophysics，58，1042-1051.)，这就是著名的BBR公式mrb=crxγrs,...(2)]]>其中mrb是信道b的r阶矩，crx是接收信号x的r阶累计量，γrs是和s有关的常数(张贤达，时间序列分析——高阶统计量方法，清华大学出版社，1996)。
我们常使用的零均值的序列x(i)的三阶累计量和四阶累计量定义如下c3x(k1,k2)=Σi=1Lxi+k1xi+k2xi...(3)]]>c4x(k1,k2,k3)=1LΣi=1Lxi+kixi+k2xi+k3xi-(1LΣi=1Lxi+k1-k2xi)(1LΣi=1Lxi+k2-k3xi)]]>-(1LΣi=1Lxi+k2-k3xi)(1LΣi=1Lxi+k1-k3xi)...(4)]]>-(1LΣi=1Lxi+k1-k3xi)(1LΣi=1Lxi+k3-k2xi)]]>其中k1，k2，k3表示时延，其取值范围是[-L+1，L-1]之间的整数。
利用公式(2)我们可以直接估计地震信号的高阶累积量来得到地震子波的高阶矩，因此在后续的方法描述中，地震子波b的r阶矩mrb是假设已知的。
发明内容
本发明的目的在于提出一种具有高可靠性的基于高可靠性高阶统计量样本的地震子波的计算机估计方法。
本发明的特征之一在于，该方法依次含有以下步骤步骤(1)向计算机输入以下设定量，并初始化迭代计数器地震观测信号x，记忆因子α∈(0，1)，误差停止条件ε，最大迭代次数T，迭代次数计数器t＝0，相邻两次迭代误差记录e＝2ε；步骤(2)向计算机输入三阶统计量样本，并从中提取靠近零时延点m3b(0，0)的样本点，用m表示，具体表示如下m3b(-k，0)，...，m3b(0，0)，...，m3b(k，0)其中2k+1≥Q，以及m3b(j，l)，j+l≤C，j≥l，其中C是远小于Q的正整数；其中b为待求的地震子波组成的列向量；步骤(3)计算机按以下步骤计算初始地震子波估计b0步骤(3.1)计算机按下式计算地震观测信号x的自相关m2x(j)，m2x(j)=Σi=1Lxixi+j,]]>-Q＜j＜Q，其中L为地震观测信号x的长度，Q为估计地震子波的长度，以下用m2x表示m2x(j)；步骤(3.2)按下式计算地震子波的自相关m2bm2b=m2xw,]]>w为一个长度为2Q+1的窗函数；步骤(3.3)按下式计算地震子波的能量谱估计Pb(ω)，Pb(ω)=|FFT{m2b}|,]]>FFT为快速傅立叶变换；步骤(3.4)按下式计算初始地震子波估计b0，b0＝IFFT{[Pb(ω)]1/2}；步骤(4)构造待求地震子波列向量b的相关矩阵Afull(b)，并从中选取与步骤(2)所述的靠近m3b(0，0)的样本点所对应的行，得到一个估计地震子波列向量b的非线性方程组A(b)b＝m；
其中 步骤(5)引入一个中间变量ba，ba为已经得到的地震子波的一个估计，初始化时为b0，建立如下的优化函数minα||A(ba)ba-A(ba)b||22+(1-α)||A(ba)b-m||22]]>其中‖·‖2为向量的2范数，从而得到一个新的地震子波估计b；步骤(6)按以下步骤更新步骤(5)得到的中间变量ba步骤(6.1)由优化函数求解得局部最小值表达式如下b＝αb+(1-α)ba；但是由于估计过程中幅度不确定的问题，必须先对ba和b归一化ba′＝ba/‖ba‖2，b′＝b/‖b‖2；步骤(6.2)ba＝αba′+(1-α)b′步骤(7)计算e，并判断e＝max(|b-ba|)＞ε，t＝t+1，并且t＜T；步骤(8)若e＞ε，则返回步骤(5)；否则输出估计的地震子波向量b，b＝ba。
本发明的特征之二在于，高阶统计量样本可以是四阶统计量样本，相应地
mfull=m4b(-Q+1,0,0)······m4b(0,0,0)······m4b(Q-1,0,0)m4b(1,2,3)···m4b(i,j,l)···,]]>选择靠近零时延点的样本构造的A(b)和m分别为 m=m4b(-k,0,0)······m4b0(0,0,0)······m4b(k,0,0)m4b(1,2,3)···m4b(i,j,l)···,]]>其中2k+1≥Q，i＜j＜l＜Q并且i，j，l都是比较小的整数。
实验仿真验证在我们的仿真实验中，使用HP-Pavilion N5425计算机。利用不同长度的i.i.d(独立同分布)信号s按x＝b*s+n生成在10dB高斯白噪声污染下的观测数据，对两个不同的地震子波(记为W＝1，2)，地震子波长度为32，在不同的实验设置下进行100次独立实验。我们使用归一化相似度来衡量估计信道ba与真实信道b的相似度，定义如下 表1，表二给出了本发明提出的方法在时域和频域实现所得到的统计结果。其中表一所选择的高阶统计量样本为三阶阶统计量零时延切片，以及非零时延样本点m3b(1，2)，m3b(1，3)，m3b(2，3)，m3b(1，4)，m3b(2，4)，m3b(3，4)，m3b(1，5)；表二所选择的高阶统计量样本为四阶统计量零时延切片。为了比较，我们同时给出一种基于累计量匹配(Lazear，G.L.，1993，Mixed-phase wavelet estimation using fourth-order cumulantsGeophysics，58，1042-1051.)的地震子波估计结果。
图2，图3给出本发明提出的方法在使用四阶统计量时域实现所得到的统计结果的图示。图中，用-o-标识估计子波的均值，-x-标识真实子波，阴影标识地震子波估计方差。N为反射系数序列长度，Q为估计地震子波长度。
可以看出，本发明提出的方法可以有效地提高地震子波的估计精度，特别是在地震数据短的情况。同时，本发明提出的方法在地震子波估计长度大于实际长度时，得到的结果变化不大。
图1、本发明所属方法的计算机程序流程图；图2、地震子波1的统计结果曲线图a，N＝1024，Q＝32 b，N＝1024，Q＝42c，N＝128，Q＝32 d，N＝128，Q＝42图3、地震子波2的统计结果曲线图a，N＝1024，Q＝32 b，N＝1024，Q＝42c，N＝128，Q＝32 d，N＝128，Q＝42其中，Q为估计地震子波长度；N为反射系数序列长度；-o-为估计地震子波的均值；-x-为真实地震子波子波；阴影为地震子波估计方差。
具体实施例方式我们提出并实现了用少数可靠的地震子波的高阶统计量样本点估计地震子波的算法。这一算法的特点是可以自由选择估计得到的地震子波的高阶统计量的样本点。我们的算法把每个样本点独立看待，按个数选择。现有的技术一般不考虑地震子波的高阶统计量的样本的估计可靠性，即使考虑地震子波的高阶统计量样本估计的可靠性(P.Petropulu and H.Pozidis，1998，Phasereconstruction from bispectrum slices，IEEE Trans.Signal Processing，46，527-530.Lu，W.，2005b，Blind channel estimation using zero-lag slice of third order momentIEEE signal processingletters，12，725-727.)，也是按切片考虑的。
·建立新的优化目标函数，实现相应的迭代优化算法来求解非线性方程组，在每一步迭代中，求解一个线性方程组。迭代算法所使用的初始地震子波可以采用利用二阶统计量估计得到的零相位子波(Lu，W.，2005b，Blind channel estimation using zero-lag slice of thirdorder momentIEEE signal processing letters，12，725-727.)，且收敛速度快。
·地震子波估计算法可以在时域实现，在利用整个高阶统计量切片时也可以在频域实现，从而可以利用快速傅立叶变化更进一步减少算法的计算量。
·可以有效地选择可靠性高的高阶统计量的样本点，使得所提出的方法突破了需要较长数据获得地震子波估计的限制，在很短的观测数据的情况下得到好的地震子波估计。
·本发明所提出的技术也可以应用于信道估计，系统辨识领域。
本发明按以下步骤实施1.输入观测信号以及由观测信号估计得的高阶统计量，地震子波的阶数(估计值)，地震子波初始估计，最大迭代次数T，迭代停止误差条件ε，记忆因子α∈(0，1)。
2.选择需要的高阶统计量样本。
3.使用高阶统计量样本和地震子波的估计构造线性方程组，求解该线性方程组(时域或频域方法)得到一个新的地震子波的估计。
4.更新地震子波的估计。
5.如果停止条件不满足，回到3。
6.输出地震子波的估计。
具体原理如下1.高阶统计量样本选择和线性方程组的构造地震子波r阶矩的定义为mrb(k1,···,kr-1)=1QΣi=1Q[bi+k1bi+k2···bi+kr-1]bi,-Q<kj<Q,1<j<r-1...(5)]]>式中，kj同样表示时延，Q是地震子波的长度。
重写公式(2)为矩阵形式如下Afull(b)b＝mfull(6)式中，b为待求的地震子波组成的列向量，Afull(b)为一个与b有关的矩阵，mfull为由高阶统计量样本mrb(k1，...，kr-1)组成的列向量，
下面以三阶统计量为例，说明如何选择构造Afull(b)和mfull。利用三阶统计量，则可以得到 mfull=m3b(-Q+1,0)······m3b(0,0)······m3b(Q-1,0)m3b(1,2)···m3b(j,k)···...(7)]]>式中，Afull(b)和mfull的行是一一对应的。
一般来说，靠近零时延点m3b(0，0)的样本点是可以比较可靠地利用地震信号估计得到，因此，只保留mfull中接近零时延点的样本点构成m，同时，也只保留Afull(b)对应的行来构成A(b)，得到估计地震子波的方程组如下A(b)b＝m (8)其中 m=m3b(-k,0)···m3b(0,0)···m3b(k,0)···m3b(j,l)···,]]>k，j，l都是比较小的整数。
为了在后面的计算中增强数值稳定性和保证每一步迭代结果的唯一性，需要保证A(b)的秩为Q。实现这一目标的最简单方法是m中至少包含样本m3b(-k，0)，...，m3b(0，0)，...，m3b(k，0)，其中2k+1≥Q。
需要指出的是，当所选用的高阶统计量样本是由1个完整的1维高阶统计量切片组成，则公式(8)可以表示为地震子波和一个由地震子波导出的向量的互相关(Lu，W.，2005b，Blindchannel estimation using zero-lag slice of third order momentIEEE signal processing letters，12，725-727.)，所以可以在频域利用快速傅立叶变换求解公式(8)表示的方程组。
利用其他阶的高阶统计量，或者不同阶的高阶统计量组合，可以按照相似的方法构造矩阵A(b)和向量m。
如四阶统计量对应的Afull(b)和mfull的形式为 mfull=m4b(-Q+1,0,0)······m4b(0,0,0)······m4b(Q-1,0,0)m4b(1,2,3)···m4b(j,k,l)···...(9)]]>选择靠近零时延点的样本构造的A(b)和m分别为
m=m4b(-k,0,0)······m4b(0,0,0)······m4b(k,0,0)m4b(1,2,3)···m4b(i,j,l)···,]]>其中2k+1≥Q，i＜j＜l＜Q并且i，j，l都是比较小的整数。
2.更新算法。
公式(8)是一个非线性方程组，直接求解是不可行的，因此引入辅助变量ba，建立如下优化函数minα||A(b)ba-A(b)b||22+(1-α)||A(b)ba-m||22...(10)]]>式中，α为预先设定的记忆因子，b为已得到的地震子波的一个估计，ba为待求的地震子波的新的估计，‖·‖2为向量的2范数，其定义是||b||2=Σi=1Qbi2.]]>利用公式(10)给出的目标函数，可以将公式(8)所给出的非线性问题转化为一个迭代的求解线性方程组的问题，具体算法流程见图1。其中初始地震子波b可以用现有的方法进行估计，如双谱切片估计(P.Petropulu and H.Pozidis，1998，Phase reconstruction from bispectrumslices，IEEE Trans.Signal Processing，46，527-530.)，累计量切片线性组合法(J.A.R.Fonollosaand V.Vidal，“System identification using a linear combination of cumulant slices，”IEEE Trans.Signal Processing，vol.41，no.7，pp.405-2412，July 1993.)等。在我们的算法实现中，我们是用信号的能量谱估计初值，具体做法是1.估计地震信号的自相关m2x(j)=Σi=1Lxixi+j,-Q<j<Q;]]>2.估计地震子波的自相关m2b=m2xw,]]>式中m2b为地震子波的自相关，m2x为地震信号的自相关，w为一个长度为2Q+1的窗函数。
3.得到地震子波的能量谱估计Pb(ω)=|FFT{m2b}|,]]>其中FFT{}是快速傅立叶变换；4.得到初始地震子波的估计b＝IFFT{[Pb(ω)]1/2}，其中IFFT{}是快速傅立叶逆变换；
相对于其他的方法(P.Petropulu and H.Pozidis，1998，Phase reconstruction from bispectrumslices，IEEE Trans.Signal Processing，46，527-530.J.A.R.Fonollosa and V.Vidal，“Systemidentification using a linear combination of cumulant slices，”IEEE Trans.Signal Processing，vol.41，no.7，pp.405-2412，July 1993.)，上述方法可以有效地减少估计初始地震子波的计算量。
表1三阶统计量试验统计结果。真实地震子波长度为33，100次独立试验。
表2四阶统计量试验统计结果。真实地震子波长度为33，100次独立试验。
权利要求
1.基于可靠高阶统计量样本的地震子波计算机估计方法，其特征在于，该方法依次含有以下步骤步骤(1)向计算机输入以下设定量，并初始化迭代计数器地震观测信号x，记忆因子α∈(0，1)，误差停止条件ε，最大迭代次数T，迭代次数计数器t＝0，相邻两次迭代误差记录e＝2ε；步骤(2)向计算机输入三阶统计量样本，并从中提取靠近零时延点m3b(0，0)的样本点，用m表示，具体表示如下m3b(-k，0)，...，m3b(0，0)，...，m3b(k，0)，其中2k+1≥Q，以及m3b(j，l)，j+l≤C，j≥l，其中C是远小于Q的非负整数，其中b为待求的地震子波组成的列向量；步骤(3)计算机按以下步骤计算初始地震子波估计b0步骤(3.1)计算机按下式计算地震观测信号x的自相关m2x(j)，m2x(j)=Σi=1Lxixi+j-Q<j<Q,]]>其中L为地震观测信号x的长度，Q为估计地震子波的长度，以下用m2x表示m2x(j)；步骤(3.2)按下式计算地震子波的自相关m2bm2b=m2xw,]]>w为一个长度为2Q+1的窗函数；步骤(3.3)按下式计算地震子波的能量谱估计Pb(ω)，Pb(ω)=|FFT{m2b}|,]]>FFT为快速傅立叶变换；步骤(3.4)按下式计算初始地震子波估计b0，b0＝IFFT{[Pb(ω)]1/2}；步骤(4)构造待求地震子波列向量b的相关矩阵Afull(b)，并从中选取与步骤(2)所述的靠近m3b(0，0)的样本点所对应的行，得到一个估计地震子波列向量b的非线性方程组A(b)b＝m；其中 步骤(5)引入一个中间变量ba，ba为已经得到的地震子波的一个估计，初始化时为b0，建立如下的优化函数minα||A(ba)ba-A(ba)b||22+(1-α)||A(ba)b-m||22]]>其中||·||2为向量的2范数，从而得到一个新的地震子波估计b；步骤(6)按以下步骤更新步骤(5)得到的中间变量ba步骤(6.1)由优化函数求解得局部最小值表达式如下b＝αb+(1-α)ba；但是由于估计过程中幅度不确定的问题，必须先对ba和b归一化ba′＝ba/||ba||2，b′＝b/||b||2；步骤(6.2)ba＝αbba′+(1-α)b′步骤(7)计算e，并判断e＝max(|b-ba|)＞ε，t＝t+1，并且t＜T；步骤(8)若e＞ε，则返回步骤(5)；否则输出估计的地震子波向量b，b＝ba。
2.根据权利要求1所述的基于可靠高阶统计量样本的地震子波计算机估计方法，其特征在于高阶统计量样本是四阶统计量样本，相应地 mfull=m4b(-Q+1,0,0)······m4b(0,0,0)······m4b(Q-1,0,0)m4b(1,2,3)···m4b(i,j,l)···]]>选择靠近零时延点的样本构造的A(b)和m分别为 m=m4b(-k,0,0)······m4b(0,0,0)······m4b(k,0,0)m4b(1,2,3)···m4b(i,j,l)···]]>，其中2k+1≥Q，i＜j＜l＜Q并且i，j，l都是比较小的整数。
全文摘要
本发明属于地震观测信号处理技术领域，其特征在于利用可靠性的三阶或四阶统计量在零时延点附近的样本点来进行地震子波的估计，同时还提出了一种计算量较少的能量谱估计法来估计初始地震子波估计，在结合根据预先设定的记忆因子而建立的优化函数，使得计算地震子波r阶非线性方程组线性化，通过迭代解法，便可求出待求的地震子波新向量。本发明提出的方法有效地提高了地震子波的估计精度，特别是在地震数据短的情况下更为明显。
文档编号G01V1/28GK1877363SQ20061008962
公开日2006年12月13日 申请日期2006年7月7日 优先权日2006年7月7日
发明者陆文凯, 张映松 申请人:清华大学