横向各向同性地层随钻多极子测井声场的分波分析与特性研究

饶博， 唐晓明,3， 苏远大,3*

1 中国石油大学(华东)地球科学与技术学院， 青岛 266580 2 中国石油大学(华东)深层油气重点实验室， 青岛 266580 3 青岛海洋科学与技术试点国家实验室， 青岛 266580

0 引言

近年来，深水、深层钻井逐渐增多，水平井及大斜度井的数量明显增加，因此对随钻测井技术的需求更加迫切.在随钻测井过程中，地层纵横波速度的获取对地层参数的反演以及指导水力压裂有着十分重要的作用，对此测井界研发了单极子纵波(Schmitt，1988；Minear et al.，1995)和四极子横波(Tang et al.，2002)的随钻多极子声波测井技术.随着以页岩油气为代表的非常规油气资源勘探开发的深入(吕庆达等，2018)，地层各向异性对随钻声波测量的影响变得十分重要.为此胡恒山和何晓(2009)详细推导了横向各向同性介质中的波动方程，求解了VTI地层中裸眼井声波的频散方程，研究了伪瑞利波、弯曲波以及螺旋波的频散特征.何晓等(2010)针对VTI地层给出了满足辐射条件的井孔声场复数解析函数的单值定义，据此进行了裸眼井多极子声场的分波计算.在随钻声波测井中，李希强等(2013)建立随钻VTI声波测井模型，通过模式分析法研究了TI地层中随钻井孔内一些模式波的传播特性，指出与裸眼井的情况类似，随钻斯通利波也可以用于反演地层水平向横波速度；许松(2015)在此基础上提出一种同时反演地层水平向和垂直向横波速度的联合反演方法，进而求取地层的各向异性；王瑞甲和乔文孝(2015)进一步研究了TI地层中各向异性对随钻四极子声波测井的影响，指出地层各向异性对地层四极子模式波的影响较大，建议在实际数据处理时，考虑采用基于数据驱动的频散校正.

综上所述，研究随钻条件下VTI地层的声传播规律具有十分重要的意义.但由于VTI地层随钻声场中钻铤模式波与地层模式波耦合复杂，难以从全波波形中直接区分.针对这一问题，本文采用分波计算的方法来深入探究地层模式波的传播特征.为此分别对单极、偶极以及四极子声源激发的随钻声场进行了分波与全波计算；通过与电缆测井的结果比较，考察了钻铤对地层模式波的影响，发现在少数参数满足δ>ε+c44/2c33的地层中，伪瑞利波、弯曲波及四极子螺旋波的截止频率速度小于真实地层垂向横波速度，应结合其他各向异性参数来确定该横波速度.计算结果为各向异性地层中随钻声波测井解释和地层评价提供了理论依据.

1 横向各向同性地层随钻多极子声场理论

图1给出了横向各向同性(VTI)地层中随钻声波测井模型示意图，分别为模型三维图(图1a)和井孔横截面图(图1b)，由内到外依次是钻铤内流体、钢钻铤、钻铤与地层之间的流体环以及井外VTI地层，图中rc代表钻铤内径，rb代表钻铤外径，ra代表井孔半径.

图1 VTI地层中随钻声波测井模型(a)及井孔截面图(b)Fig.1 Schematic diagram of acoustic model for logging while drilling in a VTI formation (a) and the cross sectional view of the model (b)

1.1 井孔内含钻铤的声场

根据柱状分层介质声场理论(唐晓明和郑传汉，2004；崔志文，2004)，钻铤位移势Φtool、tool和Γtool及其内外层流体的位移势Φin、Φout表达式为：

(1)

1.2 井孔外VTI地层的声场

TI固体的弹性性质由c11，c13，c33，c44，c66五个弹性常数来表达，若TI介质对称轴与井轴平行，则柱坐标系下的应力-应变关系为：

(2)

其中σij和eij(i,j=r,θ,z)分别是TI固体中的应力和应变.

(3)

式中，ΦVTI、VTI、ΓVTI分别表示固体中纵波、SH横波以及SV横波的势函数；jz为轴向单位矢量.

将式(3)代入弹性波动方程并且结合位移-应变关系以及上述TI介质的本构关系得到对应于三个位移势函数ΦVTI、VTI、ΓVTI的三个波动方程，具体表达式为：

(4)

