基于特征波属性参数的立体层析速度反演方法研究

基于特征波属性参数的立体层析速度反演方法研究
洪瑛;韩文功;孙小东;李振春;李芳
【摘要】研究了利用法向入射点波(NIP波)的运动学属性,实现立体层析速度反演的方法.对这些运动学参数做层析反演,以得到用于深度域偏移成像的横向非均质平滑速度模型.多次迭代过程中,沿法向射线进行动力学射线追踪获得正演参数,且拾取的参数和正演参数误差达到最小,进而得到最佳速度模型.在正演模拟过程中,利用射线扰动理论计算出Frechet导数,使得目标函数梯度最小化.该方法拾取方便、易于实现,尤其适用于低信噪比的地区.将该方法应用于二维模型数据,收到了预期的效果.%In this paper investigation is conducted to implement stereo-tomography using kinematic attributes of eigen-wave which is so called normal incidence point (NIP) wave.These kinematic parameters are utilized in tomography to achieve smooth and lateral homogenous velocity model which can be applied to migration in depth domain.Multiple iteration are carried out to minimize the picked parameters and simulated ones through which optimal velocity model can be derived.During iteration, simulated parameters are obtained from dynamic ray tracing.Subsequently Frechet derivatives of tomography matrix are calculated via ray perturbation theory to minimize the objective function.This method is facilitated to operate with the advantage of convenient picking particularly applicable to low S/N ratio area.In this paper 2D model test shows expected and promising results using this approach.
【期刊名称】《物探化探计算技术》
【年(卷),期】2017(039)003
【总页数】8页(P359-366)
【关键词】特征波;立体层析反演;动力学射线追踪;射线扰动理论
【作者】洪瑛;韩文功;孙小东;李振春;李芳
【作者单位】中国石油大学(华东),青岛,266580;中石化石油工程地球物理有限公司,北京,100010;中国石油大学(华东),青岛,266580;中国石油大学(华东),青
岛,266580;中国石油大学(华东),青岛,266580
【正文语种】中文
【中图分类】P631.4
为实现地下地质结构的清晰成像和准确定位，建立合理准确的速度场模型是至关重要的。

目前有许多速度建模的方法，它们的区别在于对速度模型估计，速度模型更新量的计算和速度模型的描述方法不同(如块状、层状、平滑等)[1]。

流行的偏移速度建模方法通常基于CIP道集的剩余时差分析，其准则就是在速度模型正确的情况下，反射成像深度与偏移距无关[2-3]；深度聚焦分析[4-5]，其准则就是基于正确的速度模型做波场延拓，会在零时间零偏移距处聚焦。

对于低信噪比和低覆盖次数的数据，基于波动方程偏移的速度建模方法受速度更新不确定性的影响，而且深层的成像分辨率也不能令人满意。

重复多次的叠前偏移也需要考虑计算量的问题。

同时采集系统的不规则也会引起额外的偏移噪音，影响成像。

另外，常用的速度建模方法就是反射层析[6-8]，做反射层析时，基于当前速度模型和反射界面进行射线追踪。

使得正演的旅行时和相应的从叠前数据中拾取的旅行时之间的误差达到最小，以此来更新速度模型。

这种方法的缺点是它假设区域内层
位是连续的，并且需要从叠前数据中拾取大量旅行时信息。

尽管后来提出了一些反射层析自动化实现的方法，但仍需在叠前数据中识别和拾取反射旅行时，工作量非常大[9]，特别是在低信噪比的叠前数据中，仅拾取反射旅行时就很困难。

Billette等[10-11]提出了一种扩展型层析速度建模方法，该方法使用叠前道集中
的反射同相轴旅行时和局部斜率，来得到一平滑速度模型。

该方法的优势是只需考虑反射同相轴的局部形态，无需确定速度模型中的界面，相应地也就无需沿层拾取旅行时了。

由于全波形反演受诸多实际采集因素及算法本身的限制，多年来一直处于理论研究和部分海上数据应用阶段。

目前工业界主要采用的是基于Kirchhoff偏移或高斯束偏移共成像点道集的反射走时层析。

相对于偏移前的道集，共成像点道集是经过偏移后的道集，有着更高的品质，对速度变化更加敏感，并且消除了成像假象的干扰。

Kirchhoff偏移和高斯束偏移均属射线类偏移方法，具有速度快、适用于大数据体速度建模的优势。

