精密测量院在地震形变理论模拟方法方面获进展

　　中国科学院精密测量科学与技术创新研究院地球物理与内部动力学研究团队与美国团队合作，在地震形变理论计算方法方面取得重要进展。该研究提出了求解地球变形的变系数微分方程组边值问题的近似解析方法，可计算任意球谐阶数的位错勒夫数，彻底解决了位错格林函数的收敛问题。  

　　“勒夫数”是英国科学家勒夫在研究地球的潮汐变形时定义的用来描述地球受引潮力作用而变形的两个无量纲的量，用h和k表示，分别描述对应于某球谐展开阶数的垂直位移与平衡潮潮高的比值和附加引力位与天体对地球引潮位的比值，后由日本科学家志田定义了l用来描述水平位移与平衡潮水平分量的比值，这三个比值一般统称为勒夫数，其中l也称为志田数。后来在研究地球受表面负载作用的变形也采用类似的三个无量纲的数，称为“负荷勒夫数”；研究地球受点源位错（地震）作用而变形时也有类似的三个“位错勒夫数”。勒夫数反映了地球对作用力的响应程度，其数值大小取决于地球介质特性，数值越大，说明单位力作用下的变形量越大。

　　通常，计算勒夫数的过程就是求解一个变系数的微分方程组的边值问题。除了少数几个简单的地球模型存在解析解之外，传统计算勒夫数的方法主要是四阶龙格—库塔（Runge-Kutta）数值积分方法。在该方法中，需要积分步长足够小以便每一层的各个地球模型参数可以合理地假设为常数。该团队的研究表明，龙格—库塔方法的积分步长不能超过500m，对于球谐阶数很大的情形，这个步长需要取更小。这样就需要将地球分成非常多的层数，从而计算速度较慢。特别地，当球谐阶数很大时，各变量之间的相关性将导致最后组成的求解勒夫数的矩阵奇异，因此特别高阶的勒夫数采用龙格—库塔数值积分方法是算不出来的。

　　该团队发展了一个近似的解析方法，同样将地球分层，但层厚最大可以达到100km。根据地球介质参数的分布特征，对每一层作一些近似，在地核中假设密度为常数，重力是球矢径的线性函数；在地幔和地壳中，重力是常数而密度是球矢径的反函数。通过比较原模型发现，基于这些假设所建立的新模型与原模型的差异非常小，对数值计算结果的影响完全可以忽略，因此这种假设是合理的，并且可以根据计算精度要求减小层厚使得二者之间的差异更小。这样，通过一个变量代换，能够使得原本的变系数的微分方程组变成常系数的微分方程组，从而存在解析解。该解析解由系数矩阵的特征值和特征向量表示。同时，该团队介绍了一种新的矩阵传播方法使得传播稳定高效，避免了当球谐阶数很高时，矩阵的特征值较大，从而普通的矩阵传播方法产生的溢出和误差放大导致的不稳定问题。这样，可以计算任意球谐阶数的勒夫数，从而为构造收敛的格林函数奠定了基础。

　　在球坐标系统下，格林函数是包含球函数级数的无穷项的和。对于计算点与源点距离较近的情形，由于级数的收敛非常慢，计算格林函数需要超高阶的位错勒夫数。如果直接计算超高阶位错勒夫数并进行级数求和进而获得格林函数，将使计算工作耗时而效率低下，因此需要求助加速级数收敛的特殊计算技巧，以便提高计算效率与精度。库玛（Kummer）变换是常用的加速格林函数收敛的技术手段，但是由于以往没有位错勒夫数的渐近值，在计算位错格林函数时库玛变换无法实施。特别对于震源较浅（甚至接近地表）的地震，如何准确模拟地表的形变是一个非常困难的问题。由于成功计算了超高阶的位错勒夫数，因此研究人员可以拟合出其渐近值，这样就使得库玛变换能够付诸实施，从而成功解决了格林函数的收敛问题。通过计算发现，只需要较低阶的位错勒夫数，再结合拟合出的渐近值就可以确保位错格林函数的收敛，增强了计算位错格林函数的效率。另一方面，通过分析传播矩阵的特征，也从理论上获得了位错勒夫数的精确渐近值。因此，该团队的研究成果可有效解决特殊情形下（如断层延续到地表、浅源爆炸等）格林函数收敛的计算问题。

　　与大地测量学家关注地表的变形不同，地震学家更加关注地震后地球内部的变形特征如何，如地震后断层周边的应力场如何变化。因此，该团队研究了地震引起的地球内部的变形特征。分析了点源位错引起的地球内部位移场、重力场和应变场的变化，构建了断层附近，甚至是断层面上收敛的格林函数。该研究一方面验证了所提出的解析方法可以有效应用于地球内部的变形特征研究，另一方面也为在球形地球框架下研究地震应力传输与触发问题创造了条件。

　　相关研究成果以系列文章的形式发表在国际地球物理期刊Geophysical Journal International。精密测量院研究员周江存为文章第一作者。

　　该研究得到国家自然科学基金和中科院战略性先导科技专项的联合资助。

　　图1.高阶地表位错Love数与渐近值拟合（震源深度100km）。当震源深度为100km时，尽管计算地表位错格林函数不需要特别高阶的位错勒夫数，但是通过计算高阶位错勒夫数，可以发现其相对于球谐阶数n的变化规律，由此可以很容易地拟合出当n趋于无穷大时的位错勒夫数渐近值，从而可以实施Kummer变换以加速格林函数的收敛。

　　图2. 利用位错勒夫数渐近值加速格林函数收敛（震源深度10km）。从图中可看出，不用渐近值加速时，计算收敛的格林函数需要20000阶的位错勒夫数，如果只计算到3000阶，格林函数存在明显的震荡导致的不收敛。采用渐近值加速后，计算到3000阶的位错勒夫数即可得到收敛的结果，甚至只需要1000阶位错勒夫数就能得到比较满意的结果。