一般来说，重力密度界面反演方法根据计算域的不同一般可以分为空间域方法和波数域方法^[18]，这两种方法都是地球物理相关工作的研究热点。空间域反演方法因其直观、方便理解的特性在该研究领域出现较早并且研究较多，一般有直接迭代法、压缩质面法和非线性反演方法等方法。然而虽然空间域算法研究成果较多，但是其大部分空间域反演方法都是离散模型，需要经过多次迭代，计算耗时较长，因此不适合用来进行大区域的界面反演。为了解决反演过程复杂、花费时间长且反演区域局限的问题，波数域密度界面正演和反演计算方法于20世纪70年代应运而生。1972年Parker首先将快速傅里叶变化引入了界面正演计算^[19]， 1974年Oldenburg在此基础上提出了密度界面波数域迭代反演法^[20]。由于该反演方法采用了 FFT（Fast Fourier Transform，快速傅里叶变换），计算速度很快，Parker-Oldenburg波数域位场迭代反演方法被广泛采用。

图1 泸定震区地形与构造简图Fig.1 Topographic and structural map of Luding M6.8 earthquake area

Parker（1973）提出了基于波数域的单界面起伏引起的重力异常表达式^[19]：

式（1）中，G是万有引力常数，Z₀指观测面，r₀是Z₀内的观测点空间坐标矢量， 为待反演二维连续界面， 为r₀在x - y平面的投影，Δρ为观测面与待反演层之间的三维剩余密度分布（或界面上下物质的密度差Δρ=ρ₁ - ρ₂）， F表示空间域(x， y)转换到波数域k的傅氏变换算子。

次年，Oldenburg（1974）通过提取式（1）右边级数表达式的第一项，得到二维界面深度的迭代反演计算表达式^[20]：

式（1）和（2）共同奠定了Parker-Oldenburg波数域快速反演算法的理论基础。

式（1）、（2）的反演实现，一般基于平面参考深度进行^[18，21]。如图2所示，将待反演界面 分解成两部分：反演界面参考深度h₀及相对参考深度的起伏 ，即：

这样可以减小界面起伏，使迭代易于收敛。

为了压制式（2）中向下延拓因子 产生的高频振荡成分，可引入正则化因子^[22]，此时式（1）、（2）可以变为：

式（4）中 为平板层h₀引起的重力变化，1 + 为正则化因子，α为正则化参数，正则化因子的最终效果相当于一个低通滤波器^[22]。式（4）、（5）构成基于平面参考深度的Parker-Oldenburg迭代反演算法。该反演算法在青藏高原周缘、长江三峡等地区均取得了较好的应用效果^[22-24]。本文采用的界面反演算法实现的原理详见文献[21]，具体实施方案参见文献[22]~[24]，这里不再赘述。

图2 三维界面反演模型图Fig.2 Three-dimensional interface inversion model diagram

布格重力异常反映了地壳内部各种密度界面和不均匀密度体产生的异常信息，通过合适的异常分离方法可以获取莫霍面起伏引起的区域重力异常并由此计算出莫霍面的分布情况^[25]。本文使用的布格重力异常数据来自全球重力场数据库（International Gravimetric Bureau）提供的全球重力场模型数据WGM2012（World Gravity Map 2012）^[26]。该重力场模型来自于地球重力场模型EGM2008和DTU10，以及从ETOPO1模型中获得的分辨率为的地形改正，本文从WGM2012重力模型中选取范围96°~108°E，22°~36°N的布格重力异常数据，并做插值处理，得到泸定震区布格重力场模型如图3所示。

图3 泸定震区布格重力异常分布图Fig.3 Distribution map of Bouguer gravity anomaly in Luding M6.8 earthquake area

重力异常由不同深度密度体信号叠加形成，其中浅部高波数信号反映局部异常，深部低波数信号对应区域异常^[27]。为提取泸定震区莫霍面相关区域场，本文采用小波变换对布格重力异常场进行多尺度分解，并通过计算功率谱斜率得到1~4阶小波细节对应的场源深度。结果表明四阶小波逼近对应的平均场源深度约40 km，可有效表征莫霍面变化^[28-29]。剔除局部高频异常（图4b），提取4阶小波逼近对应的区域重力异常，如图4a所示：区域重力异常变幅在-436~70 mgal，整体西低东高，反映了地势由西北向东南逐渐降低，莫霍面由东南向西北逐渐变深的总趋势。全区的布格重力异常以龙门山—锦屏山巨型重力梯级带为界，分为东西两大区域，其中东部的四川盆地布格重力异常值较高，等值线向西突出且分布较为稀疏，西部的松潘—甘孜地块、川滇块体北部重力场变化值较低，最低达到-436 mgal；龙门山断裂带附近重力等值线较为密集，南部的云南地区重力异常变化整体较为平缓，四川南部攀枝花地区存在南北向局部高重力异常。