方程组中前两个是P波和SV横波的耦合波动方程，第三个是纯SH横波的波动方程，与P-SV波动不耦合.

c66q2-c44k2+ρω2=0,

(5)

(q2-k2)(Uq4+Vω2q2+Wω4)=0,

(6)

其中：

式(5)只与SH横波有关，其解为：

(7)

式(6)为两个因式的乘积，第一个因式无意义，其解不代表任何波动；第二个因式代表了P波和SV横波的耦合，二者的解分别为：

(8)

将式(7)、(8)代入式(1)地层波势函数的通解，得到其解析表达式：

(9)

系数a′和b′由式(10)给出，其中ρ为地层密度：

(10)

在此基础上得到TI地层的位移、应变分量，根据式(2)求得应力分量，具体表达式见附录A.

1.3 边界条件

以上理论推导给出了图1所示模型各部分的声场，现将声场由界面处的边界条件衔接；柱坐标下固、液界面的边界条件是径向位移和应力连续和切向应力为零(唐晓明和郑传汉，2004).参考郑晓波(2017)论文，通过径向应力连续的方式将位于钻铤外表面的声压源ps引入随钻声场，以实现声场的激发.

钻铤内流体和钻铤界面(r=rc)处：

(11)

钻铤和钻铤外流体界面(r=rb)处：

(12)

钻铤外流体和地层界面(r=ra)处：

(13)

将图1模型各部分声场的应力和位移表达式代入以上边界条件，得到矩阵方程：

H·O=Q,

(14)

其中向量O和Q的具体表达式如下：

Q={0,0,0,0,0,ps,0,0,0,0,0,0}T.

矩阵H中各元素的详细形式见附录B.随钻声场中模式波的频散特征由H确定，令其行列式G为零，得到各模式波的频散方程(宋云红等，2020)：

G(k,ω)=detH(k,ω)=0,

(15)

通过求解该方程得到各模式波波数k在不同频率下的根，然后利用式(16)计算各模式波的相速度频散曲线.式(16)为：

(16)

参考Tang和Cheng(2004)等文献，钻铤外环上的声压波形可由式(17)计算：

(17)

其中S(ω)为声源函数.

1.4 分波理论基础

由于VTI地层随钻声场中钻铤模式波与地层模式波耦合复杂，难以从全波波形中区分，因此采用分波计算的方法来分析钻铤与地层各向异性对地层模式波的影响.分波计算的理论是根据柯西积分定理，将复波数域中沿实轴的积分转化为实质支点的割线积分与极点的留数之和，波数域中这些奇异点的贡献构成了波场中的各种成分，使其呈现出不同的传播形态.此方法的优点是可以有效地提取全波中的各个成分，有利于分析各成分波的传播特征.相比之下，经典的实轴积分只能计算全波波形，需进一步处理才能分析各成分波的特征；但实轴积分计算的全波可以用来验证分波的计算结果.

计算中会涉及到径向虚波数的平方，使得随钻声场函数产生多值性，径向波数为零的奇异点被称作式(17)中被积函数的支点，其中部分支点对应的径向虚波数的正负值会产生多值性，这样的支点被定义为实质支点，反之被定义为非实质支点.针对随钻声波测井模型，郑晓波(2017)指出当声源不在最内层介质时，有且只有最外层介质的支点才是整个随钻声场的实质支点.针对横向各向同性地层，何晓(2010)证明：井孔声场函数的实质支点的个数和分布受地层各向异性程度的差异的影响，但其中地层的纵横波支点始终是井孔声场函数的实质支点，其具体表达式为：

(18)

另外两个支点的表达式如下：

(19)

(20)

式中系数a、b、c由弹性常数确定，具体表达式为：

He和Hu(2009)对VTI地层各向异性与井孔声场函数的支点数的关系给出了严格的数学证明；支点的个数及分布与Thomsen(1986)参数有关：

(21)

共分为3种情况：当参数满足δ>ε+c44/2c33时，±kh和±kl均是井孔声场函数的实质支点，此时VTI地层的单值定义如图2所示；另外两种地层参数情况下井孔声场函数的单值定义参考何晓(2010)的论文.

图2 弹性TI地层井孔声场响应函数的单值定义(δ>ε+c44/2c33)Fig.2 Singled-value definition of the borehole acoustic response function for a TI elastic formation with δ>ε+c44/2c33

