三维数值流形方法研究及其在地学中的初步应用

武艳强

（中国地震局地质研究所，北京 100029）

（作者电子信箱，武艳强：chdqyw＠126．com）

伴随地震孕育、发生与震后调整过程的应力、应变问题在地震研究中占有重要地位。以GPS为代表的空间观测技术为大尺度地壳运动监测提供了高效、稳定、精确的测量结果，并为研究上述应力、应变问题提供了有效数据约束。GPS速度场和应变率场结果既可用于变形分析，又可作为数值模拟的有效约束和检验。为了获取研究区域可靠的应变率场结果，有必要开展GPS应变率场计算方法的对比研究。另一方面，由于在整个地震过程中应力、应变存在动态演化现象，因此需要研究适合于动态资料约束、大变形模拟、连续与非连续变形耦合的三维数值流形方法。基于上述分析，本文包括如下3个方面研究内容。

（1）通过模拟与实际数据分析，讨论了4种应变率场计算方法的解算精度与抗差性等问题，发展了能够较客观反映实际变形分布的最小二乘配置球面应变率解算方法，并与其他几种常用方法进行了对比分析。

首先，对大空间尺度（75°～135°E，20°～50°N）模拟数据计算结果的分析表明，采用1°×1°采样数据及其50%限定（对1°×1°数据进行50%数据量的离散化采样，并剔除2个5°×10°区域）采样数据作为输入的情况下，Delaunay三角形方法因噪声对解算结果影响太大不可取，其他3种整体应变率场解算结果具有一致性，但抗差性有所差别。通过分析理论结果与附加了不同误差的计算结果的相关系数表明，抗差性由好到坏排列如下：最小二乘配置方法、球谐函数方法、多面函数方法和Delaunay三角形方法。从输入数据的稀疏程度对不同应变率计算方法影响程度看，在数据采样率介于2°～1°网格之间时最小二乘配置方法受数据稀疏的影响最小，球谐函数和多面函数对输入数据的密集性要求较高。

其次，中等空间尺度（90°～120°E，25°～40°N）模拟数据（1°～0．5°网格）结果表明，3种整体方法在测点分布足够密的情况下均能满足实际计算需求。此时，3种方法对附加的误差敏感程度有一定差别但不显著，对比而言最小二乘配置稍强于其他2种方法。通过对不同空间采样数据的应变率计算结果的分析表明，随着输入数据越来越稀疏，多面函数和球谐函数方法计算结果与理论值的相关性减弱的幅度要快于最小二乘配置方法，表明此2种方法对数据的分布密度要求较高。

最后，通过对中国大陆1999—2004年期间应变率计算结果的分析，发现最小二乘配置方法的稳定性最高，即使仅用50%的数据作为输入，得到的应变率结果与全部数据输入基本一致，误差水平也没有明显的变化。球谐函数方法受点位分布稀疏程度的影响较大，表现为边缘效应的量值有所增大，影响范围也有所扩大。多面函数方法受测点稀疏程度的影响较大，表现为计算的不稳定性。另外，多面函数方法和球谐函数方法受点位分布的几何形状影响较大，比如中国东北地区应变率计算结果的失真。

因此总体而言，从抗差性、边缘效应、误差分布、稳定性角度来看，最小二乘配置方法最佳。究其原因在于，最小二乘配置的协方差函数是对实际观测数据统计计算得到的，能够反映数据的真实分布特征；而球谐函数的展开阶数，多面函数的光滑因子、核函数、平差结点等的选择都是通过反复试算得到的，很难做到对输入数据的最优描述。在输入数据分布密度满足的情况下，最小二乘配置方法的计算过程无需人工进行参数选择，即使不同人员进行解算也能得到一致的结果。

（2）通过单纯形积分特性研究、具体的矩阵元素表达式推导、模拟结果的理论测试等过程，从多角度对三维数值流形方法开展了研究。

首先，通过对单纯形积分公式的分析，给出了二维和三维单纯形积分的C＋＋算法实现，并结合具体实例描述了单纯形积分的求解过程。在对中空规则形状积分区域的面积（体积）和重心的计算值与理论值对比分析的基础上，对单纯形积分的精度进行了讨论。通过对比不同边长比（10－6～106）情况下计算值与理论值的差异，讨论了积分区域图形条件对积分结果的影响；通过对中空非规则形状积分区域计算结果的分析，讨论了单纯形积分的普适性特征；最后，分析了高阶多项式在中空非规则形状积分区域的积分精度。结果表明，积分项为低价或高价多项式情况下，计算值与理论值的相对误差均约为10－15～10－14，图形条件对积分结果的影响极小且不存在系统性。综合分析认为，计算结果与理论结果的差异是由计算机的计算误差造成的，单纯形积分结果具有高度的稳定性和精确性。该积分算法为本文三维数值流形方法的研究奠定了基础。

