地形和岩石圈内部密度横向分布不均匀是水平构造应力的重要来源。向文等^[1]和郭飞霄等^[2]研究了重力场和构造应力场间的转换关系，Flesch等^[3-4]构建了薄板模型近似下利用重力势能水平梯度计算偏应力的有限元方法，Ghosh等^[5-6]指出重力势能水平梯度应该与三维定义的偏应力平衡，势能受岩石圈结构和产生动力地形的径向作用力影响，可通过均衡调整消除。本文采用偏应力作为构造应力，利用Crust 1.0地壳模型^[7]、SIO V24.1重力场模型^[8]和SIO V18.1地形模型^[9]等数据计算华北地区(110~123°E，33~43°N)重力势能，采用有限元方法求解重力势能差产生的水平构造应力，并将水平构造应力最大主方向与世界应力图^[10]中震源机制解结果进行比较。

1 计算原理

薄板模型近似下重力势能和偏应力之间的平衡方程为^[6]：

(1)

式中，，τ_ij为偏应力的ij分量，h为地表地形，L为参考深度，本文取L=100 km^[4-6]。GPE为重力势能^[11]，定义为：

(2)

式中，ρ(z′)为深度z′处的密度，g为重力加速度，z和z′为深度积分变量。

偏应力满足τ_xx+τ_yy+τ_zz=0。式(1)中还有3个待求量τ_xx、τ_yy、τ_xy，偏应力解不唯一。文献[4]以偏应力第二不变量最小作为约束条件，根据Mises屈服准则，某点偏应力第二不变量数值越小，则该点变形状态越有可能处于弹性状态^[12]，在短时间内可将岩石圈形变当作弹性变形^[13]。构建泛函：

(3)

式中，S为研究区域，τ_γγ=(τ_xx+τ_yy)，τ_αβ为偏应力，λ_α为拉格朗日乘子，δ_αβ为Kronecker符号，α和β分别可取x和y。将泛函I对偏应力τ_αβ求变分，可得到区域S边界上拉格朗日乘子为零，区域内拉格朗日乘子和偏应力之间的关系为：

(4)

将式(4)代入式(1)得满足拉格朗日乘子的平衡方程：

(5)

式中，∂λ_γ/∂x_γ=∂λ_x/∂x+∂λ_y/∂y。

将下面泛函J^[4]对拉格朗日乘子求变分得到式(5)：

(6)

因此，泛函J对拉格朗日乘子取最小即能建立有限元求解方程^[4-6]，解算研究区域偏应力的唯一解。

2 数据处理 2.1 数据来源

地表地形数据采用SIO 2014-11-28发布的V18.1地形模型，其陆地地形采用ETOPO1数据，海底地形为联合船测海深和卫星测高重力异常数据构建。图 1给出了华北地区的地形图，其中活动断块边界取自文献[14]。

采用SIO V24.1重力场模型计算重力加速度，该模型陆地重力异常来自EGM2008模型，海域重力异常来自卫星测高。图 2为华北地区重力异常分布。

岩石圈密度数据采用Crust 1.0地壳模型，其空间分辨率为1°，地壳分为8层：水、冰、上沉积、中沉积、下沉积、上结晶地壳、中结晶地壳、下结晶地壳，给出了每一层的边界深度和密度(含莫霍面以下)。

2.2 数据处理

WGS-84椭球面上的正常重力为：

(7)

式中，γ₀为正常重力(单位：mGal)，B为大地纬度。

将γ₀进行空间改正归算为似地形面上的正常重力γ_Q：

(8)

式中，H^*为正常高，其数值采用V18.1地形模型结果。γ_Q加上重力异常得地表重力加速度，华北地区地表重力加速度分布见图 3所示。

根据式(2)计算重力势能，从地表积分到参考深度100 km，Crust 1.0地壳模型的水层和冰层不会产生显著的构造应力^[6]，在计算时剔除，且模型的空间分辨率1°较低，采用样条插值得到内部点数值，积分时每一层密度取常数。

