地球物理学报  2014, Vol. 57 Issue (5): 1655-1663   PDF    
引用本文
张固澜, 贺振华, 王熙明, 张建军, 张彦斌, 李飞. 2014. 地震波频散效应与反Q滤波相位补偿. 地球物理学报, 57(5): 1655-1663, doi: 10.6038/cjg20140528   
ZHANG Gu-Lan, HE Zhen-Hua, WANG Xi-Ming, ZHANG Jian-Jun, ZHANG Yan-Bin, LI Fei. 2014. Seismic wave dispersion effects and inverse Q-filter phase compensation. Chinese J. Geophys. (in Chinese), 57(5): 1655-1663, doi: 10.6038/cjg20140528   
地震波频散效应与反Q滤波相位补偿
张固澜1,2, 贺振华1, 王熙明2, 张建军2, 张彦斌2, 李飞2    
1. 成都理工大学油气藏地质与开发工程国家重点实验室, 成都 610059;
2. 中国石油集团东方地球物理公司新兴物探开发处, 河北 涿州 072751
投稿时间:2013-11-7    最后修改时间: 2014-4-11 .
基金项目：国家自然科学基金项目（41174114）资助.
作者简介：张固澜，男，1982年生，从事井中地震资料处理方法研究及井驱动的地面地震高分辨率研究.E-mail：zhanggulan@cnpc.com.cn
摘要：本文以Futterman提出的地震波振幅衰减和相速度频散表达式为基础，从井震匹配的角度出发，推导了不进行反Q滤波、仅反Q滤波振幅补偿、仅反Q滤波相位补偿，以及反Q滤波相位与振幅同时补偿四种情况下，地震波速度频散与相速度及地震记录振幅谱间的表达式，并从理论上说明了反Q滤波相位补偿的必要性.通过相同观测系统，炸药震源和可控震源分别激发采集的零偏移距VSP资料实例，验证了本文所推导的速度频散表达式的合理性；通过井震标定实例，进一步说明了反Q滤波相位补偿，可有效消除地震子波速度频散，提升地面地震资料与零偏移距VSP走廊叠加剖面的匹配度，最终提高地震资料成果的可信度.
关键词： 相速度     频散     群速度     井震标定     反Q滤波     相位补偿     零偏移距VSP    
Seismic wave dispersion effects and inverse Q-filter phase compensation
ZHANG Gu-Lan1,2, HE Zhen-Hua1, WANG Xi-Ming2, ZHANG Jian-Jun2, ZHANG Yan-Bin2, LI Fei2    
1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology, Chengdu 610059, China;
2. BGP, CNPC, Zhuozhou Hebei 072751, China
Abstract: Basing on the expression of amplitude attenuation and phase velocity dispersion proposed by Futterman, and considering the well-seismic matching, we derived the expression of seismic wave velocity dispersion between phase velocity and amplitude spectrum for the following four cases: without inverse Q-filter; with only inverse Q-filter amplitude compensation, with only inverse Q-filter phase compensation, with inverse Q-filter phase and amplitude compensation, and explained the necessity of inverse Q-filter phase compensation in theory. Through the zero-offset VSP data example of dynamite and vibroseis source with the same acquisition geometry, we verified the velocity dispersion expression derived in this paper. Through the example of well-seismic calibration, we further illustrated that inverse Q-filter phase compensation can effectively eliminate seismic wave velocity dispersion, enhance well-seismic matching, finally improve the reliability of surface seismic data.
Key words: Phase velocity     Dispersion     Group velocity     Well-seismic calibration     Inverse Q-filter     Phase compensation     Zero-offset VSP    

1 引言

地震波传播过程中，要经受地层Q吸收造成的能量衰减(Hamilton，1972;Toksoz et al., 1979)和速度频散(Sams et al., 1997;Spencer et al., 1982). 能量衰减使得地震记录的主频降低，有效频带变窄，分辨率降低.速度频散使不同频率的地震波具有不同的传播速度，造成相位畸变.特别是当地层含有油气时，能量衰减和速度频散现象更加明显.

