0 引言

页岩气甜点区的地球物理预测是页岩气勘探的重要技术环节，其中强非均质多孔介质的地球物理响应特征是关键, 为根据观测波场的散射特征来分析地下介质的非均匀性提供依据.作为岩石物理超声实验的补充手段，数字岩芯超声响应数值模拟是深入理解弹性波在非均质多孔岩芯中传播机理的重要方法.本文基于Biot双相介质方程，采用不分裂卷积完全匹配层与旋转交错网格有限差分方法，结合超声实验数据，开展龙马溪组页岩数字岩芯超声响应数值模拟研究和散射特征分析.虽然Biot双相介质方程在解释非均质多孔岩芯中弹性波传播机理还存在明显不足，如低估衰减和频散，但仍具有不可动摇的坚实地位和广泛的应用价值(Mavko et al., 2009).

作为地震记录中直达波之后的一段连续波形，尾波被广泛用于分析岩石圈小尺度非均匀介质的散射强度，反演介质的非均匀尺度(Aki, 1969; Aki and Chouet, 1975)，其中大部分研究主要基于台站观测地震数据的长波长尾波(1~30 Hz)(Sato, 1977; Wu and Aki, 1985)，而基于短波长尾波的非均质散射特征分析和应用相对较少(Fehler, 1982; Wu, 1982).基于物理模拟地震数据，高频尾波散射分析可以探测厘米级小尺度非均匀介质(Matsunami, 1991; Nishizawa et al., 1997; Sivaji et al., 2002).超声岩石物理实验可以提供更高频的观测记录，然而长期以来，超声岩石物理实验只利用直达波研究岩石的速度和吸收衰减特征, 其后续波(特别是尾波)由于其成份的复杂性一直未得到利用，近年来基于超声尾波开展的毫米至微米级小尺度非均匀介质散射分析(Grêt et al., 2006; Guo and Fu, 2007; Wei and Fu, 2014; Yoshimitsu et al., 2016)，其可靠性受到了相关研究人员的质疑.超声尾波的复杂性源于压电换能器的辐射声场特性、岩芯有限尺寸影响、岩芯边界反射干涉、顶底边界震荡干涉、波型转换影响等(Stacey and Gladwin, 1981).因此，有必要通过数字岩芯超声响应数值模拟来厘清这些因素对超声尾波的影响.本文拟根据页岩岩芯超声实验的换能器辐射声场和岩芯尺寸等相关参数，开展该页岩数字岩芯高精度有限差分数值模拟，通过吸收程度可调节的数值吸收边界，控制岩芯边界反射对超声尾波的干涉影响，进而分析页岩小尺度非均匀质散射特征.

基于孔弹方程的有限差分数值模拟得到了广泛的研究(Grêt et al., 2006; Guo et al., 2009; Guo and Fu, 2007; Wei and Fu, 2014; Yoshimitsu et al., 2016)，这些模拟研究存在的不足主要表现为：(1)双相介质参数模型局限于理想的流变模型，而非实际的数字岩芯模型；(2)模拟研究侧重算法，不考虑岩芯尺寸和超声实验环境影响，缺乏与超声实验结果的比较(Arntsen and Carcione, 2001; Gurevich, 1996; Gurevich et al., 1999); (3)对于超高频(MHz)超声频段和强非均质岩芯，算法精度有待提高.基于典型砂岩的数字岩芯模型，Zhang等(2014)利用数字岩芯像素重建双相介质参数模型，采用高精度的不分裂卷积完全匹配层旋转交错网格有限差分方法解决岩芯侧壁反射问题和600 kHz高频计算精度问题，模拟结果与该岩芯600 kHz超声实验波形数据(含超声尾波)进行了比较，目的是分析超声实验过程中岩芯侧壁反射对超声尾波的影响.本文拟开展龙马溪组页岩数字岩芯超声响应数值模拟及非均质分析，面临的挑战主要包括：(1)与砂岩的数字岩芯不同，面向页岩气勘探的龙马溪组页岩数字岩芯具有更小的微观尺度(图像分辨率达0.65 μm)，需要在更高的超声探测频率(MHz)来研究的波与微结构(矿物和孔隙)的相互作用机制，介质表现出更强的非均匀特征，挑战模拟算法的精度和数值稳定性；(2)数字岩芯只是提供岩石微观结构(矿物和孔隙)的图像灰度分布，为了确保与超声实验波形的可比较性，需要精细的孔弹性参数数值建模，确保建立的双相介质非均质数字模型基本反映实验岩芯的岩石物理特征；(3)兆赫兹超声波在岩芯侧壁引起较强的边界反射，严重干涉超声尾波，影响散射分析精度，因此需要更高精度的吸收边界模拟算法和精细的技术实施方案.

