2021, Vol. 64
2. 中国科学院地质与地球物理研究所, 中国科学院油气资源研究院重点实验室, 北京 100029;
3. 非对称性弹性波动方程联合研究组, 北京 100029;
4. 高铁地震学联合研究组, 北京 100029
2. Key Laboratory of Petroleum Resource Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. The Joint Research Group of Asymmetric Elastic Wave Equations, Beijing 100029, China;
4. The Joint Research Group of High-Speed Rail Seismology, Beijing 100029, China
根据经典连续介质力学理论，可推导传统弹性波动方程，波动方程中不包含独立的旋转自由项以及介质特征尺度参数.经典连续介质力学理论中，介质材料点的应变能密度函数只依赖于经典应变，介质被建模为连续质量体而不是离散的点.对于均匀各向同性介质，没有内部结构，内部附加自由度以及内部特征尺度，应用动量矩守恒定律，必然得到应力张量是对称的结论(赵亚溥，2016).但是，在广义连续介质力学理论中，当考虑偶应力或体力偶时，切应力互等定理不成立了，这导致应力张量的不对称性，因此，基于广义连续介质力学理论可推导非对称弹性波动方程.应用非对称弹性波动方程进行数值模拟时，在合成记录中可以观察到弹性波的传播出现了尺度效应，即地震记录中出现了新的波场成分，且随介质特征尺度的变化而变化.非对称弹性波动方程相比传统弹性波动方程，会增加独立的自由项，该自由项包含介质特征尺度参数，以描述由于偶应力的引入，在颗粒介质或者有裂隙的连续体内由于力学不对称性所产生的一种非均匀性扰动的传播.尽管传统弹性波动方程己为人类做出了巨大的贡献，可是任何从事实际工作的，油公司和服务公司处理解释人员，深知其巨大的不可知因素；多年来，学者们一直都在深化探索.经典连续介质理论和广义连续介质理论，对称弹性波动方程和非对称弹性波动方程，不是对立的关系，而是包含和相互补充的关系，结合在一起，在描述介质内多尺度的微结构相互作用对地震记录的影响规律时，可能是一种可以考虑的方法和途径.相比于经典连续介质力学理论，广义连续介质力学理论(Eringen, 1966, 1990；Mindlin and Eshel, 1968；Tomar and Gogna, 1992；Aifantis, 1999, 2011；Polizzotto，2013；Auffray et al., 2013；Singh，2002；Tomar and Garg, 2005；Madeo et al., 2015)没有将介质视为理想的连续统(ideal continuum)，而是认为介质是具有复杂内部微结构的.因此，应用动量矩守恒定律，可以得到应力张量是非对称的结论.严格上说，无论是对于自然介质还是人造介质，都是具有不同尺度的复杂内部微结构的，理想连续统的介质是不存在的.为了描述介质内部复杂的微结构相互作用，其中一个方法是在经典连续介质力学的框架下，将组成介质的微元体分别建模，随之而来的是计算复杂度的提升以及计算资源要求的增加，除此之外，在对非常复杂的构造进行数值模拟时可能还会出现病态问题(Zhu et al., 2020).另外一个方法就是基于广义连续介质力学的理论描述介质内部微结构相互作用.广义连续介质力学理论始于Cosserat提出的偶应力理论(Cosserat and Cosserat, 1909).Grioli(1960)提出一种考虑非对称应力的线弹性理论，并认为可以解决经典弹性理论遇到的一些问题.Toupin(1962)，Mindlin和Tiersten(1962)对Cosserat偶应力理论进行了约束，假设微元体相对转动为零，建立了线弹性偶应力理论.Toupin(1964)将应变的梯度引入到介质的应变能密度函数中，提出了应变梯度理论.Mindlin(1965)建立了二阶应变梯度理论，该理论包含所有高阶变形量的影响，是最一般的理论(Mindlin and Eshel, 1968).除了具有2个Lame常数外，还引入了5个非相互独立的高阶介质特征尺度参数，但是这5个高阶介质特征尺度参数很难进行实验标定，因此，为了进行实际应用，常对该理论进行简化或者近似，包括偶应力理论，应变梯度理论都可以看成是该理论的简化形式.Yang等(2002)通过引入额外的高阶力矩平衡理论，即假设连续体内偶力矩为零，提出修正偶应力理论，在该理论中，仅包含一个介质特征尺度参数.Chen等(2011)将偶应力理论推广到各向异性介质.Lam等(2003)仿照Yang等(2002)的思路，排除了曲率张量反对称部分对应变能密度函数的影响，提出了修正的应变梯度理论，在其理论中，介质应变能密度函数依赖于传统应变张量，膨胀梯度矢量，拉伸梯度张量以及曲率张量的对称部分.Aifantis(1999)基于非局部弹性理论的思路，考虑了一点的应力状态受整个介质中诸点应力影响的情况，在本构关系中引入应变的拉普拉斯算子，提出了一种单一介质特征尺度参数的应变梯度理论，可视为应变梯度理论的一种特殊形式(Askes and Aifantis, 2011).Auffray(2013)将应变梯度理论推广到各向异性介质，提出了一种适用于各向异性介质的应变梯度理论.以上这些理论与Mindlin(1965)理论的区别在于采用不同的平衡条件或者主动简化某些应变梯度量的影响，并减少了高阶介质特征尺度参数的数量，降低了在描述介质不同尺度的微结构相互作用时的能力，但是却提高了可实现性.
在广义连续介质力学理论框架下，由于考虑介质内微孔缝隙结构的相互作用导致的不均匀性，弹性波的传播出现了尺度效应.广义连续介质力学理论通过引入应力或者应变的高阶空间导数项丰富经典连续介质力学理论的内容，而这些应力与应变的高阶导数，通常附带有不同的介质特征尺度参数.这些特征尺度参数与介质微孔缝隙特征尺度有密切的关系(Baant，2002；Baant and Pang, 2007)，且与介质内部的几何结构(黄文雄和徐可，2014)以及应力扰动效应也有关.为了进一步研究分析介质内微孔缝隙结构相互作用对地震记录的影响，我们基于修正偶应力理论(Yang et al., 2002)，推导包含介质特征尺度参数的非对称弹性波动方程以及对应的非对称SH型横波波动方程，并在三层煤层模型上进行数值模拟，将检波器分别设置在地表和煤层中线，通过改变介质内微孔缝隙特征尺度参数值，合成地震记录，研究分析弹性波传播中，由介质内微孔缝隙结构相互作用所导致的尺度效应，对地震记录的影响及规律，并给出结论.
1 非对称弹性波动方程由于假设组成介质的微元体可以承受力偶的作用，则应力张量是非对称的.基于修正偶应力理论(Yang et al., 2002)，可推导包含介质微孔缝隙特征尺度的非对称弹性波动方程以及对应的SH型横波波动方程.
偶应力理论下，应变能密度函数不仅与传统应变张量有关，也与旋转矢量的梯度有关，定义旋转矢量的梯度为曲率张量，表征扭转与弯曲形变.Yang等(2002)给出了修正偶应力理论下的应变能密度表达式：
|
(1) |
其中，εij为传统应变张量，
分别与应变张量和曲率张量功共轭的传统应力张量σij和偶应力张量μij，可通过本构关系给出：
|
(2) |
小变形情况下，对于各向同性介质，应用广义胡克定律以及功共轭条件(见公式(2))，可得修正偶应力理论下的本构关系：
|
(3) |
其中，η=l2μ，是反映介质微旋转变形特性的参数，也是介质内部孔缝隙尺度效应特征的表征，l是一个综合参数，与介质的变形局部化现象(比如剪切带或剪切表面层的变形局部化)，组成介质的材料平均直径和几何形状以及应力扰动效应有关，且受其加成作用，如果不考虑组成介质的材料的几何特征或者微结构导致的应力扰动效应以及介质的变形局部化问题，其可直接视为为介质内微孔缝隙特征尺度.δij为Kronecker delta符号.
对于一个体积为Va，边界为Sa的连续体，修正偶应力理论的虚功原理可表示为(Yang et al., 2002)
|
(4) |
其中，w是单位体积的应变能密度，δ为变分符号，Fi, Mi分别是体力和体力偶.
根据公式(4)的虚功原理，则可得修正偶应力理论下的平衡方程：
|
(5) |
边界条件为
|
(6) |
|
(7) |
其中，nn，(I－nn)分别代表法向和切向方向.p(n)为表面应力矢量，m(n)为表面力偶矢量.
一般认为，真实边界上自然边界条件处处为零.
结合公式(3)和公式(5)，可推导得出修正偶应力非对称弹性波动方程：
|
(8) |
如公式(8)所示，基于广义连续介质力学理论推导的非对称弹性波动方程，包含独立的自由项以及介质特征尺度参数，描述引入偶应力后，由力学不对称特征导致的旋转运动，所产生的位移扰动传播.因为考虑η的作用，地震波传播出现明显的尺度效应，实际上，这是介质内微孔缝隙结构相互作用的宏观体现.当η=0，公式(8)即为传统弹性波动方程的表达式：
|
(9) |
即非对称弹性波动方程包含传统弹性波动方程.
以下为进一步推导修正偶应力非对称SH型横波波动方程.
令ux=uz=0，则基于公式(8)，可得到非对称SH型横波波动方程：
|
(10) |
如果不考虑微孔缝隙结构相互作用，即设置η=0，可得到传统SH型横波波动方程：
|
(11) |
本文将应用修正偶应力非对称弹性波动方程和非对称SH型横波波动方程，对不同介质孔缝隙特征尺度下的煤层模型进行数值模拟，研究分析弹性波传播的尺度效应对地震记录的影响和规律.值得注意的是，尺度效应是一个相对概念，对于岩土介质，尺度效应可出现在宏观可视的范围内.这是因为尺度效应由介质材料的微缺陷/微结构引起，而介质材料的微缺陷/微结构本来就是一个相对概念.对于岩土介质，其微缺陷/微结构通常指的是微孔缝隙，量级一般在百微米级以上(对于煤层，微纳米量级居多)，另外，尺度效应也受到Lame问题下的位移场/旋转场的空间变化梯度影响加权.对于给定的一个较小的介质特征尺度参数值，尺度效应虽然很小，但是也存在.比如，对于给定的微纳米量级的介质特征尺度，尺度效应导致的新的波场成分的数量级至少要降低3个量级，因此，需要采用较高精度的数值模拟算法以及较好的边界条件，降低数值频散噪声及边界条件所带来的干扰，以更好地提取尺度效应所导致的新的波场成分.除此之外，可以考虑介质的变形局部化现象，组成介质的材料平均直径和几何形状以及应力扰动效应等因素对介质特征尺度的加成作用，将介质特征尺度视为上述因素对介质内微孔缝隙特征尺度加成的结果(王之洋等(2020)通过对高铁实际数据的分析，提取了加成系数)，以增强尺度效应.
2 数值模拟设计一个三层煤层模型并进行数值模拟.如图 1所示，模型大小设置为nx=nz=200，网格间距为1 m×1 m.模型中间层为厚度为10 m的煤层，P波速度为2.4 km·s-1，S波速度为1.1 km·s-1、密度为1400 kg·m-3；模型上下两层分别为厚度是95 m的高速围岩层，P波速度为4.2 km·s-1，S波速度为2.0 km·s-1、密度为2600 kg·m-3.
|
图 1 煤层模型 Fig. 1 Coal-seam three-layer model |
首先，分别使用传统SH型横波波动方程(公式(11))和非对称SH型横波波动方程(公式(10))，进行二维Love型槽波数值模拟.根据Love型槽波的激发和传播特性，震源采用主频为100 Hz的Ricker子波，震源位置在(x, z)=(100 m, 100 m)，采样时间间隔为0.0001 s.检波器排列在z=100 m这条测线.图 2和图 3分别是使用传统SH型横波波动方程和非对称SH型横波波动方程获得的二维Love型槽波的波场快照和地震记录.
|
图 2 介质内微孔缝隙特征尺度为0时使用不同弹性波动方程获得的波场快照(y分量)及地震记录(y分量) (a) 使用传统SH型横波波动方程获得的波场快照；(b) 使用非对称SH型横波波动方程获得的波场快照(介质内微孔缝隙特征尺度为0)；(c) 图(a)和(b)之间的差异；(d) 使用传统SH型横波波动方程获得的地震记录；(e) 使用非对称SH型横波波动方程获得的地震记录(介质内微孔缝隙特征尺度为0)；(f) 图(d)和(e)之间的差异. Fig. 2 Snapshot (y component) of wave field and shot records (y component) using different elastic wave equations when the characteristic scale of the micropore in the media is 0 (a) Snapshot of wave field using the conventional SH-wave equation; (b) Snapshot of wave field using asymmetric SH-wave equation (the characteristic scale of the micropore in the media is 0); (c) The difference between Fig. 2 (a) and (b); (d) Shot records using the conventional SH-wave equation; (e) Shot records using asymmetric SH-wave equation (the characteristic scale of the micropore in the media is 0); (f) The difference between Fig. 2 (d) and (e). |
|
图 3 使用不同弹性波动方程获得的波场快照(y分量)及地震记录(y分量) (a) 使用传统SH型横波波动方程获得的波场快照；(b) 使用非对称SH型横波波动方程获得的波场快照(介质内微孔缝隙特征尺度为300 nm)；(c) 图(a)和(b)之间的差异；(d) 使用传统SH型横波波动方程获得的地震记录；(e) 使用非对称SH型横波波动方程获得的地震记录(介质内微孔缝隙特征尺度为300 nm)；(f) 图(d)和(e)之间的差异. Fig. 3 Snapshot (y component) of wave field and shot records (y component) using different elastic wave equations (a) Snapshot of wave field using the conventional SH-wave equation; (b) Snapshot of wave field using asymmetric SH-wave equation (the characteristic scale of the micropore in the media is 300 nm); (c) The difference between Fig. 3 (a) and (b); (d) Shot records using the conventional SH-wave equation; (e) Shot records using asymmetric SH-wave equation (the characteristic scale of the micropore in the media is 300nm); (f) The difference between Fig. 3 (d) and (e). |
如图 2可以发现，当设置介质内微孔缝隙特征尺度为0时，非对称SH型横波波动方程和传统SH型横波波动方程的数值模拟结果，不论是炮记录还是波场快照，两者都是相同的，如果将两个方程的合成记录和波场快照简单相减，如图 2(c, f)所示，两者完全相同.同时，在记录上可以清晰观察到反射波，折射波和透射波以及煤层中的Love型槽波，槽波的频散特征十分明显.进一步，将炮记录进行了频散分析，两个不同的方程得到的结果完全匹配，这说明非对称SH型横波波动方程包含传统SH型横波波动方程.
如果考虑微孔缝隙结构相互作用，设置介质内微孔缝隙特征尺度为300 nm，则使用两个方程得到合成记录是不完全一样的，无论是炮记录还是波场快照，都出现了新的波场成分，如图 3(c, f)所示.设置了一个较小的介质内微孔缝隙特征尺度值，传统方程和非对称方程在槽波模拟时得到的合成记录仍然有明显的差别，这表明，由介质内微孔缝隙结构相互作用所导致的尺度效应是可以被考虑的.
在进行比较之前，已经进行了数值频散压制以及边界处理，所以结合理论分析以及图 2和图 3的比较，可以认为，合成地震记录的改变是由弹性波传播的尺度效应产生的.
为进一步分析尺度效应对地震记录的影响，通过设置不同的介质内微孔缝隙特征尺度值，并使用传统和非对称弹性波动方程，继续在图 1所示的煤层模型上进行测试，检波器置于地面.震源位置设置为(x, z)=(100 m, 1 m)，震源采用Ricker子波，主频为100 Hz，检波器排列在z=1 m这条测线上.
首先，分别使用传统弹性波动方程和非对称弹性波动方程得到合成地震记录，如图 4和图 5所示.图 4(a, d)分别是使用传统弹性波动方程得到的合成地震记录x分量和z分量.图 4(b, e)分别是使用非对称弹性波动方程得到的合成地震记录x分量和z分量，我们将围岩微孔缝隙特征尺度和煤层微孔缝隙特征尺度都设置为零.图 4(c, f)分别是图 4(a, b)和图 4(d, e)相减得到的结果.由图 4的比较可知，当设置介质内微孔缝隙特征尺度为零时，使用非对称弹性波动方程和传统弹性波动方程，合成的地震记录是完全相同的，即不考虑介质内微孔缝隙结构的相互作用导致的不均匀性时，非对称弹性波动方程就退化为传统弹性波动方程.
|
图 4 介质内微孔缝隙特征尺度为0时使用不同弹性波动方程获得的地震记录 (a) 传统弹性波动方程(x分量)；(b) 非对称弹性波动方程(围岩的微孔缝隙特征尺度为0, 煤层的微孔缝隙特征尺度为0)(x分量)；(c) 图(a)和(b)之间的差异；(d) 传统弹性波动方程(z分量)；(e) 非对称弹性波动方程(围岩的微孔缝隙特征尺度为0, 煤层的微孔缝隙特征尺度为0)(z分量)；(f) 图(d)和(e)之间的差异. Fig. 4 Shot records using different elastic wave equations when the characteristic scale of the micropore in the media is 0 (a) The conventional elastic wave equations (x component); (b) Asymmetric elastic wave equations (the characteristic scale of the micropore in the surrounding rock is 0, the characteristic scale of the micropore in the coal seam is 0) (x component); (c) The difference between Fig. 4 (a) and (b); (d) The conventional elastic wave equations (z component); (e) Asymmetric elastic wave equations (the characteristic scale of the micropore in the surrounding rock is 0, the characteristic scale of the micropore in the coal seam is 0) (z component); (f) The difference between Fig. 4 (d) and (e). |
|
图 5 使用不同弹性波动方程获得的地震记录 (a) 传统弹性波动方程(x分量)；(b) 非对称弹性波动方程(围岩的微孔缝隙特征尺度为300 nm, 煤层的微孔缝隙特征尺度为300 nm)(x分量)；(c) 图(a)和(b)之间的差异；(d) 传统弹性波动方程(z分量)；(e) 非对称弹性波动方程(围岩的微孔缝隙特征尺度为300 nm, 煤层的微孔缝隙特征尺度为300 nm)(z分量)；(f) 图(d)和(e)之间的差异. Fig. 5 Shot records using different elastic wave equations (a) The conventional elastic wave equations (x component); (b) Asymmetric elastic wave equations (the characteristic scale of the micropore in the surrounding rock is 300 nm, the characteristic scale of the micropore in the coal seam is 300 nm) (x component); (c) The difference between Fig. 5 (a) and (b); (d) The conventional elastic wave equations (z component); (e) Asymmetric elastic wave equations (the characteristic scale of the micropore in the surrounding rock is 300 nm, the characteristic scale of the micropore in the coal seam is 300 nm) (z component); (f) The difference between Fig. 5 (d) and (e). |
当我们将围岩微孔缝隙特征尺度和煤层微孔缝隙特征尺度都设置为300 nm时，使用两个弹性波动方程合成的地震记录，不论是x分量还是z分量，都是有明显差异的，见图 5所示.应用非对称弹性波动方程进行数值模拟时，在合成记录中可以观察到弹性波的传播出现了明显的尺度效应，即地震波场由于考虑了介质内微孔缝隙相互作用而出现了新的成分.结合理论分析和试验认为，这些新的波场成分和物理频散有关，由于考虑了介质内的微结构相互作用，弹性波传播会更加复杂.
其次，研究分析不同介质内微孔缝隙特征尺度对地震记录的影响，以及地震记录对介质内微孔隙特征尺度变化的敏感性.图 6和图 7分别为使用非对称弹性波动方程，在不同介质内微孔缝隙特征尺度下获得的合成地震记录.其中，图 6(a，d)中，围岩和煤层的微孔缝隙特征尺度均设置为1000 nm.图 6(b，e)中，设置围岩微孔缝隙特征尺度为1000 nm，煤层微孔缝隙特征尺度为300 nm.图 6(c，f)分别是图 6(a，b)和图 6(d，e)相减得到的结果.由图 6可知，不同介质内微孔缝隙特征尺度对地震记录的影响是不一样的，如果将煤层和围岩设置不同的微孔缝隙特征尺度，地震记录将根据微孔缝隙特征尺度表现出不同的尺度效应.因此，不论是从x分量还是z分量，煤层的位置都能被清晰的定位出来，以与围岩作区分，见图 6(c，f).
|
图 6 不同介质内微孔缝隙特征尺度参数下使用非对称弹性波动方程获得的地震记录 (a) 围岩的微孔缝隙特征尺度为1000 nm, 煤层的微孔缝隙特征尺度为1000 nm(x分量)；(b) 围岩的微孔缝隙特征尺度为1000 nm, 煤层的微孔缝隙特征尺度为300 nm(x分量)；(c) 图(a)和(b)之间的差异；(d) 围岩的微孔缝隙特征尺度为1000 nm, 煤层的微孔缝隙特征尺度为1000 nm(z分量)；(e) 围岩的微孔缝隙特征尺度为1000 nm, 煤层的微孔缝隙特征尺度为300 nm (z分量)；(f) 图(d)和(e)之间的差异. Fig. 6 Shot records using asymmetric elastic wave equations with different characteristic scale of microporous of media (a) The characteristic scale of the micropore in the surrounding rock is 1000nm, the characteristic scale of the micropore in the coal seam is 1000 nm (x component); (b) The characteristic scale of the micropore in the surrounding rock is 1000 nm, the characteristic scale of the micropore in the coal seam is 300 nm (x component); (c) The difference between Fig. 6 (a) and (b); (d) The characteristic scale of the micropore in the surrounding rock is 1000nm, the characteristic scale of the micropore in the coal seam is 1000 nm (z component); (e) The characteristic scale of the micropore in the surrounding rock is 1000 nm, the characteristic scale of the micropore in the coal seam is 300 nm (z component); (f) The difference between Fig. 6 (d) and (e). |
|
图 7 介质内微孔缝隙特征尺度参数差异较小时使用非对称弹性波动方程获得的地震记录 (a) 围岩的微孔缝隙特征尺度为370 nm, 煤层的微孔缝隙特征尺度为300 nm(x分量)；(b) 围岩的微孔缝隙特征尺度为370 nm, 煤层的微孔缝隙特征尺度为370 nm(x分量)；(c) 图(a)和(b)之间的差异；(d) 围岩的微孔缝隙特征尺度为370 nm, 煤层的微孔缝隙特征尺度为300 nm(z分量)；(e) 围岩的微孔缝隙特征尺度为370 nm, 煤层的微孔缝隙特征尺度为370 nm(z分量)；(f) 图(d)和(e)之间的差异. Fig. 7 Shot records using asymmetric elastic wave equations when the difference in characteristic scale of microporous of media is small (a) The characteristic scale of the micropore in the surrounding rock is 370 nm, the characteristic scale of the micropore in the coal seam is 300 nm (x component); (b) the characteristic scale of the micropore in the surrounding rock is 370 nm, the characteristic scale of the micropore in the coal seam is 370 nm (x component); (c) The difference between Fig. 7 (a) and (b); (d) The characteristic scale of the micropore in the surrounding rock is 370 nm, the characteristic scale of the micropore in the coal seam is 300 nm (z component); (e) The characteristic scale of the micropore in the surrounding rock is 370 nm, the characteristic scale of the micropore in the coal seam is 370 nm (z component); (f) The difference between Fig. 7 (d) and (e). |
逐步缩小煤层和围岩的微孔缝隙特征尺度差异，直至差值只有70 nm.如图 7所示，在图 7(a，d)中，设置围岩微孔缝隙特征尺度为370 nm，煤层微孔缝隙特征尺度为300 nm.在图 7(b，e)中，围岩和煤层的微孔缝隙特征尺度均设置为370 nm；图 7(c，f)同样分别是图 7(a，b)和图 7(d，e)相减得到的结果.图 7的结果可以说明，即使两个层位的微孔缝隙特征尺度差值仅有70 nm，仍然能够在地震记录中被观察到，这证明了地震记录对介质微孔隙特征尺度参数值的变化，具有较高的敏感性，这对后续的微孔缝隙反演应用是非常有利的.
通过理论分析和数值模拟试验，可以清晰的认识到，在勘探地震的尺度内，介质内的微孔缝隙特征尺度对地震记录的影响，是可以被考虑的.地震记录对介质内微孔缝隙特征尺度具有较高的敏感性，可以在合成地震记录中观察到明显由尺度效应带来的波场变化.
3 结论尺度效应是广义连续介质力学理论框架下，考虑介质内微孔缝隙结构相互作用，所导致的介质非均匀性在弹性波传播时的宏观体现.以煤层模型为例，应用非对称弹性波动方程和非对称SH型横波波动方程进行数值模拟，设置不同的介质内微孔缝隙特征尺度参数值，研究分析弹性波传播的尺度效应.得出以下认识和结论：
(1) 非对称弹性波动方程模拟的弹性波传播表现出明显的尺度效应，并且具有不同微孔缝隙特征尺度参数值的介质具有不同的地震记录；
(2) 地震记录受介质的微孔缝隙特征尺度影响，多层介质的每层设置不同的介质内微孔缝隙特征尺度参数值，则可以在地震记录中定位具体的层位；
(3) 地震记录对介质微孔隙特征尺度具有较高的敏感性；
(4) 地震记录是多种因素共同影响的结果，介质内微孔缝隙尺度在广义连续介质力学理论框架下对地震记录的影响，是可以被考虑的.
附录基于我们已经在地球物理学报录用的论文《偶应力理论框架下的弹性波数值模拟与分析》，以及修正的偶应力理论(Yang et al., 2002)，可推导包含介质微孔缝隙特征尺度的非对称弹性波动方程.
1 平衡方程考虑一个体积为Va，边界为Sa的连续体.在该连续体中的两个微元体之间的相互作用通过二者的接触表面ds上(单位法向量为ni)的表面应力矢量pi(n)与表面力偶矢量mi(n)传递，通常使用应力张量σij和偶应力张量μij来表示：
|
(A1) |
|
(A2) |
应力张量σij和偶应力张量μij均满足线动量平衡方程(公式(A3))和角动量平衡方程(公式(A4))：
|
(A3) |
|
(A4) |
其中，Vi, Ωi分别为线动量和角动量.Fi, Mi分别是体力和体力偶.eijk为置换符号.rj是位矢, ui是位移矢量，wi是旋转矢量.
对公式(A3)和(A4)分别应用散度定理，由表面积分变换为体积分，并忽略公式(A4)中的高阶微量，可以分别得到微分形式的平衡方程：
|
(A5) |
|
(A6) |
其中，s，t=1, 2, 3.
在经典连续介质力学理论中，没有考虑到偶应力与体力偶的作用，也即μij=0, Mi=0，因此，可由公式(A6)，得到：
|
(A7) |
而在偶应力理论中，当考虑偶应力或体力偶时，切应力互等定理不成立了，这导致应力张量的不对称性，这时候，应力张量σij可以分解为对称应力张量τij与反对称应力张量rij之和：
|
(A8) |
如果将公式(A8)代入公式(A5)和(A6)，将应力张量σij用对称应力张量τij与反对称应力张量rij表示，则公式(A5)和(A6)，可以改写为
|
(A9) |
|
(A10) |
同时基于公式(A6)，建立反对称应力张量rij与偶应力张量μij的关系式：
|
(A11) |
进一步，可以将偶应力张量μij分解为偏斜部分mij和球量部分
|
(A12) |
将公式(A11)等号两边求散度，并将公式(A12)代入，可得：
|
(A13) |
将公式(A8)和公式(A13)一齐代入公式(A9)，可消除反对称应力张量rij，得到只包含对称应力张量τij以及偶应力张量的偏斜部分mij的平衡方程：
|
(A14) |
公式(A14)即为偶应力理论框架下的平衡方程，其构建了偶应力张量，应力张量，体力偶和体力的关系，如果去除了偶应力张量和体力偶，该平衡方程就为经典理论的平衡方程.
2 边界条件考虑不存在体力和体力偶的单位体积元能量守恒：
|
(A15) |
其中，n是方向矢量，δu为虚位移矢量，δω为虚旋转矢量，δw为能量改变量.
因为有：
|
(A16) |
其中，nn，(I-nn)分别代表法向和切向方向.
对公式(A16)等号右边第一项应用Hamiltonian算子，并考虑表面Sa是光滑的，则公式(A15)可写为：
|
(A17) |
由公式(A17)，可得到偶应力理论框架下的边界条件：
|
(A18) |
|
(A19) |
这里仅仅考虑微小变形|ui, j|≪1的运动学(几何方程).在笛卡尔坐标系中，参考构型中的两点A(位置为xi)和B(位置为xi+dxi)之间的相对位移为
|
(A20) |
位移梯度张量ui, j可分解为对称的无穷小应变张量εij和反对称的无穷小旋转张量ωij之和：
|
(A21) |
公式(A21)中，ui, ui0分别为A点和B点的位移矢量.
其中：
|
(A22) |
|
(A23) |
反对称的无穷小旋转张量ωij对应的旋转矢量ωi定义为
|
(A24) |
因此，旋转矢量ωi可以通过求取位移场旋度的二分之一得到.
在微小变形|ui, j|≪1假设下，经典连续介质力学理论下的旋转运动，是刚性旋转，其不造成形变(Aki and Richards, 2002)，位移梯度全部用对称的应变张量表示.偶应力理论下，引入了一种由偶应力导致的新的旋转运动，使微元线段dxi产生了扭转和弯曲形变，其由相邻两点A, B之间的相对转动dωi描述，并引入曲率张量χij描述这种形变：
|
(A25) |
这样，在微小变形中，连续体A点临近B点的形变总结起来有4种：随着A点的平动，相对A点的伸缩，相对A点的转动以及微元线段
通过构建曲率张量以及应变张量与位移的关系，我们可以推导得到几何方程，相比于经典理论的几何关系，偶应力理论框架下增加了曲率张量与位移之间的关系，几何方程见公式(A26)：
|
(A26) |
由虚功原理，在公式(A5)两端乘以虚位移δui，在公式(A6)两端乘以虚旋转δωi，并在整个体积Va上积分，可得：
|
(A27) |
|
(A28) |
对公式(A27)和(A28)应用链式求导法则以及散度定理，并相加，可得：
|
(A29) |
将公式(A1)和公式(A2)代入公式(A29)，并应用散度定理，将公式(A29)等号右边第一项的面积分变为体积分，可得：
|
(A30) |
将平衡方程(公式(A14))、几何方程(公式(A26))以及公式(A8)，代入公式(A30)并化简，可得功共轭关系：
|
(A31) |
对于各向同性介质，应用广义胡克定律以及功共轭条件(公式(A31))，可得偶应力理论框架下的本构方程：
|
(A32) |
其中，λ, μ为Lamé常数，η为反映介质旋转运动特性的参数，η=μl2，l为介质特征尺度参数，其是为了平衡无量纲的应变和有量纲的曲率张量(其量纲为长度的倒数)两项间的量纲，也为了描述偶应力张量和曲率张量之间的本构关系.
5 偶应力理论框架下的非对称弹性波动方程对于各项同性介质，将几何方程(公式(A26))和本构方程(公式(A32))以及平衡方程(公式(A14))联立，并忽略体力以及体力偶，可推导得到偶应力理论框架下的弹性波动方程：
|
(A33) |
方程分量形式表示如下：
|
(A34) |
考虑在二维各向同性介质中，上式可写为
|
(A35) |
在二维各向同性介质中，若震源沿y方向振动，可产生SH型横波，令ux=uz=0，则基于公式(A35)，可得到非对称SH型横波波动方程：
|
(A36) |
Aifantis E C. 1999. Strain gradient interpretation of size effects. International Journal of Fracture, 95(1-4): 299-314. |
Aifantis E C. 2011. On the gradient approach-Relation to Eringen's nonlocal theory. International Journal of Engineering Science, 49(12): 1367-1377. DOI:10.1016/j.ijengsci.2011.03.016 |
Aki K, Richards P G. 2002. Quantitative Seismology: Theory and Methods. 2nd ed. San Francisco: W. H. Freeman and Co.
|
Askes H, Aifantis E C. 2011. Gradient elasticity in statics and dynamics: An overview of formulations, length scale identification procedures, finite element implementations and new results. International Journal of Solids and Structures, 48(13): 1962-1990. DOI:10.1016/j.ijsolstr.2011.03.006 |
Auffray N, Le Quang H, He Q C. 2013. Matrix representations for 3D strain-gradient elasticity. Journal of the Mechanics and Physics of Solids, 61(5): 1202-1223. DOI:10.1016/j.jmps.2013.01.003 |
Baant Z P. 2002. Scaling of dislocation-based strain-gradient plasticity. Journal of the Mechanics and Physics of Solids, 50(3): 435-448. DOI:10.1016/S0022-5096(01)00082-5 |
Baant Z P, Pang S D. 2007. Activation energy based extreme value statistics and size effect in brittle and quasibrittle fracture. Journal of the Mechanics and Physics of Solids, 55(1): 91-131. DOI:10.1016/j.jmps.2006.05.007 |
Chen W J, Li L, Xu M. 2011. A modified couple stress model for bending analysis of composite laminated beams with first order shear deformation. Composite Structures, 93(11): 2723-2732. DOI:10.1016/j.compstruct.2011.05.032 |
Cosserat E, Cosserat F. 1909. Théorie des Corps Déformables. Paris: Hermann.
|
Eringen A C. 1966. Linear theory of micropolar elasticity. Journal of Mathematics and Mechanics, 15(6): 909-923. |
Eringen A C. 1990. Theory of thermo-microstretch elastic solids. International Journal of Engineering Science, 28(12): 1291-1301. DOI:10.1016/0020-7225(90)90076-U |
Grioli G. 1960. Elasticità asimmetrica. Annali Di Matematica Pura ed Applicata, 50(1): 389-417. DOI:10.1007/BF02414525 |
Huang W X, Xu K. 2014. Characteristic length in Cosserat media modelling of granular materials (in Chinese). //Proceedings of the Conference on Computational Mechanics of Granular Materials. Lanzhou: Lanzhou University Press, 279-283.
|
Lam D CC, Yang F, Chong A C M, et al. 2003. Experiments and theory in strain gradient elasticity. Journal of the Mechanics and Physics of Solids, 51(8): 1477-1508. DOI:10.1016/S0022-5096(03)00053-X |
Madeo A, Neff P, Ghiba I D, et al. 2015. Wave propagation in relaxed micromorphic continua: modeling metamaterials with frequency band-gaps. Continuum Mechanics and Thermodynamics, 27(4): 551-570. DOI:10.1007/s00161-013-0329-2 |
Mindlin R D, Tiersten H F. 1962. Effects of couple-stresses in linear elasticity. Archive for Rational Mechanics and Analysis, 11(1): 415-448. DOI:10.1007/BF00253946 |
Mindlin R D. 1965. Second gradient of strain and surface-tension in linear elasticity. International Journal of Solids and Structures, 1(4): 417-438. DOI:10.1016/0020-7683(65)90006-5 |
Mindlin R D, Eshel N N. 1968. On first strain-gradient theories in linear elasticity. International Journal of Solids and Structures, 4(1): 109-124. DOI:10.1016/0020-7683(68)90036-X |
Polizzotto C. 2013. A second strain gradient elasticity theory with second velocity gradient inertia-Part I: Constitutive equations and quasi-static behavior. International Journal of Solids and Structures, 50(24): 3749-3765. DOI:10.1016/j.ijsolstr.2013.06.024 |
Singh B. 2002. Reflection of plane waves from free surface of a microstretch elastic solid. Journal of Earth System Science, 111(1): 29-37. DOI:10.1007/BF02702220 |
Tomar S K, Gogna M L. 1992. Reflection and refraction of a longitudinal microrotational wave at an interface between two micropolar elastic solids in welded contact. International Journal of Engineering Science, 30(11): 1637-1646. DOI:10.1016/0020-7225(92)90132-Z |
Tomar S K, Garg M. 2005. Reflection and transmission of waves from a plane interface between two microstretch solid half-spaces. International Journal of Engineering Science, 43(1-2): 139-169. DOI:10.1016/j.ijengsci.2004.08.006 |
Toupin R A. 1962. Elastic materials with couple-stresses. Archive for Rational Mechanics and Analysis, 11(1): 385-414. DOI:10.1007/BF00253945 |
Toupin R A. 1964. Theories of elasticity with couple-stress. Archive for Rational Mechanics and Analysis, 17(2): 85-112. DOI:10.1007/BF00253050 |
Wang Z Y, Li Y M, Bai W L. 2020. Numerical modelling of exciting seismic waves for a simplified bridge pier model under high-speed train passage over the viaduct. Chinese Journal of Geophysics (in Chinese), 63(12): 4473-4484. DOI:10.6038/cjg2020O0156 |
Yang F, Chong A C M, Lam D C C, et al. 2002. Couple stress based strain gradient theory for elasticity. International Journal of Solids and Structures, 39(10): 2731-2743. DOI:10.1016/S0020-7683(02)00152-X |
Zhao Y P. 2016. Modern Continuum Mechanics (in Chinese). Beijing: Science Press.
|
Zhu G, Droz C, Zine A, et al. 2020. Wave propagation analysis for a second strain gradient rod theory. Chinese Journal of Aeronautics, 33(10): 2563-2574. DOI:10.1016/j.cja.2019.10.006 |
黄文雄, 徐可. 2014. 颗粒材料Cosserat介质模拟中的特征长度. //2014颗粒材料计算力学会议论文集. 兰州: 兰州大学出版社, 279-283.
|
王之洋, 李幼铭, 白文磊. 2020. 基于高铁震源简化桥墩模型激发地震波的数值模拟. 地球物理学报, 63(12): 4473-4484. DOI:10.6038/cjg2020O0156 |
赵亚溥. 2016. 近代连续介质力学. 北京: 科学出版社.
|