为描述地层的Q吸收效应，Futterman(1962)提出了能量衰减和相速度频散表达式；Kjartansson(1979)根据不同的假设条件对Q吸收效应进行了数学描述；为消除地层的Q吸收效应，有关学者利 用VSP资料(Hauge，1981; Stainsby and Worthington， 1985; Badri and Mooney， 1987;Tonn，1991; Xu C and R Stewart，2006; Gao and Yang， 2007; Gao et al., 2008; Blias，2012)和地面地震资料(Yan and Liu， 2009; Wang，2011;Zhao et al., 2013)反演地层品质因子Q，并用于反Q滤波(Hargreaves and Calvert 1991; Wang， 2002，2003，2006;Yao et al., 2003; Liu et al., 2013;Chen et al., 2014)方法，对地震波的振幅和相位进行补偿.反Q滤波振幅补偿提高地震分辨率的能力，受到业界广泛认可，得到了较好的应用；而反Q滤波相位补偿(Bano，1996)用于校正速度频散方面，实际资料处理效果鲜有发表.

地震资料采集中，受激发环境限制，经常采用多种类型震源(如气枪震源，炸药震源、可控震源、重锤等)激发，而不同类型震源激发的地震子波具有不同的频带和相位；因此即便激发点相同，接收点相同，地震记录的有效频带也会不同，速度频散也不同，这给地震资料匹配处理(如连片处理，时移地震、井震标定)带来很多难题.随着目标体勘探精度的提高，我们不得不重视速度频散和反Q滤波相位补偿的研究.

因此，本文以Futterman提出的地震波振幅和相速度频散关系式为基础，从井震联合标定的角度出发，结合与反Q滤波处理相关的四种情况，详细地分析了VSP资料与地面地震资料之间的速度频散关系，并从理论上说明了反Q滤波相位补偿的必要性；通过相同观测系统，两种不同类型震源采集的零偏移距VSP资料实例，验证速度频散关系；并通过井震联合标定实例，展示了反Q滤波相位补偿效果，验证了反Q滤波相位补偿的必要性.

2 地震波非弹性衰减与补偿 2.1 地震波的非弹性衰减

为描述地层吸收效应，Futterman提出了振幅衰减和频散方程；该方程基于Q与频率无关的假设和一维双程传播波动方程:

其中，U(z，ω)为频率ω的波在深度z的振幅谱，k(ω)为z方向波数，且.考虑距离增量与走时增量关系：

导出如下振幅和相速度频散关系式：

式中Q是品质因子；γ= ，T是传播时间，U(0，ω)是震源子波振幅谱，U(T，ω)是传播时间T处子波振幅谱，v(ω)是频率为ω的地震波的相速度，v(ωc)是参考频率为ωc的地震波的相速度.

将v(ω)/v(ωc)对ω求一阶偏导数，可得相速度频散的变化率：

由此可知，随着ω的增大，相速度的频散变化率会相应地减小，当ω大于等于ωc时，相速度频散变化率趋于零.

利用中心频率为50 Hz的带通子波，得到图 1中FFID=1的地震道；利用方程(2)，按Q=138进行Q吸收衰减正演，得到图 1中FFID=2的地震道：随传播时间增大，地震记录能量衰减，主频降低；由于地震波相速度的频散效应，子波由原来的零相位变成混合相位.

图 1 反Q滤波前后效果对比 Fig. 1 The result of synthetic tarce before and after inverse Q-filter
2.2 反Q滤波

实际资料处理中，根据不同的目的需求，反Q滤波往往只补偿振幅或相位，或同时补偿相位和振幅.由方程(2)，可得这三种情况下，反Q滤波的表达式.

如只进行反Q滤波振幅补偿，则反Q滤波的表达式为

如只进行反Q滤波相位补偿，则反Q滤波的表达式为

如同时进行反Q滤波相位补偿和振幅补偿，则反Q滤波的表达式为

利用Q=138，对图 1中FFID=2的地震道只进行相位补偿，得到图 1中FFID=4的地震道：反Q滤波相位补偿后，各时刻的记录都恢复成了零相位子波，消除了频散效应.

利用Q=138，对图 2中FFID=2的地震道同时进行反Q滤波相位补偿和振幅补偿，得到图 1中FFID=5的地震道：补偿后的记录与原来的合成记录完全匹配，不仅仅消除了频散，恢复了能量，还提高了分辨率.

