地球物理资料数字处理(第二十二十一讲)

第四节最小熵反褶积(MinimumEntropyDeconvolution)地球物理技术具有高科技密集，更新速度快的特点。每一次新技术的出现及其应用都对数字处理技术产生巨大的冲击。尤其是随着现代信息处理技术的发展和计算机技术的进步，更加促进了新技术在地震资料数字处理中的应用。最小熵反褶积就是将信息论中熵的概念引入到反褶积中的技术。与最小平方反褶积不同的是，它不需要假定反射系数是白噪，子波为最小相位的条件，从这个意义上讲，它比最小平方反褶积方法先进，从另外的角度看，它的算法较为复杂，这又是这种方法的缺点。这里讨论最小熵反褶积技术，意义不在于其方法本身。目的是把反褶积技术引入更广阔的领域。最小熵反褶积的目的是压缩地震子波，提高地震记录的分辨率。本节要点：熵的基本概念；最小熵反褶积原理；最小熵反褶积原理的步骤。一、熵的基本概念(TheBasicConceptofEntropy)研究最小熵反褶积，必须首先理解熵的概念。熵，最早并不是地震资料数字处理技术中的专有名词。它是由物理学中热学的概念借用而来。在物理热力学中，研究分子运动状态，用熵来衡量分子热运动的秩序。从分子运动的观点看，分子热运动要从有序趋于混乱，这时分子的熵发生变化，熵越大，分子运动越无秩序，规律性不强，反之熵越小，分子运动越有秩序。可见，最早熵是用来衡量分子运动秩序的物理量。在信息论中，熵用来衡量系统内信息的不确定性或不可预测性，熵越大系统内信息的不确定性就越大，反之，则不确定性就越小二、最小熵反褶积原理（ThePrincipleOFMinimumEntropyDeconvolution）根据地震褶积模型，地震记录()()()xtbtt，最小熵反褶积问题归结为1．寻找一线性（LinearDeconvolutionFactor）反褶积算子()ft2．使得地震记录()xt与反褶积算子()ft褶积后其输出为(WaveShape)波形简单，(Polarity)和延迟时间(TimeDelay)未知的尖脉冲序列。即输出为反射系数序列()t，则有下式成立()()()tftxt（5-4-1）从熵的角度分析，最小熵反褶积的输出准则是希望得到褶积输出t是由一系列极性和延迟时间未知（由地下介质的性质决定）的尖脉冲序列组成，且应使信号的秩序最好，规律性强，亦即熵最小。因此，最小熵反褶积的关键是如何使褶积的输出达到最小。为了提高反褶积算子的精度，实际确定反褶积算子()ft时，不是使用一道地震记录，而是选区一组地震道组成一个道集，如图（5-11）。这时，反褶积的输出亦为一个道集(1)isyiN，则,0mijkijKKyfx（5-4-2）其中，求取反褶积因子的道集,ijx，i为道序号(1,2,,)iN,j为记录样点采样序号((1,2)jT,m为反褶积因子()ft的长度,下面的问题是如何衡量输出道集,ijy的熵达到最小。为了确定输出波形达到熵最小，构造函数V，成为最大方差范数图5-11确定最小熵反褶积因子输入道集4121TijjiTijjyVy（5-4-3）1NiiVV（5-4-4）可以证明24210001TTijijjjyyV分析最大方差范数的特点，若取输出为尖脉冲时，即输出波形最为简单，亦即熵最小的情况，则1V，若输出为幅度相等的2个尖脉冲时，V=0.50。若输出为10个幅度相等的尖脉冲时，V=0.10。虽然，当脉冲幅度相等时，脉冲的个数越多，则最大方差范数V的值越小，且最大方差范数不受脉冲的极性和出现时间影响。对于输出波形为幅度不等的尖脉冲时，可以参看图（5-12）了解最大方差范数的特点。可以理解，最大方差范数V可以作为衡量输出信号熵大小的准则。下面讨论如何确定反褶积因子()ft。根据上面建立的最小熵反褶积准则及最大方差范数的特点，所求的反褶积因子()ft必满足使最大方差范数达到最大，也就是输出地震记录的熵最小，所以()ft必须满足下式10NiiKKVVff（5-4-5）图5-12简单脉冲序列最大方差值数V的特征所以412121TijNjTiKKijjyVffy(5-4-6)为简化求导，令21Tiijjuy，则（5-4-6）可写为123111440NTTijijiiijiijijjKKkyyVVuyVuyfff（5-4-7）由（5-4-2）式知,ijijKKyxf则（5-4-7）又可写成下式123,,,11011NTmNTiiijijKiijijKijijVufxxuyx对上式整理，又可写成下式123,,,01111mNTNTniiijijKiijijKijijfVuxxuyx0,10,1Kmm令1,,11NTxxiiijijKijRKVuxx323,11()NTxyiijijKijRKuyx将3()()xxxyRKRK和代入（5-4-7）得30()()mxxxyfRKRK（5-4-8）写成矩阵形式30313(0)(1)()(0)(1)(0)(1)(1)()(1)(0)()xxxxxxxyxxxxxxxyxxxxxxmxyRRRmfRRRRmfRRmRmRfRm（5-4-9）分析（5-4-7）式知，该方程组并不是一个线性方程组，而是一个复杂的非线性高次方程，将（5-4-7）式写（5-4-8）式制式起到使其表达式简洁的作用，因此不能用求解线性方程组的方法求解。可以采用迭代法求解。迭代求解过程为1）给定一初始线性褶积算子()ft2）根据,0mijKijKKyfx计算,ijy值。3）根据1,,11()NTxxiiijijKijRKVuxx323,11()NTxyiijijKjjRKuyx计算()xxRK和3()xyRK的值4）将计算所得的()xxRK和3()xyRK的值代入矩阵（5-4-9）并求解，可得一新的褶积算子()ft5)重复(2)~(4)步骤,直至得到满意的输出ijy为止,这时得到的()ft可做为反褶积算子。上述的迭代过程一般4~6次就可以结束。由于所解方程为非线性，这时输出的反褶积算子并不能保证最大方差范数V是唯一的极大值，但此时的反褶积算子()ft仍不失其使用意义。在迭代的过程中为确保最小熵准则的一致性，必须对每次迭代所得的反褶积算子进行归一化处理，以使每次计算的最大方差范数V是使用统一的数量标准比较。在建立最小熵反褶积准则时，我们仅考虑输出具有简单外形的条件，并没有考虑输出脉冲的极性和其时间延迟，这可能导致输出脉冲的极性与实际反射系数的极性相反，或与实际的时间延迟不一致的情况出现，这是最小熵反褶积方法的一个缺点，只要我们合理选择初始的反褶积算子()ft，我们可以避开这个问题。一般，若选线性反褶积算子的初始值为一尖脉冲时，其间脉冲出现的位置应接近实际位置。图（5-13）是利用人工合成记录进行最小熵反褶积的算例。图（5-13）（a）是原始模型，（a）(b)(c)图（5-13）最小熵反褶积算例上部是输入子波，下部是反射系数序列。图(5-13)(b)上部是由所得反褶积算子计算的最小平方逆，下部是输入子波与反射数褶积所得的纪录道集，图(5-13)(c)是最小熵反褶积的输出结果，上部是所求的线性反褶积算子。值得一提的是反褶积算子的最小平方逆与模型的输入子波形非常相似，这说明，最小熵反褶积方法的效果令人满意。第五节极大似然反褶积(MostlikelihoodDeconvolution)极大似然反褶积是建立在概率论中极大似然估计理论基础之上的一种反褶积新技术。由于反褶积准则是建立在概率基础之上，不需假设子波的最小相位特征，只需假设反射系数满足高斯-贝努力序列条件，而该条件符合实际情况，因此，极大似然反褶积方法适应性强。且所得每一个反射脉冲都比较准确地对应于反射界面，地震剖面的精度和分辨率均得到提高。由于极大似然反褶积算法比较复杂，这里只介绍基本原理。本节要点：极大似然反褶积原理；分块最优化迭代法。一、极大似然反褶积原理根据地震记录褶积模型，对于任意一地震道记录()xt，可以写成下式()()()()()()xtstntbttnt现对于波()bt和反射系数()t做如下假设1．对于任意相位的地震子波()bt和反射系数()t，均可由ARMA（自由归平均）模型表示，即子波()bt的Z变换()BZ为12121111()nnnnmmmmbZbZbZbBZZaZaZa（5-5-1）其中，12,,,maaa，12,,,nbbb，m,n为子波参数2．假设，地层的反射系数序列()t满足高斯—贝努力序列，且可用乘积模型表示，则tqtrtt（5-5-2）其中()t是与地层散射序列有关的背景散射序列，()rt是零均值的高斯白噪，它表示反射系数幅度的大小，qt是01序列，它表示反射是否存在。1qt，表示有反射，0qt，表示无反射，且其概率密度满足()11()0rqtPqtqt当当（5-5-3）其中，称为反射频次对于地震记录的褶积模型()()()()xtbttnt，若江地震记录写到向量形式，则有XBn(5-5-4)其中12XcolxxxN00001(0)00(1)(2)(0)BBBBBNBNB1,2colN1,2ncolnnnNcol表示列向量，N表示记录长度。对于(5-5-4)式，可以进一步表示成下式XBQrBn(5-5-5)1,2QdiagqqqN1,2rcolrrrN(1),2colN为使表达式简洁，对所有的变量均以列向量的形式表示，则有12,macolaaa12,nbcolbbb12,nqcolqqq,,,rnScolVVV其中，S表示统计参数列向量，V表示方差，rV表示变量rt的方差，V表示t的方差，nV表示()nt的方差，为反射频次。设地震记录X的似然函数为,,,,,|Labqrsx，根据概率理论，似然函数可以表示成下式,,,,,|(,,,|,,)LabqrsxPxqrabs(5-5-6)其中,P表示概率密度，由条件概率的乘法定理，有,,,,,|)(,,,|,,LabqrsxPxqrabs|,,,,(|,,)PXabqsPqabs(|,,)(|,,)PrabsPabs(5-5-7)若设,,ab与,qr无关，则(5-5-7)式可写为(,,,,,|)(|,,,,,)LabqrsxPXabqrs(|)(|)(|)PqsPrsPs(5-5-8)现在计算上式中的似然函数。根据概率分布性质，有2(|,,,,,)(2)exp{NnPXabqrsVXBQrB()/2}XBQrBV(|)(|,,)rPqsPqabs1(()|)NrKPqK(1)Nmqmq其中1()()NKmqqK(5-5-10)2(|)2exp/2NrrPrsVrrV(5-5-11)2(|)(2)expNPsV(5-5-12)上面各式中括号外的＂＇＂表示转置将(5-5-9)式～(5-5-1)式代入(5-5-8