图 4给出了直接采用V18.1地形模型和Crust 1.0地壳模型计算得到的重力势能，该结果含来自参考深度100 km以下地幔密度浮力产生的径向作用力的影响，为消除该影响在计算时进行均衡调整改正^[6]。文献[6]指出，Airy模式均衡调整改正可通过消除动力地形部分影响来实施，即将地壳积分区间长度减去动力地形^[15]。图 5为动力地形分布，可见华北平原断块、河淮平原断块、胶辽断块和苏沪-南黄海断块的动力地形接近零，负区域主要集中在太行山断块西北和华北平原断块内部济南和淄博附近。经动力地形改正得到的华北地区重力势能分布如图 6所示，说明势能分布和地形密切相关，即地形高则势能大，地形低则势能小。

图 7为经动力地形改正与未经改正得到的重力势能数值之差，表明在华北地区西北部前者大于后者，其余地区两者接近。利用势能计算水平构造应力时，需将势能减去参考值得到重力势能差^[16]。本文参考势能采用经过均衡调整得到的华北地区平均值。

重力势能差产生的水平构造应力采用有限元方法进行计算，为消除边界效应，将华北地区边界往外扩大3°。采用四边形格网对整个区域进行划分，格网一共有1 406个节点、1 332个单元，图 8为华北地区水平构造应力分布。

3 结果分析 3.1 重力势能分布

华北地区重力势能变化范围为1 449.412~1 511.755 MPa，平均为1 485.409 MPa。阴山-燕山断块中北部凌源、赤峰、承德和沽源附近为高势能区，最大为1 511.755 MPa；华北平原断块东北部葫芦岛、秦皇岛一带与河淮平原断块西北部和华北平原断块西南部交界区域的郑州和卫辉周围为低势能区，最小为1 449.412 MPa。

在重力势能中消除动力地形影响时，负动力地形会增大重力势能，正动力地形则会减小。图 5显示华北地区西北部动力地形为负，图 7表明经动力地形改正得到的重力势能比未经改正结果要大，而在动力地形接近零的华北地区东南部，两者接近。山西断陷带、鄂尔多斯断块东部、河套断陷带东部和阴山-燕山断块的负动力地形，使得这些地区的势能比未经动力地形改正的结果要大。

3.2 重力势能差水平构造应力分布

华北地区重力势能差水平构造应力空间分布不均匀，压缩应力和引张应力交替出现，应力大小变化范围为0.112~11.105 MPa。水平构造应力分布和重力势能分布密切相关，势能高区应力状态受引张控制，低区则为压缩。文献[6]和^[13]给出了华北地区250 km尺度的重力势能差水平构造应力分布，数值变化在2~12 MPa之间，水平构造应力分布样式与本文结果一致，但本文空间分辨率更高。华北平原断块西部衡水安阳附近最大主压方向为NWW-EES，胶辽断块西部平度周围为NWW-EES，阴山-燕山断块中北部凌源至赤峰一带重力势能差应力状态为N-S向引张。

在一定的条件下，由多个地震震源机制解的压缩轴、零轴和张力轴可以反映该区域构造应力场的最大、中等和最小主压应力方向^[17]。图 9为世界应力图震源机制解和本文计算得到的最大主方向分布，可知震源机制解最大主方向为NEE-SWW，对比发现两者存在偏差，特别是在北京、沧州、临沂、沂水、汝阳、沁水和大同等地方，震源机制解几乎与本文结果垂直。而本文计算得到的水平构造应力仅考虑了重力作用，且是岩石圈的平均响应，最大主方向与震源机制解结果不一致说明重力不是华北地区构造应力的主要来源。

4 结语

1) 华北地区水平构造应力空间分布与重力势能密切相关，高势能区构造应力以引张为主，低势能区则由压缩控制。阴山-燕山断块中北部为全区最大的高势能区，以N-S方向引张应力分布明显。华北平原东北部、西南部和河淮平原西北部为低势能区，其应力状态表现为压缩。

2) 华北地区偏应力与震源机制解的最大主方向不一致，特别在北京、沧州、临沂、沂水、汝阳、沁水和大同等地方，几乎垂直，说明重力不足以支撑华北地区的构造变形，且抵抗板块运动产生的驱动力。