图 2(a)炸药震源激发时下行子波振幅谱；(b)可控震源激发时下行子波振幅谱 Fig. 2(a)Amplitude spectrum of down-going wavelet with dynamite source; (b)Amplitude spectrum of down-going wavelet with vibroseis
3 地震波频散效应理论分析 3.1 理论模型

假设模型为一个厚度为H、品质因子为Q的均匀各向同性介质；地面地震记录与零偏移距VSP记 录具有相同的震源子波，且震源子波振幅谱为U(0，ω)； 地震波的参考频率为ωc，且其传播速度为v(ωc)；地面地震记录炮检距xi处，参考频率的地震波传播时间为ti；零偏移距VSP检波器埋置深度为H，且参考频率的地震波传播时间为tz，则有

3.2 理论分析

利用3.1节的理论模型，结合反Q滤波处理四种情况，定量分析地面地震资料与零偏移距VSP资料之间的速度频散关系.

(1)不进行反Q滤波补偿

不进行反Q滤波补偿，相当于利用方程(2)进行地震波正演.此时，零偏移距VSP资料中，地震子波的振幅谱可表示为

地面地震资料中，炮检距xi处地震子波的振幅谱可表示为

由方程(10)和(11)可得，随传播时间增大，地震波能量衰减增大，主频向低频方向移动，低频信号的相对能量逐渐增强，所占比重越来越大；传播时间无穷大时，地震波群速度必然等于某一特定频率(该频率远小于激发子波的主频)的相速度.再据方程(4)可得，对同一目的层而言，炮检距越大，速度频散越大，群速度越小.

因此，地震波的传播时间不仅和速度频散有关，还和各种频率的相对能量有关；地震波的传播时间可表示为各频率成分的地震波传播时间的加权，加权系数为各频率成分的能量在振幅谱中的百分比.对零偏移距VSP资料而言，地震波传播时间为

对地面地震资料而言，炮检距xi处地震波传播时间为

由方程(12)，(13)，有t1，02，i2，i+1.此时，地面地震反射波的传播时间不再满足双曲线方程；对于深度H的地层而言，地面地震资料经动校正和叠加后的剖面，相对于VSP走廊叠加剖面出现下拉现象，无论是传播时间还是波形都不能与VSP资料匹配.

(2)仅反Q滤波振幅补偿

仅反Q滤波振幅补偿相当于利用方程(2)正演时，只考虑速度频散.此时，零偏移距VSP资料中，地震子波的振幅谱可表示为

地面地震资料中，炮检距xi处，地震子波的振幅谱可表示为

仅考虑速度频散时，地震波的传播时间与参考频率地震波传播的时间满足线性关系；假设参考频率的地震波传播时间为单位时间，与群速度对应的地震波的传播时间为α，且α>1.针对本文模型：

对零偏移距VSP资料而言，地震波传播时间为t1，0=αt0.

对地面地震资料而言，炮检距xi处地震波传播时间为t2，i=αti，且

此时，地面地震反射波的传播时间虽然满足双曲线方程，但来自深度H的地震波的传播时间都为αt0.经动校正和叠加后的地面地震资料，与VSP走廊叠加剖面虽在时深关系匹配，但波形不匹配.

(3)仅反Q滤波相位补偿

只进行反Q滤波相位补偿，相当于利用方程(2)正演时，只考虑能量衰减.此时，零偏移距VSP资料中，地震子波的振幅谱可表示为

地面地震资料中，炮检距xi处，地震子波的振幅谱可表示为

对零偏移距VSP资料而言，地震波传播时间为t1，0=t0；对地面地震资料而言，炮检距xi处地震波传播时间为t2，i=ti.

此时，地震波传播到任意时刻的相位与震源子波相同，地面地震反射波的传播时间满足双曲线方程，且来自深度H的地震波的传播时间都为t0.经动校正和叠加后的地面地震资料，与VSP走廊叠加剖面仅在时深关系和子波相位上可完全匹配.

(4)反Q滤波相位与振幅同时补偿

同时进行反Q滤波相位和振幅补偿，相当于声波方程正演.此时，零偏移距VSP资料和地面地震 资料中，地震子波的振幅谱完全相同，且都可表示为

对零偏移距VSP资料而言，地震波传播时间为t1，0=t0； 对地面地震资料而言，炮检距xi处地震波传播时间为t2，i=ti.