高分辨率页岩数字岩芯真实地再现了岩芯内部矿物及孔隙的分布，利用图像阈值分割算法，将岩芯矿物成分主要分为石英类、黏土类、黄铁矿及孔隙四种类型，根据矿物的岩石物理参数，利用数值计算方法(Arns et al., 2002)建立精细的孔弹性非均质双相介质模型.本文采用高精度的不分裂卷积完全匹配层旋转交错网格有限差分方法模拟Biot双相介质方程，与常规的交错网格有限差分算子不同，旋转交错网格不需要平均弹性模量，有效的改善强非均质引起的数值不稳定性，提高了孔隙与矿物颗粒间弹性波的模拟计算精度(Saenger et al., 2005a).不分裂卷积完全匹配层是目前解决临界入射边界反射最有效的方法(Komatitsch and Martin, 2007; Martin and Komatitsch, 2008；张鲁新等, 2010)，高分辨数字模型的网格剖分极小，模拟频率极高，需要不断调整时间及空间步长使其算法收敛并降低频散，同时灵活控制吸收边界层数，减少边界反射对尾波的干扰以匹配超声波实验波形.岩芯数值模拟记录与超声实验记录的对比分析表明，数字岩芯超声波响应反映了弹性波在非均质多孔岩芯中的传播特征.本文进一步利用模拟波形记录的尾波成分对龙马溪组页岩散射问题进行了半定量分析，结果表明，尾波散射Q值计算的非均匀尺度与数字岩芯图像本身的非均匀尺度基本吻合.当介质的非均匀尺度和超声波的波长尺度相当时，介质引起超声波的散射最强，与散射理论预测一致.高分辨率页岩数字岩芯超声响应表明，超声波的频率尺度只有在几百兆赫兹以上才能产生强散射，即弹性波与页岩介质的相互作用表现很强的尺度效应.

1 基于数字岩芯的多孔介质模型

通常, 岩石物理实验不能定量的控制、观察，计算流体分布和岩石弹性性质，基于研究其微观因素影响的困难性，数字岩石物理学应运而生，数字岩芯技术能够真实的描述岩芯内各类矿物、孔隙的空间分布.构建数字岩芯方法主要有X射线扫描法和二维信息重建法，岩芯实验室通过观察流体流动和岩石弹性性质，获得岩石属性即弹性模量、输运性质、电性性质等(Dvorkin et al., 2003; Kewen and Ning, 2008; Liu et al., 2009; Saenger et al., 2005b).目前，数值模拟更多的应用于砂岩和泥质砂岩岩芯研究中(Andrä et al., 2013a; Madonna et al., 2013)，Liu通过X射线CT扫描得到的数字岩芯(Liu et al., 2009)，描述了岩石微孔隙结构并利用开形态算法模拟在驱油过程中孔隙内油和地层水的分布，张艳等人利用CT扫描得到的砂岩数字岩芯数据，建立初始的双相介质模型，利用尾波散射方法，探讨砂岩内非均匀性特征(Zhang et al., 2014).Josh提出了页岩储层的实验室测试技术(Josh et al., 2012).页岩具有从纳米到分米大跨度的非均匀尺度特征，对成像技术的分辨率和数值模拟算法的精度是极大的挑战，较少的学者开展非均质页岩数值模拟工作.本文试图通过X射线CT扫描技术获得较为精准的龙马溪组页岩微观结构图像，提取岩芯内各矿物的岩石物理属性，基于Biot理论，建立切合真实岩芯的二维双相介质模型，为超声波波数值模拟及散射特征分析等后续工作做准备.

1.1 X射线CT成像

本文对龙马溪组页岩岩芯进行扫描，每个X射线双晶单色成像光束提供的光子能量约为8~72.5 keV，光束尺寸为45mm×5 mm.X射线CT扫描成像可以从六个自由度对样品扫描.样品的图像重建使用有序子集期望最大化算法来提高成像精度.本研究中龙马溪组岩心成像的样本直径约1 mm，长约3 mm的大小，X射线透射成像，单色光束为20 keV，样本与探测距离为10 cm.重建图像分辨率为0.65 μm，曝光时间为0.65 s.

