SRME原理.docx

1、SRME原理4自由界面多次波预测与压制 （SRME原理）利用f-x域波场外推理论，研究在不考虑自由界面反射的情况下，有效波波 场的传播及反射规律；然后，在考虑自由界面反射的情况下，详尽推导了自由界 面多次波产生的数学模型，给出了任意阶自由界面多次波波场模拟的级数形式和 迭代形式，推导作为自由界面多次波产生的逆过程的自由界面多次波消除的数学 模型，深刻剖析自由界面多次波的产生和消除的物理机制。4.1有效波波场（无自由表面多次波）的模拟根据第三部分的理论推导，地震波在介质中的传播和反射过程可以用 f-x域波场延拓来表示82-84，其几何表示如图4-1所示。检波点坨（%） 爲源亍匕）图4-1地震波的

2、传播和反射的物理模型用系统理论可以表示为图 4-2。它把地震波场分解为五个过程：发射，下行传播，反射，上行传播，接收（和知）图4-2地面波场的正演模拟（无自由界面的反射）根据Berkhout的研究，对于地下的某一个反射界面m其对应的有效波波场（f-X 域）Po_（Zm，Zo）为108如果忽略自由反射界面，则 R Zo = o。无自由表面接收到的总的有效波波场 P（Zo） 111为代入式(4-1)，得P(fZo W Zo,Zm R Zm W Zm,Z S Zo ( 4-2)m吕进而Po Zo = X o Zo , Xo S Zo ( 4-3)MXo(Zo,Zo) = v W -(Zo.Zm)R

3、(Zm)W 久，勺) (4-4)md式中，算子Xo （Zo , Zo）为自由界面无反射的地下介质的脉冲响应;算子P 一仗0）为总的上行波场;算子P （Zo）为总的下行波场;算子PoPZo）为在假设自由界面不反射情况下，在自由界面 Zo处的上行波场;算子R （Zm）为在自由界面Zm处，对于下行入射波场的反射算子；算子S （Zo）为在自由界面Zo处的下行震源波场，包含激发方式和震源的虚 反射效应；算子W -（Zm，Zo）为上、下行传播矩阵，每一列对应于地表上一个点的传播算子，整个矩阵描述了波场从Zm到Zo的传播。对于算子的每一元素可以表示为式中，H IQ 限带微分算子，.|j为从炮点x j ,Zo

4、到反射点X|,Zm的旅行时，ik为从反射点Xk,Zm到接收点Xi ,Zo的旅行时，aij,aik球面扩散因子；算子R （Zm）为界面Zm处的下行入射波场的反射算子。无自由表面多次波的时间域表达式为m2P j Zo =7 Zo,Zm r Zm W Zm,Z。!“ S j Z。 （ 4-5）m式中“为kronecker内积。4.2自由界面多次波正演模拟图4-3 有效波场和多次波场的正演模拟，包含自由界面的反射当不考虑自由界面的反射时，在Z0界面处的上行波场为P（z。）= X0（Z,Ze）S （z。） （4-6）当考虑自由界面的反射时，上行波场 P（Zo）遇自由界面Zo后发生反射，变为下 行波场。因

5、此总的下行波场不仅包含下行的震源波场 S （zo），而且还应包含被自 由界面反射回地下的波场：R 一（Zo）p-（Zo），如图4-3所示。R （Zo）为自由界面对于上行波场的反射算子，P -(Z。)为自由界面处的上行波场。因此自由界面 Zo处的总的上行波场为：P-(Zo) = X o(Zo,Zo) S (Zo) R -(Zo) P -(Zo) (4-7)在实际地震资料中，还应该考虑检波器的特性 (例如检波器组合模式、虚反射特性等)，将检波器属性定义为算子D _(z0)，则在自由界面z0处波场为：P(Zo) = D -(Zo)P -(Zo) (4-8)将P -(Zo)代入式(4-8 )，得P(Z

6、o) = D -(Zo)X o(Zo,Zo) S (Zo) R -(Zo) P -(Zo) (4-9)可以将自由表面接收到的波场分解为有效波场 Po(Zo)和(自由界面)多次波波场M o(Zo)，即Po(Zo) = D -(Zo)Xo(Zo，Zo)S (Zo) (4-10)M o(Zo) = D -(Zo)X o(Zo,Zo)R -(Zo)P -(Zo) (4-11)由式(4-7)-( 4-11)，得P(Zo) = Po(Zo) Po(Zo) A (Zo) P(Zo) (4-12a)A(Zo) = S (Zo) 4 R -(Zo) D (Zo) (4-12b)Berkhout称A(z)为自由界