此时，任意时刻的地震波都和震源子波完全相同，任何频率的地震波都以参考频率的速度传播，且来自深度H的地面地震反射波传播时间都为t0.经动校正和叠加后的地面地震资料，与VSP走廊叠加剖面时深关系完全匹配，波形也完全匹配.

4 速度频散与反Q滤波相位补偿实例 4.1 速度频散实例

为验证地震波频散效应，并更加直观地理解速度频散现象，我们通过相同观测系统、炸药震源和可控震源激发采集的零偏移距VSP资料实例进行分析.

对零偏移距VSP记录Z分量进行初至拾取(炸药震源拾取起跳点，可控震源拾取波峰)，并进行波场分离得到下行波；在分离得到的下行波场中，选择初至附近的地震记录作为下行地震子波，并进行傅里叶变换得到振幅谱.地震子波选取原则为：尽量降低多次波的影响，使选择的地震子波的振幅谱尽可能光滑.

图 2a和图 2b分别为炸药震源和可控震源激发时，下行地震子波的振幅谱(峰值频率对应的振幅谱归一)；图 2中，纵轴为频率，横轴为道号(Trace number)，且道号随井下检波器的深度增大而增大；颜色代表能量.由于两次采集过程中，炮点位置相同，井下检波器的位置完全相同，因此相同深度的检波器记录到的地震波的传播路径相同，地层Q吸收造成的振幅衰减函数相同，速度频散函数也完全相同.

对比图 2a和图 2b中同深度处的下行地震子波的频谱可得，在道号为26-236对应的深层，炸药震源激发的震源子波较可控震源激发的震源子波而言，具有较宽的有效频带；但在道号为1-25对应的浅层，炸药震源激发的震源子波较可控震源激发的震源子波而言，具有较窄的有效频带；这主要是由震源子波差异造成.因为在VSP资料采集过程中，检波器都是从井底逐步往井口提升，而且受井下检波器级数的限制，每激发一次仅能测量某一个深度段的VSP记录，因此采集整个井段的VSP记录，需在同一个炮点多次激发，激发环境改变造成炸药震源激发震源子波主频逐渐降低，有效频带也变窄.而可控震源激发的震源子波的主频和有效频带，所受的影响较小.

图 3a中，初至1代表可控震源激发时各深度检波器对应的初至时间；图 3b中，初至2代表炸药震源激发各深度检波器对应的初至时间；图 3c和3d分别为井下相同深度的检波器对应的初至时间初至1与初至2的差值和比值.由图 3c可知，初至差值随深度增大呈增大趋势，但道号为1-25时，可控震源激发时对应的初至小于炸药震源激发时对应的初至时间，差值小于零；由图 3d可知，初至比值随深度增大逐渐增大，但变化率逐渐平缓；道号为1-25时，可控震源激发时对应的初至时间，小于炸药震源激发时对应的初至时间，初至比值小于1.

图 3 初至时间对比 Fig. 3 The comparision of real First Break

利用下行地震子波的振幅谱，按方程(12)进行初至反演，反演过程中：设fbi为深度Hi处下行波中参考频率的地震波的初至时间，vi为参考频率ωc的地震波的传播速度(单位为km/s)，Qi为炮点和检波点间地层的平均Q值，且ωc=300 Hz，vi=.

图 4中道号为1-25时，可控震源激发时反演的初至小于炸药震源激发时反演的初至，初至差值小于零，比值小于1；道号为26-236时，可控震源激发时反演的初至大于炸药震源激发时反演的初至，初至差值均大于0，初至比值均大于1，且都随深度增大呈逐渐增大的趋势，与实际情况吻合(图 2).但不同深度对应的反演初至时间受所用的平均Q值、参考频率ωc和参考频率ωc的相速度vi的准确度影响，仅大体趋势和实际资料吻合.但实例依然表明：对于同一个目的层，地震波的速度频散不仅与传播路径和时间有关，还与震源子波中各频率间相对能量关系有关，验证了公式(12)与(13)的合理性.

图 4 反演的初至时间对比 Fig. 4 The comparision of inverted First Break
4.2 井震标定实例

