SRME原理Word文档格式.docx
《SRME原理Word文档格式.docx》由会员分享,可在线阅读,更多相关《SRME原理Word文档格式.docx(24页珍藏版)》请在冰点文库上搜索。
式中“为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处
的总的上行波场为:
-(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<
mi^nnn}时,由于奇异值
-:
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)可写成
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)
其中,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*
八「i2
i斗1
因此可得
x-Xp;
=
—i2i耳1
(
4-2-7)
显然,式(4-2-7)中若取Xp=0,则有
x:
—i2
i彳
4-2-8)
由此可见,
地震记录数据的总能量等于奇异值
-i的平方和,
且奇异值愈大的
分量在地震信号中功率贡献也愈大,这就为数据压缩提供了可能。
从减小误差平方和的观点看,丢掉一些较小的奇异值,所产生的误差较小,由于「的值是单调递减的,所以选用的前p个特征图像与前p个大奇异值相对应。
因此,这样重构的地震记录Xp与实际地震记录X的误差最小[59][60]。
4.3奇异值分解滤波重构原理(K_L变换)
4.3.1奇异值分解滤波原理
其表示成
X经奇异值分解后按能量大小分成若干个特征图像,其分解式可以写成
(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)进行加权求
和,即X二^「UiViT。
这种情况下等于原数据体的