按图2所示，式(17)中计算全波的实轴积分路径可转化为复波数域中(半径很大的)半圆弧与绕过实质支点的路径之和；根据约当引理，圆弧段的贡献为零，只有竖直围道的路径对积分有贡献，对应着沿井壁传播的地层体波.以地层纵波支点为例，地层纵波支点对声场的贡献的表达式为：

(22)

其中，“R”代表积分路径在支点右侧，“L”代表积分路径在支点左侧，式中对波数k的积分则对应着体波的激发强度(本文的激发强度为井中接收声压信号的强度).

复波数域中极点的留数对声场的贡献对应着全波中的模式波，其表达式为：

(23)

其中，Ni(k,ω)则是H矩阵中第i列被声源向量Q取代后的矩阵行列式的值，式(23)中极点kpole为一阶极点.模式波的激发强度可以表示为：

(24)

2 数值计算结果及分析

本节基于上述随钻声场理论及分波计算方法，计算了横向各向同性地层中，随钻单极、偶极和四极子声源激发的各模式波的频散曲线和激发强度曲线，并在此基础上给出各极点和支点对应的时域波形.计算时所选地层及井孔参数由表1和表2给出，本文计算所用声源为余弦包络脉冲函数(郑晓波，2017).

表1 地层参数Table 1 Formation parameters

表2 钻铤和流体参数Table 2 Drill collar and fluid parameters

2.1 随钻单极声场中的分波特征

图3a、b分别给出了横向各向同性地层随钻单极声场中各模式波的频散和激发强度曲线.图3a中实线对应的是随钻情况下的频散曲线，作为对比也给出了(无钻铤)电缆条件下的频散曲线(点划线)，图中Vp、Vsh、Vsv、Vf分别表示地层纵波，地层水平向横波，地层垂向横波以及井孔流体声波的波速.如图所示，标识符号“P-R”表示的地层伪瑞利波的低频截止速度接近地层垂向横波速度Vsv，并且，钻铤的引入使声场中增加了钻铤的拉伸型模式波(C1、C2)和内斯通利波(IN-ST).另一个重要现象是钻铤对伪瑞利波和斯通利波(ST)的频散特征产生了较大的影响，使得斯通利波的相速度整体下降；而伪瑞利波的截止频率往高频移动，其频散变化与电缆测井时相比大幅度减弱，有利于地层横波速度的测量.从图3b中各极点和支点激发强度曲线中可看出，在随钻声场中，地层体波的信号十分微弱，是因为其对应的支点激发强度(图中虚线)与模式波激发强度(实线)相比，要小1～2个数量级.以下重点研究随钻条件下井孔模式波的传播特征.

根据式(22)、(23)可计算出各个极点和支点对随钻声场的贡献，给出各分波的波形.图4给出了计算所得的随钻单极子声场中各支点(体波)和极点(模式波)的分波波形，作为对比也给出了实轴积分法计算的全波波形(波形计算的声源主频为8 kHz，源距为3 m).图中第1～7道给出了这些分波的归一化波形，分波的标识符号(图左端)与图3相同，分波相对于全波的实际波幅比例在图右端给出.计算结果表明，全波中的首波是由C1钻铤模式波与纵波支点kP贡献的分波组成；由于纵波支点的激发强度远小于C1钻铤模式波，如图3b所示，地层纵波被钻铤模式波所掩盖；地层横波波群组成复杂，按照到时，包括伪瑞利波和支点ks、kh的贡献，其中伪瑞利波占主要成分，因此伪瑞利波在提取地层垂向横波速度时影响很大；斯通利波与内斯通利波混叠，组成的波包在全波波形中幅度最大，到达最晚；C2钻铤模式波截止频率较高，且在8 kHz的主频段激发强度小，因此其产生的分波没有给出；图中第8道为实轴积分法计算的全波波形(红色虚线)与各分波叠加后的波形(黑色实线)的对比，二者无论在相位还是幅度上均十分吻合，验证了分波计算结果的正确性.

图3 (a)VTI地层单极随钻声场模式波频散曲线；(b)声场响应函数中的极点和支点对应的流体声压激发强度Fig.3 (a) LWD monopole wave velocity dispersion curves for a VTI formation; (b) Pressure excitation spectra for poles and branch points of the borehole acoustic response function

图4 VTI地层中单极随钻声场的极点和支点对应的分波Fig.4 Component waves from poles and branch points of LWD monopole wavefield in a fast VTI formation