CT实验获得的岩芯数据量巨大，是由3200000张CT扫描横截面切片合成叠加而成的高分辨率三维数字图像数据体，分辨率为600×600×700像素，大小为390 μm×390 μm×455 μm.图 1a所示为截取该三维数据体中心的一块，颜色表示不同密度物质对应的X射线吸收率，红色为最强(黄铁矿颗粒)，蓝色为最弱(例如细长状的裂缝).图 1b为其中一个CT扫描XY横截面切片，为了对比不同矿物的特征，将32位彩图绘成8位灰度图，图 1c为图 1b截取的正方形区域，截取的切片尺寸大小为800 μm×800 μm.

我们根据各类矿物的X射线吸收率、泊松比、弹性模量等性质，通过阈值分割方法将岩石分成孔隙、石英类、黏土类、黄铁矿四类物质，图 2分别为图 1c对应的石英类、黏土类、黄铁矿及孔隙的分布图.

1.2 建立页岩数字岩芯双相介质模型

传统的岩石物理模型是基于经验公式或理论模型，而数字岩芯通过高分辨率图片复原方法可直接观察矿物与矿物或矿物与孔隙之间的接触关系.图 3所示为通过图 2的阈值分割方法并赋予四类物质不同颜色，建立的页岩数字岩芯岩石模型，其尺寸大小为800 μm×800 μm，分辨率为600×600个像素点，直观再现了孔隙和不同矿物分布.

为了与Biot孔弹性介质波动方程的参数对应，需要知道每一个像素点是固相或液相.根据文献(Zhang et al., 2014)的方法，为了体现流体相的等效影响，假定液相(油)是均匀分布在整个介质模型中，即每个像素点都是双相介质，既有固相参数，又有液相参数，据此将数字岩芯岩石模型转化成双相介质模型.

数值模拟需要在岩芯双相介质模型上每一个像素点定义合理的物理属性，根据矿物间的比例，计算平均模量从而简化弹性参数计算(Von, 1928)，石英的孔隙一般在0.1%~8.7%范围内(刘佑荣和唐辉明, 1999)，我们选取一个平均值5%，黏土矿物及石英类矿物的比例约为1:2，根据岩石总孔隙度为4.5%，计算得到两种矿物分布的孔隙度分别为4.25%和5%.根据总渗透率为9.768 mD，计算得到黏土类矿物和石英类矿物分布的渗透率分别为11.24 mD和6.81 mD.根据弯曲度因子v=1-r(1-1/ϕ)，而椭圆形颗粒r介于0~1之间(Mavko et al., 2009)，故石英类和黏土类矿物的弯曲度因子分别为11.6和12.3，体积模量、剪切模量及密度等基本岩石物理参数可参考岩石物理学手册(Mavko et al., 2009)，详见表 1.

合理的岩芯双相介质模型参数对于数值模拟的实际意义至关重要.本文采用不分裂卷积完全匹配层与旋转交错网格有限差分方法(公式详细推导可见附录)，对数字岩芯双相介质模型进行孔弹介质方程超声波数值模拟.

2 数值模拟及其实验比较 2.1 数值模拟程序验证

为验证有限差分程序的正确性，选取双层均匀模型，其中上层为均匀各项同性双相介质470mm×400 mm，下层为均匀各项同性弹性介质330mm×400 mm.计算空间间隔Δx=Δz=1.0 mm，时间间隔dt=0.0001 ms.吸收边界层为10个网格，完美匹配层参数m=2, R=10^-6, X_max=1, α_max=120.激发震源位于模型中心x=400 mm，z=200 mm，为固相纵波源，采用的雷克子波中心频率为120 Hz型详细参数见表 2.

图 4为孔弹模拟波在双层模型中的传播波场快照，可识别出快纵波P¹、快纵波反射P¹P¹、慢纵波P²、横波转换波P¹S、快纵波转换波P²P¹、慢纵波转换波P¹P²、慢纵波反射波P²P².纵波相对于横波有明显的振幅响应，固相中不存在横波转换波.

2.2 页岩岩芯数值模拟