随着计算机性能的提高，借助双程波波动方程的逆时偏移，可以生成更高质量的共成像点道集用于速度建模，也开始被工业界广泛采用，且在复杂构造地区应用方面显示出了一定的优势。

在最近几年里，共反射面元(CRS)技术的应用理念发生了变化。

以前仅把它作为叠加工具[12-18]，CRS叠加的副产品是CRS三参数即法向射线在地面的出射角α、两个特征波(法向入射点波NIP和法向波N)的波前曲率半径RNIP和RN[19]。

NIP波由地下反射层上一点源激发产生，初始曲率半径为0；N波由地下爆炸反射面产生，初始曲率半径与地层的局部曲率半径一致。

这些参数仅用来作为叠加过程的质量控制。

目前正逐步开发CRS参数的潜在用途[20]，与常规CMP叠加不同的是，CRS叠加的参数里面包含了地下反射层构造形态的信息(如倾角、曲率等)[21-22]，这使得利用CRS参数做立体层析速度反演成为可能[23-24]。

笔者的研究继承了共反射面元叠加的特点，通过特征波即法向入射点波(NIP波)的
运动学属性实现立体层析速度反演。

在Duveneck等[24]的NIP波属性走时层析的基础上，设计了较为复杂的、包含地下丘状反射界面的地质模型。

利用二维样条函数对其参数化，将最小平方QR分解算法应用于NIP波速度反演中，取得了较好的效果。

模型数据的测试表明，该方法对于复杂模型也具有很好的适用性，这在低信噪比地区资料的处理中将有一定的优势和应用前景。

NIP波前参数描述了共反射点CRP旅行时的二阶运动学响应。

NIP波立体层析速度反演的速度判别准则与深度聚焦分析是类似的，即在速度场正确的情况下，NIP 波波前沿法向射线反向传播后，会在零时间聚焦为一点[25-26]。

当然，对速度模型的描述方式也很重要。

相比于其他反演方法，利用NIP波做反演，是基于平滑连续的速度模型，从而无需考虑反射层位和沿层拾取分析点，易于实现。

1.1 反演问题的描述
利用NIP波曲率半径RNIP和近地表速度v0定义一个新的变量M：
其中:M代表了沿法向射线传播的NIP波旅行时的空间二阶导数。

使用M有两大好处：①对近地表速度v0的依赖可以取消，因为M是包含RNIP和v0的复合参数；②M易通过动力学射线追踪计算得到[27]。

反演的输入数据表示为：
(T,M,α,ξ)i，i=1,…,ndata
这里为了方便起见，T采用单程旅行时，ZO射线在地面坐标ξ处的出射角为α。

依据前述的准则，只要速度场正确，NIP波就会聚集。

这意味着反演计算时，固定α和ξ值，根据在地面拾取的NIP波前曲率半径或计算出的M，向地下反向传播所有的NIP波前，检测它们是否都在零时间处聚集(即T=0时RNIP=0),以此来判别速度模型是否正确[28-29]。

但这种思路难以应用于实际，因为所有的数据(T,M,α,ξ)都不可避免地受到噪音或测量误差的影响。

在噪声或测量误差不能得到有效处理的情况下，会导致反演的不稳定[30]。

采用相反的思路，如果利用动力学射线追踪方法从地下某一CRP点向上追踪正演
计算NIP波属性至地面，然后将其与拾取或观测的NIP波属性相比较，依据其匹
配程度来判断速度模型的正确性，就可以避免这种不稳定。

在这种思路下，真实的地下点位置和相应的法向射线出射角都是未知的，因此除了速度场外，这些参数也是反演的对象。

当模拟的参数组合(T,M,α,ξ)和拾取的参数组合(T,M,α,ξ)最接近时，即得到了最佳模型。

该方法类似于Billette等[10]提出的立体层析方法，既能反演出速度和地下绕射点位置，又能反演出射线角度参数。

在反演的二维模型中，每一个CRP点都有其地下坐标(x,z)和它对应的局部倾角
θ(或法向射线的起始方向)。

平滑的速度模型可以由二维 B样条函数表示：
其中：nx和nz分别是在水平和垂直方向的网格点数目；βj(x)和βk(z)是某一网格点的B样条函数。

需要对速度场的三阶导数做平滑，因此这里采用4阶B样条函
数(图2)。