7、面算子。由代数运算可得I -Po(Zo)A(Zo)P(Zo) = Po(Zo) 宀( 4-13)P(Zo) =I -Po(Zo) A (Zo)4 Po(Zo)用Neumanr级数展开得112=Po(Zo) 、Po(Zo)A(Zo) 4o(Zo) Po(Zo)A(Zo) 2Po(Zo)Po(Zo) A (Zo)2 Po(Zo) Po (Zo) A (z)n (4-14)由通项 T n 讯 Po(Zo) A (Zo)n Po(Zo)可知：当n= o时，Tn =Po(zo)A(Zo)Po(zo)，即为1阶自由界面多次波;当n=n时，T =Po(z)A包)*。)，即为n阶自由界面多次波自由界面多次波波

8、场模拟的方法可以分为以下两种。(1)截断离散级数方法P(Zo) = Po(Zo) r Po(Zo)A(Zo)nPo(Zo) (4-15)n 4级数的次数就是多次波的阶数。首先模拟出有效波场，然后叠加上多次波场，就 得到总的波场。具体展开可得当 n= 1 时，P 二 P。- PoAPo当 n = 2 时，P=P P0AP0+ P0AP0AP0当 n = 3 时，P 二 Po Po AP0 + PoAPoA R + Po APoAP。AP。当 n = n 时，P =P PoAPo + PoAPoAPo + PoAPoAPoAPo .(PoA)nPo在实际的多次波模拟中，选取有限项。(2)迭代方法由

9、 P(Zo) = Po(Zo) Po(Zo)A(Zo)P(Zo)，其迭代表达式为(4-16)；P(n41)= Po + Po AP(n)P(o)= Po其中，n为迭代次数，P(n 1)是第n+1次迭代得结果，Po是有效波波场从本质上讲，迭代算法和级数算法是一致的，只不过两者具体的实现形式不一样; 迭代形式的表达更加简洁，物理含义更加明确 (po为自由界面多次波的预测算 子，A为自由界面算子)，在计算机上更易实现，实现起来也能节约计算机资源。4.3自由界面多次波的压制(1)多次波迭代压制理论基础1o6自由界面多次波的消除过程是自由界面多次波产生过程的反过程(4-17a)P(Z) = Po(Zo)