对另一口井采集的零偏移距VSP资料，经初至拾取、波场分离、动校正和走廊叠加等一系列的处理，得到了零偏移距VSP资料走廊叠加剖面.图 5为零偏移距VSP走廊叠加剖面镶嵌于井旁地面地震资料时间偏移剖面中，用于标定地面地震资料的处理成果.图 5中道号范围为1-51和68-118时，对应井旁地面地震时间偏移数据；道号范围为55-64时，对应零偏移距VSP走廊叠加数据(一道数据，重复十道显示)；道号范围为52-54和65-67时，对应空道，并用零值填充.由图 6可知，地面地震偏移剖面与VSP走廊叠加剖面标定结果较好，但也有部分同相轴与VSP走廊叠加剖面无法对应.

图 5 原始VSP走廊叠加剖面与井旁地面地震资料标定 Fig. 5 The calibration of VSP corridor stack and surface seismic data nearby the borehole

图 6 经过叠后反Q滤波相位补偿后的标定 Fig. 6 The result of Fig. 1 after inverse Q-filter phase compensation

利用该井采集的零偏移距VSP资料，进行层Q值反演，并利用反演的层Q值，对图 5中VSP走廊 叠加数据和井旁的地面地震时间偏移数据同时进 行叠后反Q滤波相位补偿，且处理参数完全相同，补偿效果如图 6所示.对比图 5和图 6可得，若图 5中地面地震剖面的反射层与VSP走廊叠加剖面的匹配关系较好，则反Q滤波相位补偿后，这种匹配关系更好；若原始匹配关系不好，经反Q滤波相位补偿后，绝大部分匹配关系较好，即使有部分反射层对应关系还不是很理想，但匹配关系已有较大提升.

对图 5和图 6的数据分别计算井旁地面地震数据与VSP走廊叠加数据之间的互相关系数，以更加 直观地对比反Q滤波相位补偿前后地面地震数据 与VSP走廊叠加数据之间的匹配程度.图 7显示了 道号范围为40-74时，反Q滤波相位补偿前后，地面地震数据与VSP走廊叠加数据的互相关系数，互相关系数1代表反Q滤波相位补偿前地面地震数据与VSP走廊叠加数据的互相关系数，互相关系数2代表反Q滤波相位补偿后地面地震数据与VSP走廊叠加数据的互相关系数.为方便显示，将空道与VSP走廊叠加数据的互相关系数用1显示.经反Q滤波相位补偿后，道号范围为40-51和68-74的地面地震数据与VSP走廊叠加数据的互相关系数都增大，如图 7a；互相关系数比值也都大于1，如图 7b.

图 7 反Q滤波相位补偿前后地震剖面与VSP走廊叠加剖面的互相关系数对比 Fig. 7 The cross correlation coefficient before and after inverse Q-filter phase compensation of Fig. 5

综合图 5-7，反Q滤波相位补偿后，地面地震资料与VSP资料走廊叠加剖面的吻合度得到有效提升.

5 结论与建议

本文所推导的地面地震资料与零偏移距VSP资料之间的速度频散关系式，可以用于研究群速度之间的对应关系.

实际资料处理中，我们需对地面地震资料和VSP资料首先做好反Q滤波相位补偿工作，从而尽可能消除地震波的速度频散效应，有效提升地面地震资料与VSP资料的匹配度，最终提高地震资料成果可信度.