要反演的模型可做如下表述：
CRP点坐标、局部反射倾角
(x,z,θ)i，i=1,…,ndata
和速度场
vjk，j=1,…,nx，k=1,…,nz
1.2 反演问题的求解
求解反演问题：寻找一个包含式(4)和式(5)中所有模型分量的模型向量m，使得拾取点的数据向量d和模拟的dmod=f(m)达到最佳匹配。

这里非线性算子f表示在给定模型中做动力学射线追踪。

如果采用L2模，反演问题中的目标函数为：
‖
其中:Δd(m)=d-f(m)；对角矩阵CD用来对不同的数据分量加权。

考虑到计算量的大小，这里不能采用用于解决f存在非线性特征问题的全局非线性最优化方法，而
是采用局部迭代的算法，将正演算子f进行局部线性化。

因此，目标函数S的最小值可以通过利用最小二乘方法求解线性问题得到。

反演过程中，基于初始估计模型，计算一系列的Δm，更新至满足要求。

在给定模型mn附近，第n步迭代的正演
算子可以局部估计为f(mn+Δm)≈f(mn)+FΔm，其中F为包含f在mn处的frechet导数的矩阵，其中，利用射线扰动理论可以计算出Frechet导数[7]。

计算目标函数梯度，令其为“0”，得到式(7)。

如果F存在逆矩阵，就可以得到Δm的最小二乘解，但是并不是所有的模型分量
都能够充分地受到数据地约束，这导致矩阵F通常是病态的，从而稳定的反演难
以实现。

因此需要一些附加信息，使问题规则化，其中一个附加的约束条件就是认为速度模型应该尽可能是平滑的、简单的。

平滑模型保证了可以利用动力学射线追踪沿着法向射线，计算稳定可靠的NIP波属性，即波前曲率半径和角度参数，进
而描述CRP运动学响应[29]。

作为对速度场而不是速度更新量的附加约束条件，
平滑约束准则(二阶空间导数最小化)可以表示为式(8)的矩阵方程：
利用最小二乘方法求解方程(8)即可得到模型更新量Δm，使得目标函数最小化：
其中:m(v)是m中包含速度系数的部分；ε因子用于调节目标函数中的模型平滑约束条件和数据残差量约束条件之间的相对权重关系。

为了保证反演的稳定，需要依据具体数据对因子大小做合适地选择。

为了得到式(9)的解，笔者应用了LSQR迭代算法。

该算法可以高效地求解线性问题，而且不需要显式计算逆矩阵。

该算法可以利用稀疏矩阵的特点，在地震的层析反演领域被广泛地应用。

1.3 反演问题的流程
在应用反演算法之前，需要首先得到地面不同位置不同时刻的NIP波属性参数，
即波前曲率半径和出射角度。

这些参数需要从叠前地震数据中分析得到，一个办法就是利用共反射面元叠加的“副产品”参数剖面。

从CRS叠加得出的ZO模拟剖
面及相干剖面上拾取数据点，然后从相应的参数剖面上提取NIP波属性参数。

数
据点的拾取是相互独立的，无须沿层或连续的同相轴，并应尽可能地让数据点均匀分布于整个反演区域。

地震波的传播是一个渐变的过程，所以可以认为NIP波属性参数沿反射层平滑变化，且在小于子波长度的时窗内是相对稳定的。

提取过程中要首先去除野值，并进行中值滤波、平滑消除非物理波动。

参数提取时，平滑算子在时间方向上不应超过子波长度，防止串相位，在空间方向不应超过第一菲涅尔带半径，和CRS叠加孔径相一致并去除相干值低于预定门槛
值的点。

理论上，利用CRS叠加过程中得到的相干剖面，就可以自动进行数据点地拾取。

但是，实际数据的自动拾取大多会受到多次波的干扰。

考虑到拾取数据点的数量是很有限的，通常采用人机交互的方法。

在CRS叠加剖面上拾取了数据分析点后，
其对应的α和RNIP/M参数就可以从对应的CRS属性剖面上自动提取出来。

起始速度模型的建立方法如下：设定B样条在水平方向和垂直方向上的节点序列，将初始速度值赋给速度系数。

例如可以定义一个存在垂向梯度的一维速度模型，或者利用DIX公式反演出的水平层状速度模型作为初始模型。