利用已验证的有限差分程序对建立的数字岩芯双相介质模型进行超声波孔弹模拟，模型大小为800 μm×800 μm，计算空间间隔Δx=Δz=1.5385 μm，时间间隔dt=1.5×10^-11s，时间长度Nt=10000，V_t=4506 m·s^-1, 高斯子波的主频f₀=600 MHz，CMPL参数, m=2, R=10^-6, X_max=1, α_max=40，震源用星号标记(200Δx, 45Δz)，接收点用三角形标记(200Δx, 315Δz).图 5为采用不同吸收边界层数的超声数值模拟波场快照，可见当吸收边界层np=1时，存在明显边界反射，随着吸收边界层数的增加，边界反射的影响逐渐减弱.吸收边界层数增至np=5时，边界反射基本消除，极大减弱对尾波造成的干扰影响，有利于正确估计数值模拟尾波散射Q值，可利用数值模拟方法究超声波实验中边界反射对实验尾波散射Q值的影响.

2.3 实验与数值模拟对比

本文对重庆礁石坝地区下志留系龙马溪组页岩样本进行了常温常压(室温和标准大气压)下干燥岩样的600 kHz超声实验，页岩样本尺寸为100 mm×200 mm的圆柱体，层理垂直于岩石样品的轴向.根据实验波形数据，该样品密度为2650 kg·m^-3，纵波速度为4444 m·s^-1，横波速度为2777.78 m·s^-1.

图 6为超声波在铝样中传播的P波波形，其直达波可近似为数值模拟的激发子波，如图 7所示.

对比分析数值模拟和超声波实验的波形数据，有助于理解超声波在岩芯中的传播过程.图 8对比了超声实验波形与不同吸收边界层数的数值模拟P波记录.当吸收边界为1个单位，模拟波形存在明显的边界反射，导致尾波散射Q_pc值(Zhang et al., 2014; Guo et al., 2009)计算偏大，随着吸收边界层数的增加，数值模拟边界反射逐渐削弱，计算的尾波散射Q值与实验结果趋于一致.

一般而言，超声实验激发波形并非完全平面波，存在边界反射，但由于岩样常常被包裹一层橡胶皮进行岩石物理实验，一定程度削弱了边界反射强度及其对尾波的影响.基于Biot理论的孔弹数值模拟严重低估超声波在双相介质中的内耗衰减，导致数值模拟和超声实验的波形有些偏差.

3 不均匀尺度和波长尺度关系 3.1 介质不均匀尺度计算

目前常用的随机介质的类型有高斯型、指数型和冯卡曼型，高斯型介质扰动变化较圆滑，指数型相比高斯型在介质小尺度变化更加剧烈，冯卡曼随机介质的变化最为剧烈.为了探讨介质不均匀尺度和波长尺度的相互作用，我们假定实际岩芯为指数型或高斯型随机介质，计算反映岩芯非均匀程度的相关长度参数a.

在表 3中，ACF为自相关函数，P(k)为能谱密度，其中ε²=F(0)，Γ为修正的贝塞尔函数，κ是赫斯特数.自相关函数表征介质扰动的相关性统计特征，描述介质参数的相似度与滞后距离r之间的关系.

图 9为根据页岩数字岩芯数据进行自相关计算得到的自相关函数随滞后距离的变化曲线，可见曲线开始阶段变化基本满足指数型函数特征，其相关长度a约为6.276 μm，由于实际岩芯介质并非完全随机介质，曲线后段变化波动较大，偏离指数型函数特征.

3.2 散射态式图

根据页岩数字岩芯内波传播距离长度L及其自相关函数变化的相关长度参数a，假定激发超声频率范围为200~1000 MHz, 传播距离范围在0.502×10^-4~2.295×10^-4m，可计算得到该页岩岩芯的散射强度分布，位于L/a-ka散射态势图(吴如山, 1989; 尹军杰，2005)中的强散射区域，如图 10中的灰色阴影区域所示.该图反映了介质不均匀尺度和波长尺度的相关性，位于强散射区的灰色阴影区域需要数值模拟方法进行散射分析，例如图中星号表示主频为600 MHz，传播距离为2.295×10^-4m的超声波在该岩芯中传播引起的散射强度，这可以从下节数值模拟的尾波散射Q值分布得到证实.改变数值模拟中激发超声主频和接收距离，计算模拟尾波散射强度随主频的变化曲线，可评估页岩岩芯的非均质强度与超声波长尺度的相关性.