其次，在真实坐标下推导了三维数值流形方法的所有矩阵的元素表达式，具体包括单元刚度矩阵、初应力矩阵、点荷载矩阵、体荷载矩阵、惯性矩阵、速度矩阵、点位移矩阵、接触矩阵、摩擦矩阵。将所有矩阵装配起来即可实现连续、非连续三维数值流形方法模拟。由于数值流形方法不做等参变换，因此该方法的点荷载、点位移等约束可以施加到单元内任意一点。同时，在三维数值流形方法实现过程中，添加了多点位移组合约束算法，可实现相对运动约束，其求解策略包括最小二乘方法和强制约束等。

最后，通过连续变形和非连续变形具体计算实例分析，验证了固定点矩阵、点荷载矩阵、体荷载矩阵、惯性矩阵、速度矩阵、接触矩阵、摩擦矩阵等的正确性。16单元悬臂梁模拟结果与理论结果的对比分析表明程序具有高效性，同时还验证了应力逐步累积模块的有效性；多接触面剪切破坏实例模拟了弹性回跳现象，验证了开闭迭代算法的有效性，证明现有程序可以在极限平衡条件下开展静－动单次转化模拟。

通过对单纯形积分算法、矩阵元素表达式推导、关键算法描述、实例测试等研究，对三维数值流形方法的算法实现及模拟效果进行了讨论，验证了程序的有效性和模拟结果的高精度特性。

（3）汶川地震的发生与青藏高原的运动与变形密切相关，由于巴颜喀拉地块东向运动受阻，致使龙门山断裂带积累了大量的应变能。因此，开展震前区域变形特征研究有利于更加清楚地认识汶川地震孕育过程。该部分从主应力与主应变率方向差异、不同空间尺度的GPS应变率场分布特征和多步三维数值流形方法模拟等3个方面对汶川地震前地壳变形特征开展了针对性研究。

首先，利用最小二乘配置方法对中国大陆及周边World Stress Map（WSM）计划应力方向数据和GPS速度场数据进行了处理，经过数据预处理、经验协方差参数估计和拟合推估等过程，得到了中国大陆范围内应力方向和应变率场结果。数据处理过程中，对应力方向数据和GPS速度场拟合残差的无偏性和正态性进行了检验，并对应力方向和主压应变率网格结果的误差分布特征及其成因进行了分析。在充分考虑数据结果精度的基础上，讨论了中国大陆主压应力和主压应变率方向特征的异同，结果表明二者分布总体一致。在青藏地块东部存在显著差异，应力方向表现为近NE向而GPS主压应变率表现为近EW向，该特征表明青藏地块东部地区近年来东西向挤压应变积累显著；另一显著差异区域位于西域地块西部，应力方向表现为自西向东由NW向NE的转变过程，GPS主压应变率以近NS向为主。

其次，利用前文GPS应变率计算方法对比分析择优推荐的最小二乘配置球面应变率解算方法，基于1999—2007年中国大陆西部约700个GPS测点的速度场结果，分析了汶川8．0级地震前的区域应变率场分布特征。GPS主应变率场结果显示青藏高原主要表现为南北—北北东向挤压、东西—北西西向拉伸的特性。青藏中东部地区呈现由南到北、自西向东主压应变率方向逐步向东偏转的有序分布特征。青藏块体东西向GPS应变率场结果表明，块体西部（92°E以西）以东西向拉张变形为主；青藏块体东部（92°E以东、100°E以西）以东西向挤压为主，该结果反映了青藏高原物质东向流动受到了华北—华南地块阻挡的动力学特征。因此，应变能在块体边界处积累，有利于汶川地震的发生。虽然汶川8．0级地震震源区不处于剪应变率的高值区，但是处于面应变率挤压区和东西向挤压应变率高值区；发震断裂在震前存在挤压兼右旋剪切变形背景，且震源区西侧的挤压变形幅度明显大于东侧；在断裂带尺度上，应变积累较为缓慢，存在变形趋于极限现象。

最后，通过多步三维数值流形模拟方法研究了汶川地震前挤压应力特征和变形特征。结果表明，最小二乘配置GPS速度场、应变率场结果和三维数值流形方法模拟结果一致性较高，验证了三维数值流形方法模拟结果的可靠性。GPS主应变率场和三维NMM模拟结果均显示在1999—2007年间龙门山断裂带南段的应变积累速率大于中北段，中北段的挤压应变积累已趋于极限，南段依然可以继续进行应变积累。

汶川地震前地壳变形特征表明，孕震后期龙门山断裂带及其附近区域的弹性变形已接近极限，与邻近区域相比处于显著弱变形状态。龙门山断裂带南段的应变积累速率大于中北段的结果表明，龙门山断裂带中北段的闭锁程度要大于南段，更有利于地震破裂发生。因此，在研究GPS应变率场分布特征时，应该结合构造变形的背景从不同空间尺度进行分析。

总体而言，本文侧重于方法研究，其中GPS应变率场既可用于地壳变形分析又可用于检验模拟结果的正确性；三维数值流形方法的研究侧重公式推导、结果测试等，以保证算法的正确性为首要任务。文章最后给出的三维数值流形模拟实例较为简单，在地学领域进行更广泛的应用尚需开展更加深入的研究工作。