2.2 随钻偶极声场中的分波特征

图5a给出了横向各向同性地层随钻偶极声场中各模式波的频散曲线(实线)，并与电缆测井中的结果进行对比(点划线).钻铤的存在使地层弯曲波(标识为Flexural的实线)频散曲线产生崎变，高频段随频率增加，但其低频截止仍是地层垂向横波速度；与电缆测井相比，随钻偶极声场中增加了多个模式波，包括频散特征受井外地层影响较大的界面波(Surface)以及钻铤模式波(D1-D4).由于钻铤与井孔的相互作用，界面波和弯曲波的中高频部分从形态上看像是无钻铤时的地层弯曲波与无井孔时的钻铤弯曲波发生了交换.钻铤模式波(D1-D4)在高频时相互耦合的现象复杂，但在低频段对偶极声波测井的影响可以忽略.

考察随钻偶极子声场中各个分波的激发响应特征.如图5b所示，3个支点的分波激发强度远小于各模式波的强度，说明随钻偶极声场中地层体波信号依然十分微弱，常被其他模式波所掩盖.如图6给出了图5对应的随钻偶极子声场的分波波形.分波的标识符号(图左端)与图5相同，分波相对于全波的实际波幅比例在图右端给出.地层纵波最先到达，但因为其幅度极小，无法在全波中观测到；地层横波波群构成复杂，由支点ks、kh贡献的分波以及地层弯曲波和钻铤模式波(D1-D4)共同组成(D3和D4在主频段幅度太小，未画出)，其中地层弯曲波占主要成分；界面波幅度大速度慢，在其之后是偶极换能器激发的六级声场，值得指出的是分波计算时仅考虑了偶极声场，因此通过分波合成得到的全波声场中不包括其余的高阶声场(郑晓波，2012)；钻铤模式波(D1-D4)以及支点kh贡献的分波呈震荡状，但其前后部分叠加后相互抵消.各分波叠加得到合成全波(第9道黑实线)与实轴积分法得到的全波(红线)十分吻合.

图5 (a) VTI地层偶极随钻声场模式波频散曲线； (b) 声场响应函数中的极点和支点对应的流体声压激发强度Fig.5 (a) LWD dipole wave velocity dispersion curves for a VTI formation; (b) Pressure excitation spectra for poles and branch points of the borehole acoustic response function

图6 VTI地层中偶极随钻声场的极点和支点对应的分波Fig.6 Component waves from poles and branch points of LWD dipole wavefield in a fast VTI formation

2.3 随钻四极声场中的分波特征

随钻横波测量采用四极声源(唐晓明和郑传汉，2004)，且大多针对慢速地层的情况.图7给出了快速(图7a、b)和慢速(图7c、d)横向各向同性地层随钻四极声场中各模式波的频散曲线和分波的激发强度曲线，图7a、c中实线是随钻情况下的频散曲线，点划线是电缆情况下的频散曲线.快速地层中，测井频带范围内总共出现了四种模式波，分别是地层模式波(Q1、Q2)和钻铤模式波(C1、C2)；地层模式波(Q1、Q2)的低频截止速度接近地层垂向横波速度.相对于电缆情况而言，钻铤使地层一阶模式波Q1的频散增强，不利于地层横波速度的测量，但二阶模式波Q2的频散变化却大为降低.由于受频散效应的影响较小，快速地层中采用二阶四极子模式波测量垂向地层横波速度较为有利，这与Tang等(2003)对快速各向同性地层四极子声波研究得出的结论一致.相比于快速地层，慢速地层中只有单一的地层模式波(Q1)，其低频部分趋近于地层垂向横波速度，且频散程度(相对于快速地层)较弱.与快速地层类似，慢速地层中分波激发响应中支点的强度远小于极点的强度，所以测量中观察不到随钻四极子声场中的地层体波.

图7 VTI地层四极子随钻声场模式波频散曲线和声场响应函数中的极点和支点对应的流体声压激发强度(a) 快速地层频散曲线； (b) 快速地层激发强度; (c) 慢速地层频散曲线； (d) 慢速地层激发强度.Fig.7 LWD quadrupole wave velocity dispersion curves in VTI formations and pressure excitation spectra for the poles and branch points of the borehole acoustic response function(a) Fast formation dispersion curves; (b) Fast formation excitation; (c) Slow formation dispersion curves; (d) Slow formation excitation.

