基于谱法lu分解的h的快速求逆运算

基于谱法lu分解的h的快速求逆运算

1维探索中的高效长算子的提出,是使各向异性出现空间在复杂介质的地震波场模型计算中，显式有限差分法和离散方程算法都遇到了困难。由于介质的复杂性，为了确保算法的稳定性，通常有两种方法。其中之一是使用高频散分。然而，这增加了算子的长度，并增加了每一步的计算量。另一种方法是使用低阶算子，但需要在计算中使用更小的时间间隔。为了获得相同的波场，该方法需要按照更小的计算步骤和网格点的数量进行加权。尽管这两种方法最终导致计算量的增加，但隐藏波场算法需要考虑矩阵的逆，但在计算过程中，计算速度非常快，因此计算效率高。如果可以找到快速、准确的矩阵逆方法，隐式算法的优势非常明显。同样,在深度偏移方法中,隐式二维有限差分方法是一种稳健、准确的偏移方法.它能有效地处理介质横向不均匀性.但是,简单地把二维隐式方法推广到三维,在二维方法中可以通过追赶法实现的对角矩阵求逆问题,对应到三维方法中后,就是面对一个分块对角矩阵求逆问题.通常,这种矩阵的求逆将耗费大量计算时间.这一缺陷严重制约了三维隐式方法偏移在实际资料处理中的广泛应用.为了避免这一矩阵的求逆,人们转而采用显式有限差分方法.目前,显式方法在处理三维问题的有效性已经得到了认可.然而,与隐式方法不同,在介质存在横向变化时,显式方法的稳定性无法得到保证.同时,为准确处理陡倾角问题,通常需要较长的显式算子,而这种长算子不仅不能处理横向突变介质,而且需要较长的计算时间.另一种避免这类分块对角矩阵求逆的方法就是算子的分裂算法.通过把算子沿X和Y方向分解,可以避免这类矩阵的求逆.但是,这一做法又导致了传播算子在0°和90°方向的传播速度比45°方向的速度快,即出现了算子在方位角上所谓的各向异性.当然,这种各向异性可以通过增加相位矫正算子或其他方法加以消除.但所有这些方法在消除各向异性时,都不可避免地增加了计算量.因此,直接寻找快捷的矩阵求逆方法成为三维隐式方法实现的关键.最近,Claerbout提出了螺旋坐标系的观点.在螺旋坐标系上,各隐式方法中需要求逆矩阵——亥姆霍兹算子的表示矩阵可化为具有对称带状结构的正定厄密矩阵.即原本在瞬变边界条件下为分块对角的矩阵,在螺旋边界条件下可以表示为具有Toeplitz结构的矩阵.并可以由LU分解方法来实现其快速求逆.最直接的LU分解方法就是直解法,但直解法在完成几万次矩阵分解所需的时间也是不容忽视的.另一种快速的LU分解方法就是谱因式分解方法.具体地说,采用谱因式分解方法,如Kalmogoroff分解法或Wilson-Burg分解法,求得被分解矩阵的每一列ri的最小相位自相关函数ai,把这些自相关函数ai作为列依次排列而得到的矩阵作为分解后的下三角矩阵.同时,根据分解后的下三角矩阵的非零元素在结构上的特点,在计算中采用非零存储和非零计算,最大可能地避免了一些不必要的计算.这种由谱因式分解实现的LU分解方法(简称谱法LU分解)的突出优点是通过建立矩阵列的谱分解表,随后的矩阵分解就可以通过查表快速实现.它的缺陷是分解的不准确性.2波场模拟和深度偏移隐式方法中的矩阵逆问题2.1相空间波场解析地震波在介质中的传播可用如下三维标量波动方程描述Δ2u-1V2∂2∂t2u=0,(1)Δ2u−1V2∂2∂t2u=0,(1)这里V表示介质中地震波传播的速度,u表示地震波场.Δ2Δ2为拉普拉斯算子,对三维空间有Δ2=∂2∂x2+∂2∂y2+∂2∂z2Δ2=∂2∂x2+∂2∂y2+∂2∂z2.若定义v=∂u∂tv=∂u∂t为波场的广义动量,u=u为广义坐标,就得到了一个相空间(v,u),相空间中的波场由相点Z=(v,u)′表示,“′”代表转置.在相空间中,波动方程化为∂Ζ∂t=BΖ,B=[0V2Δ210].(2)∂Z∂t=BZ,B=[01V2Δ20].(2)相空间中的波场从时间t～(t+dt)的解析变换关系为Ζk+1=exp(dtB)Ζk,(3)Zk+1=exp(dtB)Zk,(3)在波场模拟中,上式可以采用隐式辛几何算法.在时间上为二阶精度的隐式辛格式为[CX4]v[CX]k+1=1/V2+dt2Μ1/V2-dt2Μ[CX4]v[CX]k+2dtΜ1/V2-dt2Μuk,uk+1=1/V2+dt2Μ1/V2-dt2Μuk+2dtΜ1/V2-dt2Μ[CX4]v[CX]k,(4)[CX4]v[CX]k+1=1/V2+dt2M1/V2−dt2M[CX4]v[CX]k+2dtM1/V2−dt2Muk,uk+1=1/V2+dt2M1/V2−dt2Muk+2dtM1/V2−dt2M[CX4]v[CX]k,(4)其中矩阵M为算子Δ2Δ2在离散网格上的表示矩阵.2.2频率空间域波场的解析变换关系在基于波动方程的深度偏移中,常把(1)式变换到频率空间域并分解成上行波和下行波.单程波方程可写为∂Ρ∂z=±i√c2-Δ2Ρ,c=ω/V,Δ2=∂2∂x2+∂2∂y2,(5)∂P∂z=±ic2−Δ2−−−−−−−−√P,c=ω/V,Δ2=∂2∂x2+∂2∂y2,(5)这里P为频率波数域的波场,ω为频率,±分别表示下行波和上行波.由(5)式,频率空间域波场从深度z～(z+dz)的解析变换关系为Ρ(z+dz,kx,ky,ω)=Ρ(z,kx,ky,ω)exp(±i√c2-Δ2dz),(6)上式可以采用隐式展开计算,若采用二阶对称Padè近似和45°展开可得Ρ(z+dz,x,y,ω)=1-4c2+icΜ1-4c2-icΜexp(±ic)Ρ(z,x,y,ω),(7)2.3螺旋测在上述各种隐式方法(4)和(7)式中,都需要对含有拉普拉斯算子的表示矩阵求逆.不失一般性,这里考虑形如亥姆霍兹算子D=αV2-Δ2(其中α=(1-4c2)/ic,α为常数)的表示矩阵求逆问题.这里以二维算子和4×2的网格为例,分析亥姆霍兹算子表示的矩阵H的特点.令网格上各点的波场组成列矩阵u=(u1,u2,u3,u4,u5,u6,u7,u8)T,对应点的速度值分别为V1,V2,V3,V4,V5,V6,V7,V8.此时在螺旋坐标h上,序列长度为8,在螺旋边界条件下,采用五点差分格式时该算子的表示矩阵为Η=1dx2×[d1-b00-1000-bd2-b00-1000-bd3-b00-1000-bd4-b00-1-100-bd5-b000-100-bd6-b000-100-bd7-b000-100-bd8],(8)其中di=αdx2V2i+m,m=2(1+b),b=dx2/dz2.矩阵H为具有Toeplitz结构的正定厄密矩阵,该矩阵的阶数即为螺旋坐标系上序列长度.在实际应用中,为模拟无限的物理空间,计算中常常采用吸收边界条件.同时采用螺旋边界条件和吸收边界条件,其计算效果与只采用吸收边界条件时的效果相同.3谱法的六安分解3.1矩阵的谱因式分解对形如(8)式的对称带状矩阵的求逆,可以通过LU分解来实现.谱法LU分解是一种快速分解方法.为讨论方便,可以令(8)式中的常数α为零,且取dx=dz,则亥姆霍兹算子的表示矩阵H的任意一列都化为r=(⋯,0,-1,0,⋯,0,-1,4,-1,0,⋯,0,-1,0,⋯)‚(9)r可以看成是某个函数a的谱函数,即r=a⨂a,这里符号“⨂”表示褶积,函数a称为r的最小自相关函数,它可由谱因式分解方法求得.(9)式中的r对应函数a=(1.791,-0.651,-0.044,-0.024,⋯,-0.044,-0.087,-0.200,-0.558),(10)如图1a所示.同样地,依次分解矩阵H的每一列ri,并把所得的函数ai按相应的位置排列,可得到一个下三角矩阵A,如图1b所示.这时便有Η=AA′,(11)其中A′为A的转置.谱法LU分解方法是近似的矩阵分解方法.它的不严格性在于其一对矩阵H中对应模型边界点的列的处理是不严格的.例如,在图1b中H的前4列不具有(9)式对称形式,但在矩阵A的对应列上仍采用(9)式谱分解结果.其二是当矩阵H的各列非零元素值存在差异时,(11)式的等式是不严格的.3.2谱法的吕分解有助于提高计算效率的原则3.2.1谱因式分解的速度求解函数r的最小相位自相关函数a的过程称为谱因式分解过程,它可由Kalmogoroff分解法或Wilson-Burg分解法实现.两种方法在计算量上基本一致,本文采用Wilson-Burg方法.谱因式分解的速度取决于(9)式中非零元素的分布位置.对H而言,由于速度模型在z方向上的网格点数nz决定了非零元素的分布位置,因而谱分解速度与nz相关.如图2a所示,谱因式分解时间随nz增大而近似线性地增大.具体地说,在速度模型网格为800×800时,按螺旋坐标系展开后,H为640000阶的方阵.对应的r是该方阵一个列元素,含有640000个元素.据图2a,分解列r约需0.27s(在主频为350M的微机上).3.2.2基于谱法健全的同步控制型用谱法LU分解一个矩阵,可分两步.第一步采用谱因式分解方法分解矩阵的一列.第二步把分解的结果按相应的顺序构造出矩阵LU分解的下三角矩阵.由于这第二步仅仅是一个排序的过程,其占用的时间可以忽略.因此,用谱法LU分解法分解一个矩阵的时间主要就是第一步所占用的时间.首先来考虑一种特殊情况,即矩阵H的各个列具有相同的非零值.如在均匀介质条件下,其螺旋坐标系下的表示矩阵H如图1b所示,矩阵的每一列具有(9)式中的非零值.因此,只需要针对(9)式做一次谱因式分解.这样,若以上述网格数为800×800的速度模型对应的640000阶矩阵H为例,完成LU分解的时间仅为做一次谱因式分解的时间,即0.27s.其次,考虑一般情况,如果采用谱因式分解方法分解矩阵H,由于各列具有不同的非零值,必须在第一步中依次采用3.2.1节中的谱因式分解方法分解矩阵的每一列.若仍以上述网格数为800×800的速度模型为例,分解每一列的时间是0.27s,由于矩阵H为含有640000个列元素,需要做640000次谱因式分解,这将花费170000s.即使在不考虑分解精度仅考虑时间的情况下,对一般矩阵采用谱法LU分解似乎也是一种非常笨的算法.因为在相同条件下,采用快速直解法只需0.93s就能完成该矩阵H的分解,无论H是否具有相同非零元素或不同非零元素.但是,在实际地震波成像应用中,矩阵H的各个列是具有相似的非零值,即非零值的变化范围是有限的,都是由该列的对应点的速度值决定的,如(8)式所示.即具有相同速度值的点对应的矩阵列的非零值相同.由于速度变化范围有限,因此非零值的变化范围也是很有限.因此,没有必要对H的所有列分别做谱因式分解.只需要对矩阵H中具有不同的非零值的列分别做谱因式分解.这个分解量是非常有限的,如对均匀介质,只需要做一次分解,对层速度模型介质,只需要对每层介质速度分别做一次分解.3.2.3谱法lu分解算法的建立在地震波场模拟和成像的计算中,往往面临的是对一系列具有相同结构和相似非零值的矩阵的分解.同样,由于介质速度值的变化范围是有限的,所有这些矩阵中的非零值的变化范围也是非常有限的.在采用谱法LU分解这些矩阵时,我们分两步来实现矩阵的分解.第一步就是针对所有可能出现的具有不同非零值的列(即针对所有可能出现的介质速度值)分别做谱因式分解,并把谱因式分解的值(即自相关函数的值)存在一个共查询的表中.这一步所占用的时间是:具有不同非零值的列的数目×谱因式分解一列所需要的时间.第二步就是构造每一个矩阵的下三角矩阵A.这个构造过程可以通过查询第一步中建立的表来实现.这一步几乎不需要占用时间.所需要的时间即为建立针对具有不同非零值的列对应的谱因式分解值的表的时间.因此,一旦建立起分解表,后面的矩阵分解计算速度就变得非常迅速了.图2b表示谱法LU分解所占用的时间与需要分解的矩阵数目的关系图.对应的速度模型网格数为800×800,谱法LU分解首先建立含有1000列自相关函数表.从图中可以看出,采用谱法LU分解矩阵几乎与需要分解的矩阵数目无关.因此该方法非常适用于巨大数目的矩阵的分解.直解法在处理大量矩阵的分解时只能是一个矩阵一个矩阵地分解.对分解次数不是很多时,它比谱因式分解快.但分解的次数很多时,它就不如谱法LU分解快了.图2b为快速直解法分解所占用的时间与需要分解的矩阵数目的关系图.可见直解法分解花费的时间与被分解矩阵的数目成正比.在被分解矩阵数目小于300时,直解法分解的速度比谱因式分解快.但对分解次数大于300次的计算,后者就比直解法快得多了.4谱法中lu分解的误差、误差分布和波场计算的影响4.1的速度网格若令(8)式中V=2000m/s,α=1,dx=dz=10m,则对应的亥姆霍兹算子表示H的所有列都具有相同的非零元素.对10×10的速度网格,矩阵H为一个100阶方阵,如图3a所示,由谱法LU分解得到的下三角矩阵如图3b.图3c和3d分别为误差矩阵H-AA′和它前20列的放大图.可见,谱法LU分解在前10列存在较大的误差,而这10列对应着速度网格左边界点.这说明谱法LU分解不能很好地处理矩阵中对应着模型网格左边界点的列的计算.在通常采用的吸收边界条件的计算中,这一误差不会影响整个波场的计算.4.2谱法lu分解误差若令(8)式中α=1,dx=dz=10m,速度模型为V=(2000+10×z)m/s如图4a所示,呈线性变化.网格为20×10,对应的矩阵H如图4b所示.谱法LU分解的误差矩阵H-AA′如图4c所示.图4d为误差矩阵前42列的放大图.谱法LU分解的误差不仅出现在模型边界点对应的列上,几乎所有的列上都存在误差(图4d中的白点和灰点分别表示误差的值为负和正).值得注意的是,所有白点出现在网格上部边界点对应的列上,在吸收边界条件下的这些误差是不影响波场计算的.但是误差还出现在网格中间点对应的列上(白点间的灰线),这些误差对波场计算的影响可是致命的.图4(e～g)分别为网格左边界上点A、网格中间点B和网格上部边界点C对应列的误差灰度图.由图4h可见,左边界点A对应列的误差随比率的增大变化不大,网格中间点B对应列的误差以及网格上部边界点C对应列的误差均随比率的增大迅速增大.4.3谱法lu分解误差同样,在20×10的二维网格,网格间距dx=dz=10m,位于深度z=110m层面之上的介质速度为2000m/s,层下速度为6000m/s,如图5a所示.由(8)式可知亥姆霍兹算子表示矩阵H中对应2000m/s速度点的列与对应6000m/s速度的列中的非零值的差异很大,如图5(b～d)分