10、 Po(Zo) A (Zo) P(Zo)(4-17b)A(zo)=S (zo)- R -(z)D -(zo) J经过矩阵变换：Po(Zo)= P(Zo) I A (Zo) P (Zo)二I P(Zo) A (Zo)J P (Zo) =P(Zo )-p(Zo AZo j)P(Zo )+ p(Zo A(Zo 丫P(Zo) 9(Z)A(z。P(z)+(4-18) 由于原始数据中存在强多次波，该级数有能量不收敛，从而导致消去多次波的反 演过程的不稳定。自由界面算子A(Zo)A(Zo) = s (Zo)R -(Zo)D -(Zo) J (4-19)在多次波的消除过程中，起着自适应滤波作用。只有经过 A(

11、Zo)自适应滤波后， 预测出的多次波P(Zo) 2、P(Zo) 3，才能与数据中实际存在的多次波在振幅和 相位上相匹配。正是因为 A(Zo)算子的存在，才使得该算法是自适应的。由等式P(Zo) = Po(Zo) - Po(Zo) A(Zo) P(Zo) (4-20)可得Po(Zo) = P (Zo) - Po(Zo) A (Zo) P(Zo) (4-21)写成迭代形式为POn 1)(Zo P(Zo) - Po(n) (Zo)A n 1 (Zo)P(Zo) (4-22)式中，n为迭代次数，P(zo)为原始的地震输入数据(含多次波)，POn) (Zo)为第 n次无多次波记录的估计值，POn 1)(

12、zo)为第n+1次迭代的结果。对于每次迭代， 我们可以使用不同的自由界面算子 A(n1)(Zo)。自由表面多次波压制的Neuma nr级数表达式为Pozo 二 P ZO -1 n4 P Zo A Zo F P Zo (4-23)nmqQ令 F zo - 7 -1 n4 P zo A zo 片，贝Un =1Po Zo = P Zo -F Zo P Zo 二 I F Zo P Zo (4-24)根据多道 Wiener预测滤波理论，F z。可以看作基于模型的预测算子。在算子F Zo中，只需估计自由表面算子 A Zo。从自由表面算子A Zo定义来看，如果1假设R Zo - -1，则A Zo=D 。如果

13、假设地震记录与炮检对方向特性无关，则算子A Zo就成对角矩阵。在地表一致的假设前提条件下，算子A Zo的对 角元素将相同，即单位矩阵的倍数，对角元素序列表示 Fourier域的反子波。根据Verschuur的工作io8112，估算算子A zo的目标函数为EP Zo -FN Zo P Zo U (4-25)N式中，FN Zo 1 nd P Zo A Zo n，E为地震记录的所有频率成分的平均响应。该目标函数为非线性的。当N：时，Fn Zo A Po Zo A Zo。这意味着自由表面多道多次波预测算 子为反褶积后的有效波响应。为了简化计算，作以下假设：1震源子波不发生改变；2自由界面的反射系数为-

14、1 ；3在数据采集过程中，检波器保持稳定。则自由界面算子简化为A(Zo) = -s ( ) JI 二 A( ) (4-26)为一个只与频率有关的对角矩阵，使算子 F n Zo待求系数大大减少。因此，式(4-24)变为Po(n d)(Zo P (Zo)- A (n1)( JPo(n)(Zo) P(Zo) (4-27)将式(4-27)改写为迭代算法迭代形式P0n+)(zo P - A(n41) P0n) P (4-28a)(4-28b)P0o)二 Po 式中，P为原始波场（含多次波）；P0n）P为预测出的多次波波场；Po（n）为多次波的预测算子；K 为自由界面算子，起着自适应滤波算子的作用。（2）

15、 理论验证现在假设平面波在水平层状介质中传播，并且只有一个反射界面 （Z 1），则Zo处的上行波场为：式中I； = 2乙-Zo /c，r为反射系数。因此AiQ卜S-因此有预测算子F 的级数表示为（4）算法具体实现步骤消除与自由界面有关多次波的具体实现步骤是：1预测出与自由界面有关的多次波波场，数据选排方式如图 4-4。为了增加覆盖次数（Kirchhoff 积分孔径宽度）还需对称拷贝数据（即模拟中间放炮的 结果）m S,R, - 八 ps x, Pr x, （4-29）x式中，Ps Ft JtPs X,t IF为傅立叶变换，V为水速。图4-4 f -X域数据选排方式(平面图)2作关于频率参量的反

16、傅立叶变换，把 M返回到时间空间域m S, R,t A t F. Jm S,R, 1 (4-30) )式中S,R分别为检波点坐标和震源坐标，t为时间变量。3在每个炮集记录上，根据能量最小 L2模，估计一个滤波算子an1 t，使得从输入数据中减去预测出的多次波数据后所得到的能量最小，该公式可以写成EdJm lp(t,Xr,x )-ai(n+)(tm(n+it,R,S(i (4-31)t,Xr其中i表示炮数，*表示褶积。这一步可用标准的最小平方方法 (维纳滤波)进行求解。4奇异值分解(SVD )滤波重构的原理Beltrami和Jordan二位学者是奇异值分解的主要创始人： Beltrami于187

17、3年发表了奇异值分解的第一篇论文50，一年后Jordan发表了自己对奇异值分解 的独立推导51。后来,Auto nne把奇异值分解推广到复正方矩阵52 ； Eckart与You ng 又进一步把它推广到一般的长方形矩阵53。尽管对奇异值分解方法的理论研究始于 19世纪70年代，但其真正被广泛应用却是在1970年以后。电子计算机的飞速发展推动了数值计算方法的研究，Golub等人首先提出了可以在计算机上实现的奇异值分解算法， 这就为奇异值分解方法的广泛应用提供了可能网屈。现在，奇异值分解（包括各种推广）已经是 数值线性代数的最有效的工具之一， 它在统计分析、信号与图像处理、系统理论和控制中被广泛地

18、应用。4.1矩阵的奇异值分解奇异值分解最早是由 Beltrami在1873年对实正方矩阵提出来的。Beltrami从双线性函数f x, y = xT Ay， A Rn n出发，通过引入线性函数x = U ， y = V将双线性函数变为式中S=UTAV （ 4-1-1）Beltrami观测到，如果约束U和V为正交矩阵，则它们的选择各存在 n? 一 n个 自由度。他提出利用这些自由度使矩阵S的对角线以外的元素全部为零，即矩阵S = = diag（, ，:n）为对角矩阵。于是用U和V分别左乘和右乘式（4-1-1），并利用U和V的正交性，立即得到4-1-2)(4-1-3)4-1-4)A二UVT （这就

19、是Beltrami于1873年得到的实正方矩阵的奇异值分解。4.1.1矩阵的奇异值分解及其解释任何一个m n矩阵A有正交分解A 二 HRKt其中Rn 00R1为kk矩阵，秩为k ； H、K为正交矩阵现在我们在此基础上做进一步分解。 首先对满秩矩阵Rii做分解。矩阵RIr, 是对称的正定矩阵，因此由它的特征值和特征向量可得到分解RRi =ViiDiiV|1 （ 4-1-5）其中Vii为正交矩阵，Dii为对角矩阵，对角线上的元素为正的并要求为不减的。设为对角矩阵，它的对角线上的元素为 Di,的相应的元素的正平方根。于是有对R,i做分解其中4-1-8)Uii = R1V11S1由公式（4-1-5）和