重力反演计算中涉及的模型参数主要包括莫霍面平均深度∆h、壳幔密度差∆ρ、正则化参数α以及反演迭代控制的阈值δ、ε等参数^[14，21-22]。根据前人的地震测深结果可以得知，青藏高原东北缘莫霍面深度在50~62 km，四川盆地莫霍面深度约40~45 km，云南地区的整体地壳厚度在36~48 km^[6，30-31]，于是选取48.4 km为界面平均深度参与反演；壳幔密度差根据人工地震测深结果和Nafe-Drake波速-密度经验转换关系（式6）计算得出^[17，22]，震区地震测深结果显示：上覆下地壳P波速一般在6.7~7.0 km/s之间，下伏上地幔P波速在7.8~8.2 km/s之间，相应的莫霍界面的密度差一般在0.3~0.4 g/cm³。反演时，密度差以0.3~0.4范围分别按0.01步长变化，取反演结果最优的密度差为最终界面密度差，结果为0.39 g/cm³；正则化参数一般取值在0~1之间，以0.01为步长变化，在确保反演收敛且满足反演精度要求下尽可能小，最终取α=0.3；反演精度主要考虑：模型正演重力异常与观测重力异常之差的均方根小于给定阈值，莫霍面深度值与已知约束数据的最大偏差的相对均方根误差小于给定阈值，这里取δ=0.001，ε=0.01。

图4 布格重力异常信号分离图Fig.4 Signal separation diagram of Bouguer gravity anomaly

图5显示的是泸定震区莫霍面深度分布情况，其整体深度在27~76 km，由西北到东南莫霍面逐步抬升，西北部的青藏高原内部莫霍面埋深较深，在60~72 km，向SE方向莫霍面逐渐变浅，等值线向SE方向呈U型凸起；东南部莫霍面分布较浅，四川盆地莫霍面深度大约在36~42 km，川滇菱形块体自南东向北西莫霍面逐渐变深，约40~60 km，滇西地区自北向南莫霍面逐渐向上抬升，大约30~48 km，攀西地区有轻微隆起；龙门山断裂带、安宁河断裂带和金沙江断裂带上莫霍面深度变化梯度最大，与重力异常情况相呼应，反映该处莫霍面深度变化最剧烈。

将上述结果与前人利用重力反演得到的泸定震区莫霍面深度结果^[9-17]进行对比可知（见表1），本文莫霍面埋深与前人结果的基本特征相同，四川盆地深度范围与表中结果具有空间分布一致性，青藏高原内部主体深度区间与前人结果偏差2~3 km，可能与地形校正精度相关；滇西南地区偏深4~5 km，推测由于红河断裂带横向密度差异未能被常密度模型充分表征。

图5 泸定震区莫霍面深度等值线图Fig.5 Contour map of Moho depth in the Luding earthquake area

表1 泸定震区莫霍面的不同重力反演结果Table 1 Statistical results of Moho model in Luding earthquake area obtained by different gravity inversion

为了厘清差异成因，本文整合了地震测深剖面数据进行交叉验证：王椿镛等利用深地震测深剖面得到了松潘甘孜地块的莫霍面深度在52~62 km^[5]，与本文结果吻合良好；张恩会等利用深地震宽角反射/折射剖面得到了一条云南西部的二维地壳速度模型，沿测线莫霍面由西南向东北逐渐变深，深度在35~43 km^[32]，比我们的结果稍浅4~5 km；嘉世旭等布设了自四川盆地至川西北高原总长500 km的深地震探测剖面，得到四川盆地的莫霍面埋深为41~43 km^[6]，较本文结果稍深2~3 km；徐涛等利用“丽江—攀枝花—清镇“剖面的人工地震数据获得了攀西地区的莫霍面深度，莫霍面沿剖面自西向东逐步抬升，整体深度在38~53 km^[33]，在攀枝花附近有上隆特征，这与我们的反演结果相一致。重力反演对横向密度变化敏感度不足、深大断裂带处高频信号衰减过快，且Parker-Oldenburg反演算法对界面起伏较大模型具有局限性综合导致两者结果在深浅构造强烈地区深度差异较大。