图8a、b分别给出了快速和慢速地层中的分波计算结果.快速地层中，全波列首波是由Q2四极子模式波与来自支点ks和kh的分波共同组成，其中Q2模式波占主要成分，这与Tang等(2003)认识一致.采用Q2模式而非Q1模式来提取沿井轴方向的横波速度，是因为Q1模式波的频散大，高频时其速度远小于地层垂向横波速度，这也是全波列中Q2模式波比Q1模式波到时超前的原因.虽然Q1模式波的激发强度大于Q2模式波，但两者无量级上的差异，因此利用Q2模式波来测量地层垂向横波速度是可行的.在随钻测井中，钻铤模式波可能严重干扰Q2模式波的波形，对此常通过在钻铤上周期性刻槽来削减钻铤模式波的干扰.除此之外，钻铤模式波频散特征与钻铤尺寸密切相关：当钻铤外径不变，内径减小时，频散曲线向高频偏移，反之则向低频偏移，因此在实际仪器设计中，也应选择合适的仪器结构来尽量避开钻铤模式波的干扰；慢速地层中，Q1四极子模式波振幅最大，因此可利用该模式波来测量地层垂向横波速度.与单极和偶极的分析一样，将各分波叠加合成后的波形(黑线)与实轴积分计算结果(红色虚线)吻合.全波中的高频震荡是一阶钻铤模式波C1的干扰，测井时采用该钻铤模式波截止频率以下的频段便可消除这一干扰(唐晓明和郑传汉, 2004).

图8 VTI地层中四极随钻声场的极点和支点对应的分波(a) 快速地层； (b) 慢速地层.Fig.8 Component waves from poles and branch points of LWD quadrupole wavefield in VTI formation(a) Fast formation; (b) Slow formation.

2.4 VTI地层随钻四极声波低频波速的讨论

虽然四极子声波适宜于随钻横波测量(唐晓明和郑传汉，2004)，但对于VTI地层各向异性多种复杂变化情况对测量的影响加以讨论很有必要.Thomsen(1986)对实际地层岩石的TI各向异性数据做了大量调研，总结出了56组典型的TI岩石参数.He和Hu(2009)根据地层参数差异对这56组TI岩石进行分类，以式(21)的Thomsen参数及其组合来衡量，发现有接近90%的岩石各向异性参数满足δ≤ε+c44/2c33；考察这一普适条件对随钻四极子横波测量的影响，参照Thomsen(1986)论文，选取其中6组典型TI地层岩石参数(见表3快地层1～3、慢地层1～3)，依次是Berea sandstone、Oil Shale、Green River shale、Aluminum-lucite composite、Pierre shale以及Gas sand-water sand.

表3 地层各向异性参数Table 3 Anisotropic parameters of formations

计算结果如图9所示，实线为四极子模式波(Q1、Q2)的相速度，点划线为群速度频散曲线.图9a的快速地层中，4条水平虚线自上而下依次代表Berea sandstone、Oil Shale以及Green River shale

图9 满足条件δ≤ε+c44/2c33 的VTI地层随钻四极子模式波频散曲线(a) 快速地层; (b) 慢速地层.Fig.9 LWD quadrupole wave dispersion curves in VTI formations with δ≤ε+c44/2c33(a) Fast formation； (b) Slow formation.

岩石的垂直向横波速度，图9b的慢速地层中，4条水平虚线自上而下依次代表井内流体速度，Aluminum-lucite composite、Pierre shale以及Gas sand-water sand岩石的垂直向横波速度；结果表明在满足条件δ≤ε+c44/2c33的快速和慢速TI地层岩石中，地层四极子模式波的低频截止速度都接近于真实地层的垂向横波速度；并且，与图7类似，快速VTI地层的Q1四极子模式波为强频散模式，而Q2四极子模式波的频散相对较弱.

另外，He和Hu(2009)对Thomsen数据的统计也表明，约有10%岩石的TI各向异性满足条件δ>ε+c44/2c33.针对这一类各向异性岩石计算随钻四极子模式波的频散曲线，参数分别选自Thomsen(1986)论文中典型的Mesaverde laminated siltstone和何晓(2010)论文，由表3快地层4和慢地层4给出.