20、（4-1-6）可得因此，U,为k k矩阵。由式（4-1-4）和（4-1-7）可知把式（4-1-10）代入（4-1-3）中可得由上式不难得到下面的定理:令A Rm n，则存在正交矩阵U Rm m、V Rn n和对角矩阵a使得4-1-12)A =U VT式中瓦七杲 ( 4-1-13)0 0 i =diag(、i，、2,), 其对角元素按照顺序 d260 ,rank A排列。这就是矩阵的奇异值分解定理，称:.(i =1,2，r 1,r)为矩阵A的奇异值，A二U Yt称为矩阵A的奇异值分解式5657。下面是关于矩阵的奇异值和奇异值分解的几点解释：(1)矩阵A的奇异值分解式(4-1-12)可以改写成向量

21、表达形式：rA :iUivH ( 4-1-14)i 4这种表达式称为A的并向量(奇异值)分解(dyadic decomposition)。(2)当矩阵 A的秩r = r a (rA)k m inn n时，由于奇异值-：r彳=：r .2二=-F = 0，h = min，m, n，故奇异值分解公式(4-1-12)可以简化为A=UrrVrH ( 4-1-15)式中Ur= U1,U2, ,u，Vr= V1,V2, ,v， r = diagC 1,、2, , r)式(4-1-15)称为矩阵A的截尾奇异值分解(truncated SVD)；与之形成对照，式(4-1-12)则称为全奇异值分解(full SV

22、D)。(3)令A Rm n(m n)的奇异值为1 一、2 - -、n -0 (4-1-16)则min/lE f ： rank(A K)k-1, k=1,2； ,n ( 4-1-17)并且存在一个满足|Ek|F =6的误差矩阵E使得4-1-18)rank(A Ek)二 k -1, k =1,2, , n由以上分析可知，当原n n矩阵A有一个零奇异值时，该矩阵的秩 ra nkA)空n-1，即原矩阵A本来就不是满秩矩阵。因此，如果一个正方矩阵具 有零奇异值，则该矩阵必定是奇异矩阵；一个正方矩阵只要有一个奇异值接近于 零，那么这个矩阵就接近于奇异矩阵。 推而广之，一个非正方的矩阵如果有奇异值为零，则这

23、个长方矩阵一定不是列满秩的或者行满秩的。 这种情况称为矩阵的秩亏缺，它相对于矩阵的满秩是一种奇异现象。 总之，无论是正方矩阵还是长方矩阵，零奇异值都刻画了矩阵的奇异性。这就是矩阵奇异值的内在涵义。4.1.2矩阵奇异值的性质(1)奇异值与特征值的关系由式(4-1-12)可得AAh 二U 2U h ( 4-1-19)这表明，m n矩阵A的奇异值是矩阵乘积AAh的特征值打(这些特征值是非负的)的正平方根,即；打=.j ， (i =1,2, , r -1, r )， i (i = 1,2/ , r -1,r )。(2)奇异值与范数的关系昭矩阵A的谱范数等于A的最大奇异值，即Aspect ( 4-1-2

24、0)根据矩阵的奇异值分解定理，并注意到矩阵 A的Frobenius范数A F是正交不变的，即|uhav|f TIAIf，故有-m n 2 f2IAf十送aj( 4-1-21)7”八F也就是说，任何一个矩阵的Frobe nius范数等于该矩阵的所有非零奇异值平方和 的正平方根。考虑矩阵A的秩k近似，并将其记作A，其中k ： r二rank(A)。矩阵A定义如下：kAk - 、山2 , k : r ( 4-1-22)14则A与秩为k的任一矩阵B之差的l1和Frobenius范数分别为min a_Bi 二 A-人 i -宀 ( 4-1-23)rank (B)主和22 2 2 2min |A-BF =|

