有限频层析成像方法

有限频体波层析成像简介徐小兵赵亮传统射线走时层析成像基于高频射线理论。

在地球内部的速度异常变化不大，且异常体的尺寸要远大于地震波的特征波长的情况下，可以根据费马原理确定地震波沿走时最短的路径，此时可以认为到时差主要取决于射线路径上的速度结构。

然而，实际的地震波具有有限频率，当地震波经过尺寸较小或相近（相对地震波波长）的速度异常体时会发生绕射现象。

随着传播距离的增加，地震波波前面由于非均匀介质导致的超前或滞后特征能逐渐愈合甚至消失（Nolet and Dahlen, 2000；Hung et al., 2001）。

地震波的这种有限频率特性表明地震到时差并不主要取决于射线路径上的速度结构。

针对上述问题，Dahlen 和Hung 等人推导了地震到时在有限频地震波传播波场下的灵敏场函数（或根据其几何形状被简称为“香蕉—甜面圈”理论）（Dahlen et al.,2000；Hung et al., 2000），进而发展了有限频层析成像方法。

相对传统射线走时层析成像，该方法由于更真实地反映地震波的传播特征，因此能够更好地重建地球深部结构速度图像。

一、有限频层析成像方法原理地震波的有限频率特性决定了不同频率的地震波的走时对射线周围三维速度结构具有不同的敏感程度。

Dahlen 等（Dahlen et al., 2000）结合体波传播理论和线性Born单次散射近似来模拟通过地球内部任意一点的散射波到达台站的波形扰动，利用不同频率的体波与未经散射的直达体波之间的相互干涉结果，得出有限频率走时变化与中心几何射线附近第一菲涅尔带内的三维速度扰动之间的关系为：()()()x d x c x c w x K T ⎰⎰⎰=3/,δδ (1)式中，T δ是经过波形互相关处理得到的散射波与直达波的走时差；()()x c x c /δ为空间点x 处的体波速度扰动量；()w x K ,是该点三维灵敏度核，代表了某个散射点x 处的速度异常对走时变化的贡献。

从上式可以看出，敏感场K 不仅是空间位置x 的函数，同时还是频率w 的函数，即不同频率的地震波对不同空间范围的速度结构进行了采样。

在射线理论下地震波走时表示为慢度沿射线路径的线积分，而在有限频率灵敏核理论中地震波走时可以更精确地表示为灵敏度核函数在射线周围一定空间内的体积分，这是有限频率理论与射线理论最本质的差别。

因此，确定在一定速度模型下的有限频率灵敏核函数是进行有限频层析成像的关键所在。

二、有限频率地震体波灵敏度核函数及其性质获得有限频率地震体波灵敏度核函数是开展有限频层析成像的基础。

图1. 有限频体波走时在球对称介质下的几何路径示意图（Hung et al 2004）()()()()[]()∑∑⎰⎰∞∙∞∙---+--+⨯⎪⎪⎭⎫ ⎝⎛ℜℜℜΩ--='''022023,''',2''''''sin '''121rays rays r c N N N c d m d M M M T T T m c c K ωωωωπωωωπρρ……（2） 公式（2）表明有限频敏感度核函数是对时间窗内所有单点散射波进行积分求和的形式。

由于对所有满足条件的单点散射波进行射线追踪是一个非常耗时的过程，因此在实际成像过程很难实现。

针对此情况，Dahlen 等（2000）在计算敏感场函数时引入旁轴近似，即认为走时残差只对第一菲涅尔带内的速度扰动比较敏感，从而来有效的减少射线追踪的运算量。

采用的旁轴近似的公式为：()()()⎰⎰∞∞∆⎥⎦⎤⎢⎣⎡ℜℜℜ-=022023sin '''21ωωωωωωωπd S d T S c c K syn syn r (3)图2. 公式2得出的直达P波灵敏度核函数（Hung et al 2000）图3. 根据旁轴近似（公式3）给出的直达P波灵敏度核函数（Hung et al 2000）从图2和图3的对比中我们可以看出通过旁轴近似由公式3得出的灵敏场可以很好地对由公式2得出的理论敏感场进行近似。

图5. 灵敏度核函数与频率的关系（Hung et al., 2004）由高、中、低三个不同频带内的灵敏度核函数（图5）的形态可以看出，地震波的频段越高其灵敏度核越窄且强度越大，最接近于传统的射线；而频带范围越低则灵敏度核越宽，这时的灵敏度核强度也会降低。

理论分析表明，对宽频带波形分频带处理，不同频带内的分析结果可以分别约束中心射线路径周围不同范围处的速度结构。

三、有限频理论与射线理论的比较图6. 冰岛下方射线理论与有限频方法的灵敏场函数分布对比。

（a）射线方法的灵敏场；（b）有限频方法的灵敏场。

图件来自Hung et al., 2004。

图6展示了射线方法和有限频率方法的灵敏场分布：相对射线方法，有限频率方法对空间的采样更加密集、平滑，且空间范围更大，因此能够充分利用宽频带资料中的信息，减少数据不均匀覆盖的影响，进而能更加充分对速度异常体进行重建。

如同很多对比实验所示，有限频层析成像能够提高反演结果分辨率。