图10分别给出了快速(图10a)和慢速(图10c)VTI地层四极子模式波的相速度(实线)和群速度频散曲线(点划线)，结果表明在满足δ>ε+c44/2c33的条件下，无论快速还是慢速VTI地层的四极子模式波截止频率处的波速均小于真实地层垂向横波速度，趋于实质支点kh对应的速度vh=ω/kh(蓝色虚线)，这与上述满足条件δ≤ε+c44/2c33的各向异性地层很不相同.图10b、d进一步给出了快速和慢速VTI地层的分波计算结果；图10b中自上而下分别是主频为8 kHz的Q2四极子模式波的分波，主频为3.5 kHz的Q1模式波分波，以及主频为8 kHz的全波.作为对比，全波波形中画出了Q2模式波分波(红色虚线)，说明此时首波仍是Q2四极子模式波；图中绿点代表地层垂向横波到时，可见快速地层中四极子模式波到时滞后于沿井孔传播的地层横波.图10d中主频为3kHz的地层四极子模式波Q1到时也同样出现滞后于地层横波的现象.值得指出的是，单极伪瑞利波和偶极弯曲波也有类似现象，但是在前文的算例中，Vh与Vsv接近，所以这一现象不太明显.上述现象与He和Hu(2009)对裸眼井的理论分析结果一致，其产生的原因是，在参数满足条件δ>ε+c44/2c33的TI地层中，实质支点ks右侧将会出现新的实质支点kh，使四极子模式波的截止频率波速不是趋于地层垂向横波速度，而是趋于支点kh对应的速度vh.测量该速度是有意义的，因为根据式(19)，vh是地层各向异性参数的组合，可以用所测波速结合这些参数来求取地层垂向横波速度.

图10 满足条件δ>ε+c44/2c33的VTI地层随钻四极子模式波频散曲线和分波分析(a) 快速地层频散曲线; (b) 快速地层分波分析； (c) 慢速地层频散曲线; (d) 慢速地层分波分析.Fig.10 LWD quadrupole wave dispersion and component wave analysis in VTI formations with δ≤ε+c44/2c33(a) Fast formation dispersion； (b) Fast formation complement wave； (c) Slow formation dispersion； (d) Slow formation complement wave.

3 结论

本文通过对横向各向同性地层随钻声场的分波分析和计算，将声场全波分解成来自波数域中的奇异点(包括支点和极点)的分波的贡献，对比研究了单极、偶极和四极子声场模式波的传播特征，得出以下结论：

(1)与电缆测井结果相比，快速地层中钻铤的存在使得斯通利波和地层一阶四极子模式波(Q1)的频散变化增加，不利于地层横波速度的测量；钻铤和地层的弯曲波频散特征和形态也发生了较大变化；伪瑞利波和地层二阶四极子模式波(Q2)的频散则大幅度减弱，有利于地层横波速度的测量.

(2)快速VTI地层随钻多极子声场的横波波群构成复杂，单极的横波波群由伪瑞利波与支点ks和kh贡献的分波共同组成，其中伪瑞利波占主要成分；偶极的横波波群由弯曲波和一阶钻铤模式波(D1)与支点ks和kh贡献的分波共同组成，其中弯曲波占主要成分；四极的横波波群由地层二阶四极子模式波(Q2)与支点ks和kh贡献的分波共同组成，其中地层二阶四极子模式波(Q2)占主要成分.慢速VTI地层中，横波波群由地层四极子模式波和支点ks贡献的分波共同组成，其中地层四极子模式波占主要成分.因此，多极子声场中的伪瑞利波、弯曲波和四极子螺旋波可用于测量VTI地层的垂直向横波速度.

(3)在各向异性参数满足条件δ>ε+c44/2c33的VTI地层中，伪瑞利波、弯曲波以及四极子螺旋波的截止频率处的波速均小于垂直向地层横波速度，将这一波速与地层各向异性参数结合可以获得垂直向地层横波速度.

附录A

井孔外VTI地层径向位移，径向、环向及轴向应力分量表达式为：

附录B

H矩阵元素为(所有未列出的矩阵元素都为0)：

H21=ρfω2In(frc)，

H68=ρfω2In(frb)，

H69=ρfω2Kn(frb)，

H10,8=ρfω2In(fra)，

H10,9=ρfω2Kn(fra)，

+qpKn+1(qpra))，

+qsvKn+1(qsvra))，

-qpKn+1(qpra))，

+qsvKn+1(qsvra))，