3.3 尾波Q值

超声波在岩芯中传播时的衰减包括散射衰减和吸收衰减，本文主要讨论散射衰减.采用完全匹配层吸收边界技术，削弱边界反射对尾波的干扰.改变超声子波主频和接收距离，开展页岩岩芯超声波数值模拟研究，计算模拟尾波的散射Q值，绘制散射Q值随主频的变化曲线，据此评估页岩岩芯的非均质性，为页岩储层声学非均质预测提供依据.如图 11所示为采用不同超声子波主频模拟的P波记录，由于尾波振幅相对直达波振幅极弱，所以我们保持直达波振幅不变，同时调节尾波增益至初始的10倍, 从而凸显尾波的存在，其中图中竖线标记为直达波和尾波的分界线.根据文献(Zhang et al., 2014; Guo et al., 2009)截取模拟记录尾波的起止时间，采用Sato散射模型(Sato, 1977)计算尾波散射Q值，图 12为尾波散射逆品质因子值随子波主频变化的曲线图.可见，当主频波长尺度和介质的相关长度参数相当时，散射强度达到最大.

4 结论

页岩气甜点区多孔介质的强非均质是制约地球物理预测的关键因素，深入理解强非均质多孔介质的地球物理响应特征是页岩气地球物理勘探的基础，为根据观测波场的散射特征来分析甜点区的非均匀性提供依据.超高频(几百兆赫兹)超声数值模拟微米至纳米尺度的马溪组页岩数字岩芯及其强非均质性严重挑战数值模拟算法的精度和数值稳定性，作为岩石物理实验的辅助手段，本文基于Biot双相介质方程，采用不分裂卷积完全匹配层与旋转交错网格有限差分方法，结合超声实验数据，开展龙马溪组页岩数字岩芯超声响应数值模拟研究和散射特征分析.主要研究成果和结论如下：

(1) 合理的岩芯双相介质模型参数对于数值模拟的实际意义至关重要，本文利用图像阈值分割算法将龙马溪组页岩数字岩芯主要成分分解为石英类、黏土类、黄铁矿及孔隙四种类型，假定液相(油)均匀分布在整个介质模型中，根据数字岩芯的孔隙度、渗透率和各类矿物的岩石物理参数，建立精细的非均质双相介质模型.

(2) 基于Biot孔弹性方程的旋转交错网格有限差分方法(空间8阶和时间2阶)虽然严重低估双相介质的内耗衰减，但能精确模拟超声波在页岩岩芯中传播的散射衰减.

(3) 采用不分裂卷积完全匹配层吸收边界能精确控制数值模拟的边界反射量及其对尾波的干涉，通过与实验尾波散射Q值的比较，可了解超声实验中的边界反射量.

(4) 通过不同超声子波主频数值模拟试验，600 MHz子波主频对应的散射衰减达到最大，更低和更高子波主频的散射衰减逐渐减弱，结合L/a-ka散射态式图分析表明，本文采用的龙马溪组页岩数字岩芯(分辨尺度为800×800 μm)的非均质强度与600 MHz波长尺度相当.

总之，开展页岩岩芯超声波数值模拟研究，计算模拟尾波的散射Q值，绘制散射Q值随主频的变化曲线，据此评估页岩岩芯的非均质性，为页岩储层声学非均质预测提供依据.

附录 “孔弹介质方程及旋转交错网格有限差分方法” 附录A “Biot方程”

Biot方程，设u_i和ω_i分别是固体骨架位移和流相相对位移，则忽略震源项的各项同性流体饱和孔隙介质Biot方程

(A1)

(A2)

其中ρ_b和ρ_m由固体颗粒密度ρ_s和流体密度ρ_f确定，其中，ν是描述孔隙结构特征的弯曲度因子，各项同性孔隙介质中固体体积模量为K_s，液相体积模量为K_f；干岩石骨架的体积模量为K_d，剪切模量为μ，则biot系数β，耦合模量M，饱和岩石的拉梅系数λ_u，和黏滞系数b由下列式子给出，

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

其中，η是孔隙流体的黏滞系数，κ是渗透率.对均匀孔隙弹性性质，其特征频率为

(A9)