对每一个拾取的数据点(T,M,α,ξ)，首先基于初始模型向地下进行射线追踪，就得
到相应的CRP点的起始坐标和局部倾角(x,z,θ)，然后反过来从这些CRP点出发，向上做动力学射线追踪直到射线到达地面。

在这过程中，计算Frechet导数，起
始阶段，Δd仅包含观测数据和模拟数据的残差，确定更新后的模型
mn+1=mn+λΔm，其中0<λ≤1。

基于更新后的模型再做动力学射线追踪得到新的模拟数据，如果目标函数增大，就减小λ值，重新计算目标函数；否则，继续
更新模型，计算frechet导数，继续下一次迭代计算。

在不断迭代过程中，随着模型越来越精确，ε可以逐渐减小，从而可以反演出速度
场更多的细节信息。

如果迭代次数达到预定值、数据残差小于预定精度或目标函数达到最小值(即减小λ不会导致目标函数值的减小)就终止迭代，输出最终反演结果。

需要说明的是，笔者仅对模型数据做了测试，以验证该方法的正确性，所采用的NIP波属性参数是事先基于动力学射线追踪正演得到，以此作为反演算法的输入。

在处理实际资料时，NIP波属性应该从实际的叠前数据中基于CRS叠加分析得到。

设计一平滑速度模型如图2(a)所示，模型的维数为323×240。

根据前述的NIP波立体层析反演的算法流程，从图2(a)中设定的起始点进行射线追踪(图2(b))，得到射线达到地面处的时间T、属性参数M、出射角度参数α(图2(c)、2(d)、2(e))。

将正演得到的这些参数T 、M 、α作为NIP波速度反演的输入，构建初始模型，设定控制点处的速度值，做B-样条平滑构建的初始速度模型(图3(a))。

NIP波速
度反演中，第1次迭代计算射线路径图、T、M、α(图3(b)、3(c)、3(d)、3(e))。

基于初始速度模型图3(a)，利用正演得到的参数 T 、M 、α与标准参数T 、M 、α之间的差异计算出控制节点处的速度更新量。

根据模型控制点，做B-样条平滑
构建的平滑速度模型，作为第一次迭代得出的速度模型(图4(a))。

基于该模型做第2次迭代计算射线路径图、T、M、α(图4(b)、4(c)、4(d)、4(e))。

经过23次迭代，最终反演出的速度模型如图5(a)所示，图5(b)为其与标准速度模型之差。

由模型试算的结果可以看到，笔者利用NIP属性参数进行的立体层析反演方法，
可以很好地重建速度场，得到的反演结果误差主要集中在模型的底层且并不是很大。

由图5(b)所示，3 500 m处的反演结果误差不大，但随着深度的增加，由于越往
深层NIP波射线越稀疏，相对模型的约束效果降低，导致结果误差逐渐变大。

NIP波立体层析速度反演方法存在一些优势：
1) 与传统算法相比，本算法可以高效拾取数据分析点，能够克服基于反射波的层
析反演算法存在的效率低等“瓶颈”问题。

本算法不是在原始叠前数据中拾取数据
分析点的，而是在提高信噪比的CRS叠加剖面上进行拾取的。

无需像其他层析方法一样在叠前数据中逐个偏移距拾取旅行时，NIP波立体层析速度反演方法利用在某一数据分析点拾取的NIP属性，即可算出CRP点对应的不同偏移距的旅行时，解决了前者计算工作量过大的问题。

而且本算法也无需沿层连续拾取观测点，因为基于平滑速度模型，并不一定要求拾取的反射、绕射点位于速度或层位分界面。

2) 共反射面元(CRS)叠加方法，是将从叠前数据中以纯数据驱动的方式得到的出射角度和NIP波波前曲率作为反演的输入，能够在不影响分辨率的情况下，在菲涅尔带范围内进行地震数据叠加，极大地提高了地震资料的信噪比。

NIP波属性的拾取是在增强信噪比后的CRS剖面中进行的，因此本算法在处理低信噪比地区资料时具有很强的适应能力。

3) 零偏移距斜率层析/立体层析在射线旅行时的基础上，考虑了地层倾角/射线出射角度，从而扩展了正反演的模型参数空间。

更多的参数数据体用于反演可以提高反演的可靠性，但如何从低品质的地震数据中获得高质量的参数数据体用于层析算法的输入，在实用中是一个挑战。

