SRME原理.docx
《SRME原理.docx》由会员分享，可在线阅读，更多相关《SRME原理.docx（24页珍藏版）》请在冰豆网上搜索。
SRME原理
4自由界面多次波预测与压制（SRME原理）
利用f-x域波场外推理论，研究在不考虑自由界面反射的情况下，有效波波场的传播及反射规律；然后，在考虑自由界面反射的情况下，详尽推导了自由界面多次波产生的数学模型，给出了任意阶自由界面多次波波场模拟的级数形式和迭代形式，推导作为自由界面多次波产生的逆过程的自由界面多次波消除的数学模型，深刻剖析自由界面多次波的产生和消除的物理机制。
4.1有效波波场（无自由表面多次波）的模拟
根据第三部分的理论推导，地震波在介质中的传播和反射过程可以用f-x域
波场延拓来表示［82-84］，其几何表示如图4-1所示。
检波点坨（%）爲源亍匕）
图4-1地震波的传播和反射的物理模型
用系统理论可以表示为图4-2。
它把地震波场分解为五个过程：
发射，下行
传播，反射，上行传播，接收
（和知）＜―
图4-2地面波场的正演模拟（无自由界面的反射）
根据Berkhout的研究，对于地下的某一个反射界面m其对应的有效波波场
（f-X域）Po_（Zm，Zo）为[108]
如果忽略自由反射界面，则RZo=o。
无自由表面接收到的总的有效波波场P°—（Zo）[111]为
代入式（4-1），得
P（fZo—W~Zo,ZmRZmWZm,Z°SZo（4-2）
m吕
进而
PoZo=XoZo,XoSZo（4-3）
M
Xo（Zo,Zo）=vW-（Zo.Zm）R（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的传播。
对于算子的每一元素可以表示为
式中，HIQ]限带微分算子，.|j为从炮点xj,Zo到反射点X|,Zm的旅行时，・ik为
从反射点Xk,Zm到接收点Xi,Zo的旅行时，aij,aik球面扩散因子；
算子R（Zm）为界面Zm处的下行入射波场的反射算子。
无自由表面多次波的时间域表达式为m2]
PjZo=7Zo,ZmrZmWZm,Z。
!
“SjZ。
（4-5）
m
式中“为kronecker内积。
4.2自由界面多次波正演模拟
图4-3有效波场和多次波场的正演模拟，包含自由界面的反射
当不考虑自由界面的反射时，在Z0界面处的上行波场为
P°—（z。
）=X0（Z°,Ze）S（z。
）（4-6）
当考虑自由界面的反射时，上行波场P"（Zo）遇自由界面Zo后发生反射，变为下行波场。
因此总的下行波场不仅包含下行的震源波场S（zo），而且还应包含被自由界面反射回地下的波场：
R一（Zo）p-（Zo），如图4-3所示。
R—（Zo）为自由界面对
于上行波场的反射算子，P-（Z。
）为自由界面处的上行波场。
因此自由界面Zo处
的总的上行波场为：
P°-（Zo）=Xo（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（Zo）=D-（Zo）Xo（Zo,Zo）[S（Zo）R-（Zo）P-（Zo）]（4-9）
可以将自由表面接收到的波场分解为有效波场Po（Zo）和（自由界面）多次波波场
Mo（Zo），即
Po（Zo）=D-（Zo）Xo（Zo，Zo）S（Zo）（4-10）
Mo（Zo）=D-（Zo）Xo（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）]4R-（Zo）[D"（Zo）]（4-12b）
Berkhout称A（z°）为自由界面算子。
由代数运算可得
[I-Po（Zo）A（Zo）]P（Zo）=Po（Zo）宀
（4-13）
P（Zo）=[I-Po（Zo）A（Zo）]4Po（Zo）
用Neumanr级数展开得[112]
=Po（Zo）、Po（Zo）A（Zo）4o（Zo）Po（Zo）A（Zo）2Po（Zo）
Po（Zo）A（Zo）2Po（Zo）Po（Zo）A（z°）n（4-14）
由通项Tn讯Po（Zo）A（Zo）}nPo（Zo）可知：
当n=o时，Tn={Po（zo）A（Zo）}Po（zo），即为1阶自由界面多次波;
当n=n时，T={Po（z°）A包）}*®。
），即为n阶自由界面多次波
自由界面多次波波场模拟的方法可以分为以下两种。
（1）截断离散级数方法
P（Zo）=Po（Zo）r[Po（Zo）A（Zo）]n}Po（Zo）（4-15）
n4
级数的次数就是多次波的阶数。
首先模拟出有效波场，然后叠加上多次波场，就得到总的波场。
具体展开可得
当n=1时，P二P。
-PoAPo
当n=2时，P=P°P0AP0+P0AP0AP0
当n=3时，P二PoPoAP0+PoAPoAR+PoAPoAP。
AP。
当n=n时，P=P°PoAPo+PoAPoAPo+PoAPoAPoAPo.…（PoA）nPo
在实际的多次波模拟中，选取有限项。
（2）迭代方法
由P（Zo）=Po（Zo）Po（Zo）A（Zo）P（Zo），其迭代表达式为
（4-16）
；P（n41）=Po+PoAP（n）
P（o）=Po
其中，n为迭代次数，P（n1）是第n+1次迭代得结果，Po是有效波波场
从本质上讲，迭代算法和级数算法是一致的，只不过两者具体的实现形式不一样;迭代形式的表达更加简洁，物理含义更加明确（po为自由界面多次波的预测算子，A为自由界面算子），在计算机上更易实现，实现起来也能节约计算机资源。
4.3自由界面多次波的压制
（1）多次波迭代压制理论基础[1o6]
自由界面多次波的消除过程是自由界面多次波产生过程的反过程
（4-17a）
P（Z°）=Po（Zo）Po（Zo）A（Zo）P（Zo）
（4-17b）
A（zo）=[S（zo）]-R-（z°）[D-（zo）]J
经过矩阵变换：
Po（Zo）=P（Zo）[IA（Zo）P（Zo）]」二[IP（Zo）A（Zo）]JP（Zo）=
P（Zo）-{p（ZoAZoj）P（Zo）+{p（ZoA（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（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）
写成迭代形式为
POn1）（Zo^P（Zo）-Po（n）（Zo）An1（Zo）P（Zo）（4-22）
式中，n为迭代次数，P（zo）为原始的地震输入数据（含多次波），POn）（Zo）为第n次无多次波记录的估计值，POn1）（zo）为第n+1次迭代的结果。
对于每次迭代，我们可以使用不同的自由界面算子A（n1）（Zo）。
自由表面多次波压制的Neumanr级数表达式为
Pozo二PZO-1n4PZoAZoFPZo（4-23）
nm
qQ
令Fzo-7-1n4PzoAzo片，贝U
n=1
PoZo=PZo-FZoPZo二IFZoPZo（4-24）
根据多道Wiener预测滤波理论，Fz。
可以看作基于模型的预测算子。
在算子
FZo中，只需估计自由表面算子AZo。
从自由表面算子AZo定义来看，如果
1
假设R~Zo--1，则AZo[=[D。
如果假设地震记录与炮检对方向特性无
关，则算子AZo就成对角矩阵。
在地表一致的假设前提条件下，算子AZo的对角元素将相同，即单位矩阵的倍数，对角元素序列表示Fourier域的反子波。
根据Verschuur的工作[io8][112]，估算算子Azo的目标函数为
E〈PZo-FNZoPZoU（4-25）
N
式中，FNZo1ndPZoAZon，E「｝为地震记录的所有频率成分的平均
响应。
该目标函数为非线性的。
当N「：
时，FnZoAPoZoAZo。
这意味着自由表面多道多次波预测算子为反褶积后的有效波响应。
为了简化计算，作以下假设：
1震源子波不发生改变；
2自由界面的反射系数为-1；
3在数据采集过程中，检波器保持稳定。
则自由界面算子简化为
A（Zo）=-[s（）]JI二A（）（4-26）
为一个只与频率有关的对角矩阵，使算子FnZo待求系数大大减少。
因此，式（4-24）变为
Po（nd）（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）理论验证
现在假设平面波在水平层状介质中传播，并且只有一个反射界面（Z1），则Zo
处的上行波场为：
式中I；=2乙-Zo/c，r为反射系数。
因此AiQ卜S'-
因此有
预测算子F••的级数表示为
（4）算法具体实现步骤
消除与自由界面有关多次波的具体实现步骤是：
1预测出与自由界面有关的多次波波场，数据选排方式如图4-4。
为了增加覆
盖次数（Kirchhoff积分孔径宽度）还需对称拷贝数据（即模拟中间放炮的结果）⑹
mS,R,-八psx,•Prx,•（4-29）
x
式中，PsFtJtPsX,tI
F为傅立叶变换，V为水速。
图4-4f-X域数据选排方式（平面图）
2作关于频率参量的反傅立叶变换，把M返回到时间空间域
mS,R,tA—t"F.JmS,R,•1（4-30））
式中S,R分别为检波点坐标和震源坐标，t为时间变量。
3在每个炮集记录上，根据能量最小L2模，估计一个滤波算子an1t，使得从
输入数据中减去预测出的多次波数据后所得到的能量最小，该公式可以写成
EdJmlp（t,Xr,x^）-ai（n+）（t^m（n+it,R,S（i^（4-31）
t,Xr
其中i表示炮数，*表示褶积。
这一步可用标准的最小平方方法（维纳滤波）进行
求解。
4奇异值分解（SVD）滤波重构的原理
Beltrami和Jordan二位学者是奇异值分解的主要创始人：
Beltrami于1873
年发表了奇异值分解的第一篇论文［50］，一年后Jordan发表了自己对奇异值分解的独立推导［51］。
后来,Autonne把奇异值分解推广到复正方矩阵［52］；Eckart与Young又进一步把它推广到一般的长方形矩阵［53］。
尽管对奇异值分解方法的理论研究始于19世纪70年代，但其真正被广泛应
用却是在1970年以后。
电子计算机的飞速发展推动了数值计算方法的研究，
Golub等人首先提出了可以在计算机上实现的奇异值分解算法，这就为奇异值分
解方法的广泛应用提供了可能网屈。
现在，奇异值分解（包括各种推广）已经是数值线性代数的最有效的工具之一，它在统计分析、信号与图像处理、系统理论
和控制中被广泛地应用。
4.1矩阵的奇异值分解
奇异值分解最早是由Beltrami在1873年对实正方矩阵提出来的。
Beltrami
从双线性函数
fx,y=xTAy，ARnn
出发，通过引入线性函数
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二U'VT（
这就是Beltrami于1873年得到的实正方矩阵的奇异值分解。
4.1.1矩阵的奇异值分解及其解释
任何一个mn矩阵A有正交分解
A二HRKt
其中
Rn0
「0
R1为k^k矩阵，秩为k；H、K为正交矩阵
现在我们在此基础上做进一步分解。
首先对满秩矩阵Rii做分解。
矩阵RIr,"是对称的正定矩阵，因此由它的特征值和特征向量可得到分解
R^Ri=ViiDiiV|1（4-1-5）
其中Vii为正交矩阵，Dii为对角矩阵，对角线上的元素为正的并要求为不减的。
设％为对角矩阵，它的对角线上的元素为Di,的相应的元素的正平方根。
于
是有
对R,i做分解
其中
4-1-8）
Uii=R1V11S1
由公式（4-1-5）和（4-1-6）可得
因此，U,,为kk矩阵。
由式（4-1-4）和（4-1-7）可知
把式（4-1-10）代入（4-1-3）中可得
由上式不难得到下面的定理:
令ARmn，则存在正交矩阵U•Rmm、VRnn和对角矩阵a使得
4-1-12）
A=U'VT
式中
瓦七杲（4-1-13）
00
'i=diag（、i，、2,………）,其对角元素按照顺序d^>2>……>6^0,
rankA排列。
这就是矩阵的奇异值分解定理，称:
.（i=1,2，…，r—1,r）为矩阵A的奇异值，
A二U'Yt称为矩阵A的奇异值分解式[56][57]。
下面是关于矩阵的奇异值和奇异值分解的几点解释：
（1）矩阵A的奇异值分解式（4-1-12）可以改写成向量表达形式：
r
A^:
iUivH（4-1-14）
i4
这种表达式称为A的并向量（奇异值）分解（dyadicdecomposition）。
（2）当矩阵A的秩r=ra（rA）k-：
r彳=：
r.2二…=-F=0，h=min，m,n，，故奇异值分解公式（4-1-12）可以简化
为
A=Ur'rVrH（4-1-15）
式中
Ur=U1,U2,,u」，Vr=V1,V2,,v」，'r=diagC1,、2,,'r）
式（4-1-15）称为矩阵A的截尾奇异值分解（truncatedSVD）；与之形成对照，
式（4-1-12）则称为全奇异值分解（fullSVD）。
（3）令ARmn（m•n）的奇异值为
1一、2--、n-0（4-1-16）
则
^min/lEf：
rank（AK）^k-1』,k=1,2；,n（4-1-17）
并且存在一个满足||Ek|F=6的误差矩阵E使得
4-1-18）
rank（AEk）二k-1,k=1,2,,n
由以上分析可知，当原nn矩阵A有一个零奇异值时，该矩阵的秩rankA）空n-1，即原矩阵A本来就不是满秩矩阵。
因此，如果一个正方矩阵具有零奇异值，则该矩阵必定是奇异矩阵；一个正方矩阵只要有一个奇异值接近于零，那么这个矩阵就接近于奇异矩阵。
推而广之，一个非正方的矩阵如果有奇异
值为零，则这个长方矩阵一定不是列满秩的或者行满秩的。
这种情况称为矩阵的
秩亏缺，它相对于矩阵的满秩是一种奇异现象。
总之，无论是正方矩阵还是长方
矩阵，零奇异值都刻画了矩阵的奇异性。
这就是矩阵奇异值的内在涵义。
4.1.2矩阵奇异值的性质
（1）奇异值与特征值的关系
由式（4-1-12）可得
AAh二U'2Uh（4-1-19）
这表明，mn矩阵A的奇异值「是矩阵乘积AAh的特征值打（这些特征值
是非负的）的正平方根,即；打=.'j，（i=1,2,,r-1,r），■i（i=1,2/,r-1,r）。
（2）奇异值与范数的关系昭
矩阵A的谱范数等于A的最大奇异值，即
Aspect（4-1-20）
根据矩阵的奇异值分解定理，并注意到矩阵A的Frobenius范数AF是正交
不变的，即|uhav|fTIAIf，故有
-mn2f2
IAf十送aj
（4-1-21）
7”
八F
也就是说，任何一个矩阵的Frobenius范数等于该矩阵的所有非零奇异值平方和的正平方根。
考虑矩阵A的秩k近似，并将其记作A，其中k：
：
：
r二rank（A）。
矩阵A定义
如下：
k
Ak-'、山2「,k:
:
:
r（4-1-22）
14
则A与秩为k的任一矩阵B之差的l1和Frobenius范数分别为
mina_Bi二A-人i-宀（4-1-23）
rank（B）主
和
22222
min|A-BF=||A—AkF=朿半*642+…+6（4-1-24）
rank（B）_k
这一重要结果是许多概念和应用的基础，它就为用一个低秩矩阵逼近原始矩
阵提供了理论依据。
4.1.3矩阵奇异值分解的展开式
根据奇异值分解定理我们不难证明，在A的奇异值分解A二U'7T中的U的
列和V的行分别为AAt和ATA的特征向量。
这是因为
AAT二U'VTV'UT二U'2UT（4-1-25）
即
（AAT）U二Udiag（'1,'2,,r,0,,0）（4-1-26）
将U按列分块，并记U-（u1,u2/,un4,un）
代入式（4-1-26）可得
（AAT）Ui=iUii=1,2/,n—1,n（4-1-27）
这就说明ui是AAt的属于打的特征向量。
同理，因为
ATAMUTU'VTW2VT（4-1-28）
所以，y是ATA属于打的特征向量。
由上式可知协方差矩阵aat和ata的秩是相等的，所以非零特征值也相等。
由于所以式（4-1-12）可写成
r
A八（4-1-29）
iJ
我们称UiViT为A的第i个特征图像。
同时由上式的组成可以看出特征图像在
重构A中的作用是正比于其奇异值「的大小，而:
.是按大小顺序排列的，所以最
前面几个特征图像在重构A中所占的比重最大。
将上式进行图像分解如下：
我们注意到，UivT是秩为1的简单的mn矩阵，需要储存单元为m•n个，由上
式可知，按这种方法储存A需要的储存单元为r（m•n），而原矩阵需要的储存
单元为mn个，当r很小时，r（mn）比mn小的多。
因此用最佳逼近矩阵A代替原矩阵可以极大地节省储存单元。
但是，更为重要的是最佳逼近矩阵A反映了原矩阵中最重要的信息，所以矩阵的奇异值分解方法在二维数据处理中得到了广泛地应用［58］。
4.2矩阵的奇异值分解可用于地震数据处理的理论依据
设有N道地震记录，每道采样点数为M，那么，这NM个数据可以构成一个二维图像，用矩阵表示为Xnm。
现在从数学上说明奇异值分解可用于地震资料处理的理论依据。
首先让我们来证明
X-Xp||；八J（4-2-1）
式中Xp为奇异值6p十,―七,…，q取值为0时的重构矩阵；|x|F为矩阵X的F范
任意一个NM矩阵X的F范数定义为
MN
XF=Xi2）12（4-2-2）
yjq
显然，地震记录数据的总能量为
X；=tr1XX丁；=trlxTX?
（4-2-3）
N
其中，tLA；二為aii称为方阵A的迹，它等于方阵A的对角线元素之和，迹的另i4
一个性质是
trlAB—tr「BA?
（4-2-4）
取最初的p个特征图像重构X可以表示为
Xp=UEpVt（4-2-5）
其中，Ep二diag（、i,、2,,、p,0,,0）
由式（4-2-4）和（4-2-5）可得
22
X-Xp|L二UEVT_UEpVTF
TI2
=U（E—Ep）V|f
=trU（E-Ep）VTV（E-Ep）U”
=tr°（E-Ep）（E-Ep）UT?
（4-2-6）
=tr"E_Ep）UTU（E_Ep）?
=tr1（E-Ep）（E-Ep*
r
八「i2
i斗1
因此可得
x-Xp；=
r
—i2i耳1
（
4-2-7）
显然，式（4-2-7）中若取Xp=0，则有
r
x：
—i2
i彳
（
4-2-8）
由此可见，
地震记录数据的总能量等于奇异值
-i的平方和，
且奇异值愈大的
分量在地震信号中功率贡献也愈大，这就为数据压缩提供了可能。
从减小误差平方和的观点看，丢掉一些较小的奇异值，所产生的误差较小，由于「的值是单调递减的，所以选用的前p个特征图像与前p个大奇异值相对应。
因此，这样重构的地震记录Xp与实际地震记录X的误差最小[59][60]。
4.3奇异值分解滤波重构原理（K_L变换）
4.3.1奇异值分解滤波原理
其表示成
X经奇异值分解后按能量大小分成若干个特征图像，其分解式可以写成
r
（4-3-2）
X八：
iUivT
i=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。
这种情况下等于原数据体的
升级会员 