附录B “速度应力格式”

波动方程数值模拟中通常将二阶方程改写为一阶格式，如速度—应变格式、速度-应力格式、等本文采用速度应力格式，另固体速度，液相相对速度，总应力为，孔隙流体压力为p，则Biot方程可表示:

(B1)

(B2)

(B3)

(B4)

(B5)

附录C 完美匹配层吸收边界

引入复数因子将直角坐标系扩展为复数坐标：

(C1)

其中是衰减因子，ω是角频率，则对应的微分算子即，其中，

(C2)

不分裂的卷积完全匹配层吸收边界首先给出经改造的扩展函数的形式：

(C3)

其中χ_x≥1，α_x≥0，然后分别介绍用迭代卷积法和微分方程法推到迭代格式CPML吸收边界的过程，最后给出CPML吸收边界中参数设置的一般原则.

下面只介绍微分方程法

本小节对文献(Qin et al., 2009)中公式推导过程进行了部分修正.将代带入公式中

(C4)

其中，可写成如下形式：

(C5)

对上式做逆Fourier变换，有

(C6)

其中ψ_x是的逆Fourier变换，形如的方程具有以下通解：

(C7)

其中，c为热力常数，将上式进行离散，得到

(C8)

通过迭代格式可将常数c消去，则有

(C9)

有上式可知，ψ_x可有下式迭代求解：

(C10)

采用迭代格式的CPML在实现过程只需要一个额外的辅助变量ψ_x，可以极大地提高计算效率，既避免了场分解的众多方程，又不在时间域显式的进行卷积运算，对空间导数做如下的替换可将波动方程过渡到PML区域波动方程：，参数设置，对d_x作如下设置：

(C11)

(C12)

其中，l(0≤L≤M)为PML内的计算点到PML区域内边界的距离；L为该PML区域的厚度；m为多项式的阶数，一般取为2或3；R为理论反射系数，本文取10^-6.

令χ_x和α_x分别为

(C13)

(C14)

当α_x=0，χ_x=1时，s_x与常规PML一致，PML区域内α_x在0~α_max间线性变化，通常取α_max=f₀，其中f₀是子波频率.

附录D 旋转交错网格有限差分不分裂卷积完全匹配层吸收边界的实施

RSG技术沿着空间网格的对角线方向计算空间导数，在将计算结果线性内插得到坐标轴方向的空间导数，得到网格指教坐标系下坐标轴方向的空间导数，令，则对角线方向x和z可表示为

(D1)

(D2)

直角坐标下的坐标轴方向的一阶空间导数为

(D3)

(D4)

下面给出2D孔隙弹性介质速度-应力格式方程的旋转交错网格与不分裂卷积完全匹配层吸收边界的实施的详细过程：

固液相各速度及应力分量的旋转交错网格有限差分的离散格式及实施细节

(D5)

(D6)

(D7)

(D8)

(D9)

(D10)

(D11)

其中c_n是差分系数，L表示该差分近似是L阶的，对上式求取各参数的辅助变量，然后可求出复数坐标下的空间导数(即施加CMPL)，再代入公式B，可得到

(D12)

(D13)

(D14)

(D15)

(D16)

(D17)

(D18)

(D19)

附录E 尾波散射Q值的计算公式

尾波散射Q值计算采用Aki和Chouet(Aki and Chouet, 1975)提出的单次散射模型研究尾波振幅随时间的衰减，单次散射模型是描述尾波衰减最早也是最简单的模型，当离逝时间比平均自由程小得多的时候，将尾波振幅和散射衰减联系起来.该模型用于描述三维空间中散射能量密度的时间依赖性以及当尾波离逝时间小于平均自由程的情况下散射和尾波振幅的关系.由单次散射模型，在一定频率下，尾波振幅与时间的函数关系可写为

(E1)

式中K(f)是频率f下的尾波震源因子.与时间和辐射模式无关，Q_c(f)是尾波的品质因子，α是几何扩散参数，对体波、面波和散射波分别取值为1.0, 0.5和0.75(Sato, 1977; Wu and Aki, 1985)由于尾波属于逆散射体波，在公式E1中取α=1，并取对数

(E2)

Q_c可由最小二乘拟合得到的lnA(f, t)与t的线性关系式得出.

致谢

谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.感谢两位审稿专家对论文提出宝贵的修改意见.