NIP波立体层析速度反演可以借助CRS叠加获取高质量CRS参数的优势，保证了速度反演的效果，同时由于考虑了NIP波的曲率(二阶项)，相对于仅考虑射线出射角度(一阶项)的斜率层析/立体层析，对速度反演结果有更好的约束作用。

该方法也存在一些不足：
1) 基于NIP波属性的立体层析反演方法，在应用于地下速度场复杂区域存在较大误差，这是因为该方法是沿法向射线计算NIP波属性的。

此外，为了确保双曲近似是合理的，利用CRS叠加求NIP波属性需要确定剩余叠加的CMP道集数量以及偏移距允许的最大值，但也不能把叠加孔径设置的太小，会降低NIP波属性分析的准确性。

2) CRS叠加在计算时距曲线公式中可将地表高程信息考虑进来，直接处理未做静
校正的数据，进而解决了起伏地表对采集数据的影响。

但由于实际资料中常存在近地表强横向变速的问题以及双曲近似和球形波前假设，基于NIP波属性的立体层
析反演方法，无法较好地处理起伏地表地区的速度反演问题。

笔者讨论了利用NIP波属性进行立体层析速度反演的方法。

需要反演的模型参数
有平滑速度场、反射点位置及起始角度。

反演无须考虑层位，沿层连续拾取，只需随机按照相干值高的点独立拾取即可，容易实现。

NIP波属性的立体层析速度反演，在处理低信噪比资料时比其他反演方法更具优势。

笔者将NIP波立体层析反演方
法应用于一横向变速介质的模型中，得到了预期的结果，从而验证了该方法的正确性。

需要指出，该方法基于高质量的叠后CRS参数。

基于NIP波属性的立体层析反演方法，能快速稳健地得到较为合理的速度模型，可以作为后续精细的波动方程速度建模和全波形反演的高质量的初始模型，这对提高实际地震资料处理的效果有着借鉴作用。