参考文献
[1] Badri M, Mooney H. 1987. Q measurements from compressional seismic waves in unconsolidated sediments.   Geophysics, 52(6):772-784.
[2] Bano M. 1996. Q-phase compensation of seismic records in the frequency domain.   Bulletin of the Seismological Society of America, 86(4):1179-1186.
[3] Bickel S H, Natarajan R R. 1985. Plane-wave Q deconvolution.   Geophysics, 50(9):1426-1439.
[4] Blias E. 2012. Accurate interval Q-factor estimation from VSP data.   Geophysics, 77(3): 149-156.
[5] Chen Z B, Chen X H, Li J Y, et al. 2014. A band-limited and robust inverse Q filtering algorithm. OGP (in Chinese), 49(1):68-75.
[6] Futterman W I. 1962. Dispersive body waves.   Journal of Geophysical Research, 67(11): 5279-5291.
[7] Gao J H, Yang S L. 2007. On the method of quality factors estimation from zero-offset VSP data. Chinese J.   Geophys (in Chinese), 50(4): 1198-1209.
[8] Gao J H, Yang S L, Wang D X. 2008. Quality factor using instantaneous frequency at envelope peak of direct waves of VSP data. Chinese J. Geophys (in Chinese), 51(3): 853-861.
[9] Hamilton E L. 1972. Compressional-wave attenuation in marine sediments.   Geophysics, 37(4): 620-646.
[10] Hargreaves N D, Calvert A J. 1991. Inverse Q filtering by Fourier transform.   Geophysics, 56(4): 519-527.
[11] Hauge P S. 1981. Measurements of attenuation from vertical seismic profiles.   Geophysics, 46(11): 1548-1558.
[12] Kjartansson E. 1979. Constant Q wave propagation and attenuation.   Journal of Geophysical Research, 84(B9): 4737-4748.
[13] Liu C, Feng X, Zhang J. 2013. A stable inverse Q filtering using the iterative filtering method. OGP, 48(6): 890-895.
[14] Sams S M, Neept J P, Worthington M H, et al. 1997. The measurement of velocity dispersion and frequency-dependent intrinsic attenuation in sedimentary rocks.   Geophysics, 62(5): 1456-1464.
[15] Spencer T W, Sonnadt J R, Butlers T M. 1982. Seismic Q-stratigraphy or dissipation.   Geophysics, 47(1): 16-24.
[16] Stainsby S D, Worthington M H. 1985. Q estimation from vertical seismic profile data and anomalous variations in the central North Sea.   Geophysics, 50(4): 615-626.
[17] Toksoz M N, Johnston D H, Timur A. 1979. Attenuation of seismic waves in dry and saturated rocks: 1. Laboratory measurements.   Geophysics, 44(4): 681-690.
[18] Tonn R. 1991. The determination of the seismic quality factor Q from VSP data: A comparison of different computational methods.   Geophysical Prospecting, 39(1): 1-27.
[19] Wang X J, Yin X Y. 2011. Estimation of layer quality factors based on zero-phase wavelet. Progress in Geophys.   (in Chinese), 26(6): 2090-2098.
[20] Wang Y H. 2002. A stable and efficient approach of inverse Q filtering.   Geophysics, 67(2): 657-663.
[21] Wang Y H. 2003. Quantifying the effectiveness of stabilized inverse Q-filtering.   Geophysics, 68(1): 337-345.
[22] Wang Y H. 2006. Inverse Q-filter for seismic resolution enhancement.   Geophysics, 71(3): 51-60.
[23] Xu C D, Stewart R R. 2006. Seismic Attenuation (Q) Estimation from VSP Data.   CSEG National Convention, Extended Abstracts, 268-277.
[24] Yan H Y, Liu Y. 2009. Estimation of Q and inverse Q filtering for prestack reflected PP-and converted PS-waves.   Applied Geophysics, 6(1): 59-69.
[25] Yao Z X, Gao X and Li W X. 2003. The forward Q method for compensating attenuation and frequency dispersion used in the seismic profile of depth domain.   Chinese Journal of Geophysics (in Chinese), 46(2): 229-230.
[26] Zhao J, Gao J H, Wang D X, et al. 2013. Estimation of quality factor Q from pre-stack CMP records. Chinese J. Geophys. (in Chinese), 56(7): 2413-2428.
[27] 陈增保, 陈小宏, 李景叶等. 2014. 一种带限稳定的反犙滤波算法.   石油地球物理勘探, 49(1): 68-75.
[28] 高静怀, 杨森林. 2007. 利用零偏移VSP 资料估计介质品质因子方法研究.   地球物理学报, 50(4): 1198-1209.
[29] 高静怀, 杨森林, 王大兴. 2008. 利用VSP 资料直达波的包络峰值处瞬时频率提取介质品质因子.   地球物理学报, 51(3): 853-861.
[30] 刘财, 冯晅, 张瑾. 2013. 稳定的迭代法反Q滤波.   石油地球物理勘探, 48(6): 890-895.
[31] 王小杰, 印兴耀. 2011. 基干零相位子波地层Q值估计.   地球物理学进展, 26(6): 2090-2098.
[32] 姚振兴, 高星, 李维新. 2003. 用于深度域地震剖面衰减与频散补偿的反Q滤波方法.   地球物理学报, 46(2): 229-230.
[33] 赵静, 高静怀, 王大兴等. 2013. 利用叠前CMP资料估计介质品质因子.   地球物理学报, 56(7): 2413-2428.