华北地区的实际观测结果也表明，有限频方法获得的异常体分布尽管大体与射线方法的结果一致，但重建得到的异常体扰动幅度更大、分辨率更高（图7）。

基于层析成像反演理论，射线走时层析成像和有限频走时层析成像最终都是求解模型参数与观测数据所构成的方程组，两者的不同之处在于正演计算模型参数和观测数据的（敏感场泛函）之间的不同。

(a)(b)图7. 采用射线方法和有限频率方法得到的层析成像结果比较。

(a)采用射线方法得到的150公里深度切片，左图为P波结果，右图为S波结果；(b)有限频方法结果。

四、反演过程应用网格节点对模型空间参数化后，可将有限频率地震层析成像的反演问题写成下面的矩阵方程：Gm d =式中，d 为观测走时向量；G 为灵敏度核函数矩阵，该矩阵第i 行第j 列的元素G ij 是第i 个观测走时的灵敏度核函数在三维网格节点j 处一定空间范围内的体积分；m 为要求求解的速度扰动向量。

反演问题转化为求解线性方程组b Ax =的问题，即是求解满足最小二乘条件2min b Ax -的解。

A 为n m ⨯阶矩阵；x 为n 维矢量，称为模型空间矢量；b 为m 维矢量，称为数据空间矢量。

一般数据空间矢量的维度m 大于模型空间矢量的维度n 。

A 为大型稀疏矩阵，如果采用传统的求解广义逆的方法来求解，计算机的内存需求很大，计算效率很低。

所以反演过程采用LSQR 方法。

LSQR 算法运用了三个主要的技术：(1) Lanczos 迭代： 这个迭代方法产生模型空间的正交基矢量()k v ,1γ以及数据空间的正交基矢量()k u ,1γ。

令数据空间初始矢量为b u =11β，则相应的模型空间初始向量为111u A v T =α，利用Lanczos 迭代构建正交基，把模型按照正交基展开，得到()k k ⨯+1阶的下三角矩形阵()k B ；如果令模型空间初始向量为b A V T =11θ，则相应的数据空间初始向量为111v A p =ρ，同样利用Lanczos 迭代构建正交基，把模型按照正交基展开，得到k k ⨯阶上双对角方阵()k R 。

由于所选的正交基是相同的，所以这两个双对角方阵一定存在特定的关系。

(2) Givens 变换运用一系列非常有规律的Givens 平面旋转正交变换实现正交变换，建立这两个双对角矩阵之间的联系。

由()k B →()k R 的过程即是通过Givens 变换对()k B 的QR 分解。

这个过程实现了最小二乘和QR 分解的完美结合（这也是被称为LSQR 的原因）。

(3) 模型迭代通过迭代运算直接计算模型矢量()k x ~，逐步逼近最小二乘解。

LSQR 算法程序设计的框架流程如下：1）初始化()()().,,0~,,,111101*********αρβφααββ========x v w u A u A v d d u T T2）循环计算 for k=1,2,3,...重复(3)-(6).3）Lanczos 算法，双对角化k k k T k k kk k k k v u A v u v A u 111111++++++-=-=βααβ4）正交变换kk k k k k k k k k k k k k k kk k k k k s c c s s c φφαρφφαθρβρρβρρ=-=---------====+=+++++++1111112121+k ρ与1+k φ是下一次迭代的起始值。

5）更新x ~和w ()()()()kk k k k k k k k k w v w w x x ρθρφ1111~~++++-=+=6）检验是否收敛如果达到收敛条件，停止迭代。

LSQR 算法的优点是，计算机的内存要求低，易于实现并行算法，迭代收敛快，求解精度高。

缺点是给不出显示的广义逆，也就给不出显示的模型分辨矩阵和协方差矩阵。

如下是层析成像反演流程图：NY开始 选取研究区域 拾取相对走时残差 建立初始速度模型 选用LSQR 反演方法调整上一次模型拟合观测和理论相对走时残差 得到最终的速度扰模型选相对走时残差是否小于设定阈值 结束参考文献1.Dahlen, F. A., S.H., Hung, G. Nolet, Frechet kernels for finite-frequencytraveltimes -I. Theory[J].Geophys J Int,141:157-174 , 20002.Hung, S.H., F. A., Dahlen, G. Nolet, Fréchet kernels for finite-frequencytraveltimes -Ⅱ.Example, Geophys J Int ,141:175-203, 2000.3.Hung, S. H., Y. Shen, L.Y., Chiao, Imaging seismic velocity structure beneath theIceland hot spot : a finite frequency approach. J Geophys Res,109:B08305,doi:10.1029/2003JB002889 , 20044.Zhao, L., R. M. Allen, T. Zheng, R. Zhu., High-resolution body-wave tomographymodels of the upper mantle beneath eastern China and the adjacent areas,Geochem. Geophys. Geosyst., 13, Q06007, doi: 10.1029/2012GC004119, 2012. 5.Zhao, L., R.M. Allen, T.Y. Zheng, S. Hung (2009). Reactivation of an Archeancraton: Constraints from P- and S-wave tomography in North China, Geophy. Res.Lett., 36, L17306, doi:10.1029/2009GL039781.。