【相关文献】
[1] DIX, C. H.. Seismic velocities from surface measurements[J].Geophysics, 1955(20): 68-86.
[2] SPRINGER-VERLAG.DEREGOWSKI,mon-offset migrations and velocity analysis[J].First Break,1990,8(6):224-234.
[3] LIU,Z..An analytical approach to migration velocity
analysis[J].Geophysics,1997(62):1238-1249.
[4] JEAN-PIERRE FAYE ,JEAN-PAUL JEANNOT.Prestack migration velocities from focusing depth analysis[C].SEG Technical Program Expanded Abstracts,1986:438-440.
[5] MACKAY, S., Abma, R.. Imaging and velocity estimation with depth-focusing analysis[J]. Geophysics, 1992(57): 1608-1622.
[6] BISHOP, T., BUBE, K., CUTLER, R., et al. Tomographic determination of velocity and depth in laterally varying media[J]. Geophysics, 1985(50): 903-923.
[7] FARRA,V.,Madariaga,R..Seismic waveform modeling in heterogeneous media by ray perturbation theory[J].Journal of Geophysical Research,1987,92(B3):2697-2712.
[8] STORK, C.. Reflection tomography in the postmigrated domain[J]. Geophysics,
1992(57):680-692.
[9] BILOTI, R., SANTOS,L.T., TYGEL, M.. Multiparametric traveltime inversion[J]. Studia Geophysica et Geodaetica, 2002(46): 177-192.
[10]BILLETTE, F., LAMBAR′ e, G.. Velocity macro-model estimation from seismic reflection data by stereotomography[J]. Geophysical Journal of International, 1998(135): 671-690. [11]BILLETTE, F., Le B′ egat, S., PODVIN, P., et al. Practical aspects of 2D stereotomography[J]. Geophysics, 2003(68): 1008-1021.
[12]LI Z C, SUN X D, ZHAO J Z, et al. Common reflection surface stack based on complex top surface topography[C]. BGP Technical Seminar, 2004.
[13]SUN X D, LI Z C, TENG H H. Common reflection surface stack using dip decomposition for rugged surface topography[J]. Applied Geophysics, 2007, 4(1): 45-50 .
[14]李振春, 孙小东, 刘洪. 复杂地表条件下共反射面元(CRS)叠加方法研究[J]. 地球物理学报, 2006, 49(6): 1794-1801. LI Z C, SUN X D, LIU H. Common reflection surface stack for rugged surface topography[J]. Chinese Journal of Geophysics, 2006, 49(6): 1794-1801. (In Chinese) [15]孙小东, 李振春, 丁仁伟. 基于最优化方法的共反射面元叠加[J]. 中国石油大学学报, 2008, 32(5): 46-49. SUN X D,LI Z C, DING R W. Common reflection surface stack based on optimization method[J]. Journal of China University of Petroleum, 2008, 32(5): 46-49. (In Chinese) [16]王华忠, 杨锴, 马在田. 共反射面元叠加的应用理论—从共反射点到共反射面元[J]. 地球物理学报, 2004, 47(1): 137-142. WANG H Z, YANG K, MA Z T. An applied theory on common reflection surface stack-From common reflection point to common reflection surface[J]. Chinese Journal of Geophysics. 2004, 47(1): 137-142. (In Chinese)
[17]杨锴, 王华忠, 马在田. 共反射面元叠加的应用实践[J]. 地球物理学报, 2004, 47(2): 327-331 YANG K, WANG H Z, MA Z T. The practice on common reflection surface stack[J]. Chinese Journal of Geophysics, 2004, 47(2): 327-331. (In Chinese)
[18]谭未一, 杨长春, 李瑞忠, 等. 共反射面元叠加的实现途径及流程[J]. 地球物理学进展, 2004,
19(2): 325-330. TAN W Y, YANG C C, LI R Z, et al. Common reflection surface stack algorithm and processing flow[J]. Progress in Geophysics, 2004, 19(2): 325-330. (In Chinese)
[19]孙小东, 李振春, 滕厚华,等. 共反射面元技术及其在偏移成像中的应用[J]. 石油物探, 2008, 47(5): 465-469. SU X D, LI Z C, TENG H H, et al. Common reflection surface stacking and its application in migration imaging[J]. Chinese Journal of Geophysics, 2008, 47(5): 465-469. (In Chinese)
[20]SUN X D, LI Z C. Application of common reflection surface technology from stack to Prestack time migration[C]. 73rd EAGE Conference, 2011.
[21]J¨AGER, R., MANN, J., H¨ OCHT, G., et mon reflection-surface stack: Image and attributes[J]. Geophysics, 2001, 66, 97-109.
[22]MANN, J., JA¨GER, R., MU¨ LLER, et al. Common-reflection surface stack—A real data example[J].Journal of Applied Geophysics,1999,42:3, 4,283-300.
[23]BERKOVITCH, A., GELCHINSKY, B.. Inversion of common reflecting element (CRE) data (migration combined with interval velocity determination)[c]. 59th Annual International Meeting, Society of Exploration Geophysicists, Expanded Abstracts, 1989.
[24]DUVENECK, E., HUBRAL, P.Tomographic velocity model inversion using kinematic wavefield attributes[C]. 72nd Annual International Meeting, Society of Exploration Geophysicists, Expanded Abstracts, 2002.
[25]MACKAY,S.,ABMA,R..Imaging and velocity estimation with depth-focusing
analysis[J].Geophysics,1992,57:1608-1622.
[26]GELCHINSKY,B.,BERKOVITCH,A.,KEYDAR,S.Multifocusing homeomorphic imaging, part 1—Basic concepts and formulae[J].Journal of Applied Geophysics,1999,42(3):4,229-242.
[27]李芳. 共反射面道集速度分析[M].东营:中国石油大学(华东),2009. LI F. Velocity model estimation based on common reflection surface gather[M]. Dongying: China University of Petroleum (East China), 2009. (In Chinese)
[28]H¨OCHT,G., DE BAZELAIRE, E., MAJER, P., et al. Hyperbolae and curvatures[J]. Journal of Applied Geophysics, 1999, 42(3):4, 261-281.
[29]SHAH, P. M.. Use of wavefront curvature to relate seismic data with subsurface parameters[J]. Geophysics,1973(38): 812-825.
[30]ERIC DUVENECK.Velocity model estimation with data-derived wavefront
attributes[J].Geophysics, 2004,69(1):265-274.
[31]孙小东，李振春. 特征波立体层析速度反演方法及应用研究[C]. 中国地球物理学会第二十九届年会. 2013. SUN X D, LI Z C.Eigen-uave stereo-tomography velocity inversion and its application[C]. 29th Chinese Geophysical Society Seminar. 2013. (In Chinese)。