25、A Ak F =朿半*642+6 ( 4-1-24)rank(B) _k这一重要结果是许多概念和应用的基础， 它就为用一个低秩矩阵逼近原始矩阵提供了理论依据。4.1.3矩阵奇异值分解的展开式根据奇异值分解定理我们不难证明，在A的奇异值分解A二U 7 T中的U的列和V的行分别为AAt和AT A的特征向量。这是因为AAT 二U VTVUT 二U 2UT (4-1-25)即(AAT)U 二Udiag(1,2, , r,0, ,0) ( 4-1-26)将U按列分块，并记U -(u1,u2/ ,un 4,un)代入式(4-1-26)可得(AAT)Ui= iUi i =1, 2/ , n 1, n ( 4

26、-1-27)这就说明ui是AAt的属于打的特征向量。同理，因为ATA M UTU VT W 2VT ( 4-1-28)所以，y是ATA属于打的特征向量。由上式可知协方差矩阵aat和ata的秩是相等的，所以非零特征值也相等。由于 所以式(4-1-12)可写成rA 八( 4-1-29)i J我们称UiViT为A的第i个特征图像。同时由上式的组成可以看出特征图像在重构A中的作用是正比于其奇异值的大小，而:.是按大小顺序排列的，所以最前面几个特征图像在重构A中所占的比重最大。将上式进行图像分解如下：我们注意到，UivT是秩为1的简单的m n矩阵，需要储存单元为m n个，由上式可知，按这种方法储存 A需

27、要的储存单元为r (m n)，而原矩阵需要的储存单元为m n个，当r很小时，r (m n)比m n小的多。因此用最佳逼近矩阵A 代替原矩阵可以极大地节省储存单元。但是，更为重要的是最佳逼近矩阵A反映 了原矩阵中最重要的信息，所以矩阵的奇异值分解方法在二维数据处理中得到了 广泛地应用58。4.2矩阵的奇异值分解可用于地震数据处理的理论依据设有N道地震记录，每道采样点数为 M，那么，这N M个数据可以构成 一个二维图像，用矩阵表示为 Xn m。现在从数学上说明奇异值分解可用于地震 资料处理的理论依据。首先让我们来证明X -Xp|；八 J (4-2-1)式中Xp为奇异值6 p十,七,，q取值为0时的

28、重构矩阵；|x|F为矩阵X的F范任意一个N M矩阵X的F范数定义为M NX F = Xi2)12 (4-2-2)y jq显然，地震记录数据的总能量为X ； = tr 1XX 丁 ； = tr lx T X ? (4-2-3)N其中，tLA；二為aii称为方阵A的迹，它等于方阵A的对角线元素之和，迹的另 i 4一个性质是trlABtrBA? (4-2-4)取最初的p个特征图像重构X可以表示为Xp =UEpVt (4-2-5)其中，Ep 二diag(、i,、2, ,、p, 0, ,0)由式(4-2-4)和(4-2-5)可得2 2X -Xp|L 二 UEVT _UEpVT FT I 2= U(EEp

29、)V |f=tr U(E -Ep)VTV(E -Ep)U”=tr (E -Ep)(E -Ep)UT? ( 4-2-6)=tr E _Ep)UTU(E _Ep) ?=tr 1(E -Ep)(E -Ep*r八i2i斗1因此可得x-Xp； =ri2 i耳1(4-2-7)显然，式（4-2-7）中若取Xp =0，则有rx：i2i彳(4-2-8)由此可见，地震记录数据的总能量等于奇异值-i的平方和，且奇异值愈大的分量在地震信号中功率贡献也愈大， 这就为数据压缩提供了可能。从减小误差平 方和的观点看，丢掉一些较小的奇异值，所产生的误差较小，由于的值是单调 递减的，所以选用的前p个特征图像与前p个大奇异值相对

30、应。因此，这样重构 的地震记录Xp与实际地震记录X的误差最小5960。4.3奇异值分解滤波重构原理(K_L变换)4.3.1奇异值分解滤波原理其表示成X经奇异值分解后按能量大小分成若干个特征图像，其分解式可以写成r(4-3-2)X 八：i Ui vTi=1式中上角T表示转置号，r为X的秩，Ui为XXT的第i个特征向量，Vi为XTX的 第i个特征向量，为X的第i个奇异值，称UivT为X的第i个特征图像。显然，当N道地震数据为线性无关时，它的秩r = N。此时所有的均不为零，因此要完整地重构X就需要把所有的特征图像UivT (i =1,2,，N)进行加权求N和，即X二UiViT。这种情况下等于原数据体的