重力梯度是重力位的二阶导数，在地球物理学和大地测量学领域有重要的科学意义，在资源勘探、飞行器导航、导弹定轨等方面有不可替代的应用价值[1-3]．重力梯度张量包含6个独立分量，反映浅层地壳信息比重力更敏感，可更丰富地展现地下物质组成和地质结构情况[4]．传统重力梯度测量通过对同一方向上两个邻近相对重力仪的读数求差实现，精度较低，易受地面干扰因素影响．近年来，航空重力梯度测量技术的快速发展，使得区域全覆盖的高精度重力梯度全分量测量成为可能，其测量成本低、方式灵活，还可填补地面重力和卫星重力测量空缺频段信息，具有很强的优越性[5-10]．北美和欧洲多国已相继研发了Falcon[11]，Air-FTG[12]和HD-AGG[13]等航空重力梯度测量系统，我国也在大力推进航空重力梯度测量装置国产化进程[14]．为检验和评估航空重力梯度仪的测量精度和稳定性，须建设我国自主的航空重力梯度标定场，利用高分辨率地形数据和地面重力加密测量数据，对航空重力梯度仪测量结果进行检核和校验，以标定仪器精度．地形和地面重力数据可解算空中任意位置的航空重力梯度全分量，解算方法有直接法和等效源法两类．直接法有FFT变换法、斯托克斯积分法、径向基样条函数法及最小二乘配置法等[15-18]，由给定地面重力与航空重力梯度的关系式计算；等效源法是先建立合适的地下密度异常体模型，再据此正演计算航空重力梯度．常规等效源仅考虑点质量、多边形异常体或单层密度异常假设，不适用于地质条件复杂区域，多层等效源法丰富了地下密度随深度变化信息，相比而言更接近真实情况[19-20]．建立多层等效密度结构最直接的方法即利用三维重力反演技术[21]直接反演，此法步骤简单，可直接应用于非格网的地面重力数据，能够提供高精度的航空重力梯度全分量解算结果，密度反演结果还可进一步用于区域地质解释，但较少见到此法直接用于航空重力梯度标定的相关研究．这里围绕基于三维重力反演技术实现航空重力梯度全分量的解算展开，介绍了三维重力反演的基本原理以及解算的技术流程，通过构建理论模型证明此法的有效性和可靠性，并应用于五大连池尾山地区航空重力梯度特征研究．1 原理与方法1.1　三维重力反演三维重力反演的基本思想为：将地下空间划分为三维规则网格，每个网格视为均匀长方体，由观测面(地面)重力值反演各长方体的密度．长方体的重力公式为gz=Gρxln(y+r)+yln(x+r)-z arctan[xy/(zr)]x1x2y1y2z1z2, (1)式中：G为万有引力常数；ρ为密度；x，y，z为积分体元坐标；r=x2+y2+z2为坐标点与原点间距；x1，x2，y1，y2，z1，z2分别为x，y，z方向的积分下限和上限，此6个值可以确定长方体的空间位置和大小．三维长方体网格在一系列观测点引起的重力响应，可写为如下线性方程Lm=d，(2)即提取式(1)中的ρ构成未知向量m，剩余部分构成设计矩阵L和d为观测向量，由此建立反演观测方程．采用直接反演法，即对于观测方程式(1)构建反演目标函数ϕ=Wd(d-Lm)2+λ2Wmm2，(3)式中：Wd为观测数据权重矩阵；Wm为模型权重矩阵；λ为正则化因子．模型权重矩阵由平滑权重和深度权重两部分组成，表达式可写为Wm=αsIWr+αxBxWr+αyByWr+αzBzWr，(4)式中：αs为全局光滑因子；αx，αy，αz 为3个方向的光滑因子；Bx，By，Bz为相应的光滑权阵，这里光滑因子均取1；Wr为深度权阵，是由深度加权函数w(z)构成的对角阵，具体为w(z)=(z+z0)-β/2，(5)其中，z0为观测面高度，β为衰减因子，对于重力通常取2，重力梯度通常取3．在目标函数最小条件下求解观测方程，不考虑初始条件，正则化解的形式可写为m=(LTL+λ2WmTWm)-1LTWdTd．(6)三维重力反演先天具有解不唯一的特性，无约束条件下，正则化因子决定解的优劣性，通常采用L曲线法确定合适的正则化因子．1.2　航空重力梯度解算流程航空重力梯度的测量值包含区域地形起伏、区域浅层密度异常、地球内部质量分布这3个部分的影响．通过构建区域高分辨率地形模型、区域地下浅层三维密度结构、利用全球静态重力场模型低阶成分恢复背景场，实现航空重力梯度的理论解算．高分辨率地形模型主要用于：地面重力异常测量结果的地形改正、航空重力梯度地形效应的正演计算．参考YANG等[21]的方法，综合全球地形数据和实际高程测量数据，建立标定场区域的高分辨率剩余地形模型(DTM)．地面重力异常原始数据经地形改正和背景场改正，获得残余重力异常，利用三维重力反演技术进行无约束反演，获得标定场区域地下浅层三维密度结构．根据如下长方体重力梯度公式正演计算恢复残余航空重力梯度扰动：Γxx=Gρarctan[yz/(xr)]x1x2y1y2z1z2；Γxy=Gρlnz+rx1x2y1y2z1z2；Γxz=Gρlny+rx1x2y1y2z1z2；Γyy=Gρarctan[xz/(yr)]x1x2y1y2z1z2；Γyz=Gρlnx+rx1x2y1y2z1z2；Γzz=Gρarctan[xy/(zr)]x1x2y1y2z1z2．由此结合全球地形模型、重力场模型、实测地形数据及实测地面重力数据，可解算航空重力梯度全分量结果．将解算流程分为四步：地面重力数据获取，地形建模，三维密度结构的建立，空中重力梯度的计算．2 模型试验为检验解算方法的有效性和可靠性，设计如下模型试验：设计20 km×20 km范围的研究区域，模型由1个中心点正下方较深位置的大长方体，以及4个浅表位置的小长方体组成，模型信息在表1列出．以模型在0 km的重力值添加5%的白噪声，作为模拟地面重力观测，如图1所示．10.13245/j.hust.220914.T001表1复杂模型参数信息序号密度/(g∙cm-3)长方体尺寸/km中心位置/km11.02×2×2(-5，-5，-3)2-1.02×2×2(5，-5，-3)3-1.02×2×2(-5，5，-3)41.02×2×2(5，5，-3)50.510×10×8(0，0，-12)10.13245/j.hust.220914.F001图1模拟地面重力异常理论值与含5%噪声模拟观测图(色标单位：mGal)对模拟观测数据进行无约束三维反演，多次试验确定正则化因子为1×103，反演均方根残差值为0.029 mGal，平均残差水平为1.02%．在y=-5 km和y=5 km处(图1(b)虚线位置)作深度剖面，剖面处模型密度和反演结果对比如图2所示．由密度反演结果正算500 m高度的重力梯度全分量Γinv，与真实值Γmodel进行比较，如图3所示，差值统计信息在表2中列出．10.13245/j.hust.220914.F002图2虚线位置模型密度与反演结果深度剖面(色标单位：g/cm3)10.13245/j.hust.220914.F003图3模型重力梯度全分量真实值(左)与解算结果(右)对比图(色标单位：E)10.13245/j.hust.220914.T002表2模型解算重力梯度误差统计变量最大正误差最大负误差均方根误差Γxxinv-Γxxmodel2.63-2.960.044Γxyinv-Γxymodel1.78-1.600.025Γxzinv-Γxzmodel3.32-3.030.051Γyyinv-Γyymodel2.72-3.360.047Γyzinv-Γyzmodel4.19-4.300.054Γzzinv-Γzzmodel6.0-4.670.074E图2显示：由于重力反演的趋肤性，反演结果仅清晰呈现浅表异常体特征．但对航空重力梯度的解算影响较小，图3显示：解算结果与理论结果形态基本一致，表2显示最大均方根误差为0.074 E(Γzz分量)，平均误差水平为0.95%，低于反演的平均残差水平．重力梯度各分量解算误差由大至小的顺序为ΓzzΓyzΓxzΓxxΓyyΓxy．理论模型试验结果证明：三维重力反演技术具有良好的去噪和平滑效果，可保证航空重力梯度解算的准确性．在实际应用中，进一步分离低频背景场，可提高解算精度．3 实例分析选择五大连池尾山地区作为航空重力梯度解算测试区域．五大连池位于黑龙江中北部，地质构造复杂，五大连池火山群拥有14座新生代火山锥体，其中老黑山和火烧山在1719—1721年间喷发，距今仅有300 a，仍处于休眠状态，火山口下方可能有残余岩浆成分．尾山位于老黑山和火烧山的东北侧，张森琦等[23]在尾山地区利用重力勘探、磁法勘探、大地电磁测深、天然地震背景噪声成像四种手段开展综合地球物理勘察，地震背景噪声成像结果与大地电磁测深结果均体现，尾山下方存在正在冷却的岩浆囊[24-25]．使用的重力异常加密测量资料，基本测量比例尺为1:50 000，尾山山体附近进行1:25 000的加密测量，覆盖16 km×16 km区域．测试区域地形采用SRTM30全球数字地形模型[26]，分辨率达1'' (≈30 m)，计算重力和重力梯度地形效应考虑测区外1.5°(约166 km)的地形影响．背景场效应由EGM2008模型2~70阶球谐系数计算得到．经地形改正和背景场改正后获得的残余重力异常如图4所示，呈现尾山山体所在中心负异常，东、南、西大片正异常的特征．对残余重力异常进行三维反演，为降低边界效应影响，反演水平范围外扩至24 km×24 km区域，最大反演深度设置为地面以下15 km，正则化因子取1×103．表4列出反演残差统计信息，反演均方根残差为0.11 mGal，平均残差水平为2.58%．三维网格密度反演结果分别沿两个水平方向做成切片图，如图5所示，密度负异常体是尾山下方低密度岩浆囊的反映，周围的密度正异常体对应高密度变质岩围岩，与张森琦等[23]的结论一致．10.13245/j.hust.220914.F004图4五大连池尾山地区残余重力异常图(色标单位：mGal)10.13245/j.hust.220914.T003表3残余航空重力梯度扰动全分量解算结果统计分量最大正误差最大负误差平均值标准差δΓxx4.18-13.13-2.625.06δΓxy5.15-4.890.192.43δΓxz14.70-3.026.623.66δΓyy7.36-11.84-2.094.74δΓyz14.09-6.803.545.14δΓzz20.56-10.824.718.13E10.13245/j.hust.220914.T004表4地面重力实测数据与EGM2008模型解算航空重力梯度结果差值统计分量最大正误差最大负误差标准差Γxxin-situ-ΓxxEGM9.64-5.083.10Γxyin-situ-ΓxyEGM2.73-6.632.43Γxzin-situ-ΓxzEGM7.26-10.963.66Γyyin-situ-ΓyyEGM14.48-10.254.76Γyzin-situ-ΓyzEGM8.21-11.314.45Γzzin-situ-ΓzzEGM15.09-22.685.70E10.13245/j.hust.220914.F005图5五大连池尾山反演结果密度切片图(色标单位：g/cm3)划定中心10 km×10 km，高度500 m的区域为航空重力梯度模拟测区，计算100 m间距的残余重力梯度扰动全分量δΓ，如图6所示，统计信息见表3．图6显示：中心密度负异常体特征显著，地面重力异常的高频噪声被削弱，结果较为理想．残余航空重力梯度加上航空重力梯度地形效应和背景场效应，获得完全航空重力梯度扰动解算结果Γin-situ，与由EGM2008模型2-2160阶与地形模型计算航空重力梯度结果ΓEGM比较，如图7所示，统计信息见表4．图7显示：航空重力梯度扰动主要体现高频地形特征，但相对低频的浅表密度异常对解算结果也有一定影响，由表3可见：Γzz分量最大有20.56 E的贡献．表4显示：基于地面重力实测数据+RTM的解算结果，与EGM2008+RTM的解算结果相比，RMS最大值为5.70 E(Γzz分量)，最小值为2.43 E(Γxy分量)．主要原因为：五大连池地区浅表地质结构复杂，水平向密度变化明显，EGM2008模型最大有效水平分辨率仅10 km，无法有效进行小测区短间距的航空重力梯度解算．此外，对五大连池航空重力梯度解算误差进行估计．地面重力测量误差、地形建模误差、三维反演误差分别为0.05，1.00和0.11 mGal，根据误差传递关系计算总误差为1.21 mGal，转换为航空重力梯度估计解算精度不超过5.8 E，满足航空重力梯度标定误差要求．10.13245/j.hust.220914.F006图6五大连池尾山残余航空重力梯度扰动全分量结果图(色标单位：E)10.13245/j.hust.220914.F007图7地面重力实测数据+RTM(左)与EGM2008模型+RTM(右)解算五大连池尾山地区航空重力梯度扰动全分量结果对比图(色标单位：E)本研究介绍了将三维重力反演技术应用于航空重力梯度解算的方法和流程，并通过模型试验和实际区域应用验证了方法的可靠性和有效性．理论模型试验结果表明，基于三维反演获得的地下浅层等效密度异常结果用于航空重力梯度解算，其解算误差低于反演残差水平，同时反演的趋肤性对解算精度影响不大．选择五大连池尾山地区作为实际测试区域，利用地面重力加密测量数据、SRTM数字地形模型数据、EGM2008重力场模型数据，解算了测区上方500 m高度的航空重力梯度全分量，总误差低于5.8 E，满足航空重力梯度标定误差要求．五大连池地区地质结构复杂，采用地面实测数据与仅采用EGM2008模型解算的航空重力梯度扰动结果，最大相差5.70 E (Γzz分量)，最小相差2.43 E(Γxy分量)．在五大连池开展航空重力梯度标定工作，须采集更密集的地面重力和地面高程数据，以保证标定精度．五大连池地区正在建设我国首个民用航空重力梯度标定场，是国家重大科技基础设施建设项目“精密重力测量研究设施”(PGMF)的重要组成部分．本文基于三维重力反演解算航空重力梯度的技术路径，以及在五大连池尾山地区的测试应用，可为国产航空重力梯度测量系统的标定任务提供方法和算例参考．