航空地球物理测量数据的精度评估主要包括内符合和外符合精度评估两种方式[1-2]．内符合精度评估一般是指重复测线和交叉点不符值的评估[3-5]；外符合精度评估是在地面有实测数据的情况下，将空中测量数据向下延拓到地面，与地面实测值进行比较[6-8]，或将地面实测数据向上延拓到空中，与空中测量数据进行比较[9-11]，从而获得精度评估结果．向下延拓是一个典型的不适定问题[12-17]，随着延拓高度的增加，延拓模型的病态性将对空中测量数据中的高频干扰信号起着显著的放大作用，从而严重影响向下延拓的数据质量．因此，对于航线上测量数据的外符合精度评估一般采用向上延拓这种模式，可以比较真实地反映航空地球物理测量系统的实际性能．这里主要研究航空重力矢量测量数据的外符合精度评估方法．航空重力矢量包括东向分量、北向分量和垂直分量等三分量，其中垂直分量的外符合精度评估比较容易实现，可以采用泊松积分公式，将地面实测重力垂直分量直接向上延拓，然后将延拓值与空中数据进行比较后获得精度评估结果．重力矢量东向和北向两个水平分量，由于其并不是调和函数，因此经典的泊松积分公式无法使用．另外，将地面上的重力矢量水平分量向上延拓，代价较高，这是因为水平分量实测数据获取目前只能通过天文方法测量的垂线偏差转换得到，而这种测量手段对于天气要求特别苛刻，并且实际操作费时费力，若大规模开展，则势必须要投入大量的人力物力．因此，重力矢量水平分量的外符合精度评估，成为横亘在航空重力矢量测量系统性能评估面前的一道难题．在目前已有的地面重力数据向上延拓方法中，泊松积分法应用非常广泛[2,11,18-21]，但是由于其仅能从地面一种重力场元积分得到空中对应的重力场元(如由地面重力异常积分得到空中重力异常)，因此一般用于航空标量重力测量数据的外符合精度评估．基于斯托克斯积分方法可以由地面重力异常积分得到空中重力矢量三分量，而且相比垂线偏差测量，地面重力异常的测量相对简单、易于实现，并且采用斯托克斯积分方法向上延拓的理论也非常成熟[22-25]，故可用于航空重力矢量数据的外符合精度评估．对于经典的斯托克斯积分模型，在积分中心区，当计算点接近于边界面时，积分核函数存在奇异性，导致空中重力矢量无法正常计算，因此须要对经典斯托克斯赋值模型的奇异性进行处理．这里从物理大地测量的基本边值问题理论出发，推导了斯托克斯非奇异的严密直接赋值模型，并采用某物理大地测量试验区的实测地面重力数据，实现了对我国首套航空重力矢量测量系统的外符合精度评估．1 经典赋值模型由边值问题求解扰动位的基本公式为[18-21]T=R4π∬σΔgS(r,ψ)dσ，(1)式中：T为扰动位；Δg为地面重力异常；S(r,ψ)为斯托克斯积分核函数；r=R+H为目标点的地心向径；R和H分别为地球的平均半径和计算点的高程；ψ为目标点到流动面元dσ的球面角距．将式(1)在3个坐标轴方向求偏导，即可得到重力矢量三分量的基本计算公式为δr=∂T∂r；δφ=1r∂T∂φ；δλ=1rcosφ∂T∂λ，式中：δr为重力矢量垂直分量；δφ为重力矢量北向分量；δλ为重力矢量东向分量．由ψ与φ，λ的关系可得空中重力矢量的经典赋值模式——Stokes-Pizzetti公式[18]，具体为δr=R4π∬σΔg∂S(r,ψ)∂rdσ;δφ=-R4πr∬σΔg∂S(ρ,ψ)∂ψcosαdσ;δλ=-R4πr∬σΔg∂S(ρ,ψ)∂ψsinαdσ, (2)式中：S(r,ψ)=2Rl-3Rlr2+Rr-5R2cosψr2-3R2cosψr2lnr-Rcosψ+l2r; (3)l2=R2+r2-2Rrcosψ；cosψ=sinφ'sinφ+cosφ'cosφcos(λ'-λ)；sinψcosα=cosφsinφ'-sinφcosφ'cos(λ'-λ)；sinψsinα=cos φ'sin(λ'-λ)，其中，l为计算点至流动面元dσ的空间距离，α为计算点到流动面元dσ的方位角，{r,φ,λ}分别为计算点的地心向径、地心纬度和地心经度，{φ',λ'}为流动面元的地心纬度和地心经度．在实际计算中，通常用离散求和代替式(2)的积分形式，因此重力矢量直接赋值的实用公式为δr=R4π∑i∑jΔgij∂Sij(r,ψ)∂rΔφi'Δλj'cos φi';δφ=-R4πr∑i∑jΔgij∂Sij(r,ψ)∂ψΔφi'Δλj'cos αijcos φi';δλ=-R4πr∑i∑jΔgij∂Sij(r,ψ)∂ψΔφi'Δλj'sin αijcos φi'. (4)由式(3)可以看出：对于中心区，当计算点位于边界面时，H=0．因此，在积分中心区，当ψ→0时，则l→0，重力矢量的经典赋值公式(2)是奇异的，须要对奇异点进行处理．2 非奇异的严密直接赋值模型由物理大地测量的基本边值公式[18-21]Δg=-∂T/∂r-2T/r，可得∂(r2T)/∂r=-r2Δg(r)．(5)将式(5)在[∞,r]区间进行积分，则有r2T=-∫∞rr2Δg(r)dr．(6)利用重力异常的向上延拓公式Δg(r)=R24πr∬σr2-R2l3-1r-3Rr2cosψΔgdσ．(7)将式(7)代入式(6)，并令F(r)=-r2Δg(r)=R24π∬σ-r3-rR2l3+3Rrcosψ+1Δgdσ, (8)则有T=1r2∫∞rF(r)dr．重力矢量的垂直分量为δr=∂T∂r=1r2F(r)-2r3∫∞rF(r)dr．(9)若以计算点为原点进行极坐标ψ,α分块，则式(8)可表示为F(r)=R24π∫0π∫02π-r3-rR2l3+3Rrcosψ+1⋅Δgsinψdψdα=R24π∑i∑jΔg¯ij∫ψiψi+1∫αjαj+1-r3-rR2l3+3Rrcosψ+1sinψdψdα=∑i∑jFij(r)Δg¯ij,         (10)式中：Δg¯ij为第i，j块的平均重力异常；Fij(r)的计算公式为Fij(r)=R24π(αj+1-αj)⋅r2-R2Rl+3R2rsin2ψ-cosψψiψi+1．同理，根据式(10)，式(9)右边第二项中的积分可表示为∫∞rF(r)dr=∑i∑j∫∞rFij(r)Δg¯ijdr=∑i∑jKijΔg¯ij，式中Kij=∫∞rFij(r)dr=R24π(αj+1-αj)l2R(r+3Rcosψ)-3R2sin2ψlnr-Rcosψ+lr-rcosψψiψi+1∞r ．(11)当r→∞时，有limr→∞l=limr→∞r1-Rrcosψ+1r20+⋯=r-Rcosψ. (12)将式(12)代入式(11)，则其大括号中的值趋向于r22R-3R2sin2ψln2-3R2cos2ψψiψi+1，因此有Kij=R24π(αj+1-αj)l2R(r+3Rcosψ)-3R2sin2ψlnr-Rcosψ+lr-rcosψ-r22R+3R2sin2ψln4+3R2cos2ψψiψi+1,将其代入式(11)经整理后可得δr=∑i∑jBijΔg¯ij，(13)式中Bij=1r2Fij-2r3Kij=R24πr2(αj+1-αj)(r+R)HRl+3Rrsin2ψ⋅12+lnr-Rcosψ+l2r-lRr(r+3Rcosψ)-3Rrcos2ψ+cosψψiψi+1. (14)对于计算点，ψ=0,l=H，因此有limψ→0sin2ψlnr-Rcosψ+l2r=0．因此，式(13)在计算点周围不存在奇异问题，而且从理论上讲也不存在积分离散误差，是计算空中重力矢量垂直分量的严密公式．对于空中重力矢量的水平分量，仍用经典的赋值公式(2)计算，离散化后表示为：δφ=∑i∑jΔg¯ijAi(sinαj-sinαj+1)；(15)δλ=∑i∑jΔg¯ijAi(cosαj+1-cosαj)，(16)式中Ai=-R4πr∫ψiψi+1∂S(r,ψ)∂ψsinψdψ⋅∂S(s,ψ)∂ψ=-t2sinψ2D3+6D-8-31-tcosψ-DDsin2ψ-3ln1-tcosψ+D2,其中：t=R/r；D=l/r=1-2tcosψ+t2．式(15)和(16)的积分用辛普生数值积分法计算．从式(15)和(16)不难看出：由于sinα和cosα的存在，使得重力异常对重力矢量两个水平方向的影响具有相对坐标轴的反对称特性；因此，中心区对重力矢量水平方向的贡献为零，其影响不用计算，利用上述公式可不考虑重力矢量水平分量的奇异性．在应用式(13)、式(15)和式(16)实施赋值的过程中，与经典的赋值模型一样，也是应用不同分辨率的数据，分别对外空重力矢量进行赋值，合理的数据组合依次是5°×5°，1°×1°，20'×20'，5'×5'和1'×1'．它们对外空场的影响分别记为δr5D，δr1D，δr20F，δr5F和δr1F，对于水平分量也有类似的表示．于是地球外部空间点上的重力矢量可表示成：δr=δr5D+δr1D+δr20F+δr5F+δr1F；δφ=δφ5D+δφ1D+δφ20F+δφ5F+δφ1F；δλ=δλ5D+δλ1D+δλ20F+δλ5F+δλ1F．3 数值实验与结果分析3.1　非奇异的严密直接赋值模型的有效性验证3.1.1　离散误差比较在重力矢量的实际计算中通常采用离散求和的方法来代替全球积分，但会产生一定的计算误差，称为离散误差．为了估计离散误差的大小，在不考虑平均重力异常数据的情况下，在同样的积分区域内，分别计算了经典直接法和非奇异的严密直接法赋值中重力矢量垂直分量的积分核函数离散求和与其严格理论积分值的差异．其计算公式分别如下所示．根据式(4)可得经典赋值模型的积分核函数离散求和公式为Sr=R4π∑i∑j∂S(r,ψ)∂rΔφiΔλjcosφ¯．根据式(13)，新推导的非奇异的严密直接赋值模型的积分核函数离散求和公式为Jr=∑i∑jB¯ij．根据式(14)，严格理论值的积分核函数离散求和公式为J0=R22r2(r+R)HRl+3Rrsin2ψ12+lnr-Rcosψ+l2r-lRr(r+3Rcosψ)-3Rrcos2ψ+cosψψiψi+1.离散误差所占百分比分别为：PS=Sr-J0J0×100%；PJ=Jr-J0J0×100%．为了分析离散误差对重力矢量计算结果的影响，计算了纬度0°～60°、经度70°～150°区域不同高度处的离散误差百分比．离散求和及积分算法所利用的数据见表1所示，计算得到的离散误差百分比结果见表2．10.13245/j.hust.220911.T001表1重力异常数据分辨率数据范围重力异常总数/1031°×1°60°×80°4.8030′×30′30°×40°4.805′×5′10°×11°15.8410.13245/j.hust.220911.T002表2不同高度下的重力矢量离散误差统计表高度/km重力矢量离散误差/%PSPJ16.680.9256.390.83105.360.88506.520.761006.710.715006.900.62从表2可以看出：非奇异的严密直接赋值模型的积分核函数离散化误差明显小于经典赋值模型，并且随着计算高度的增加，其离散化误差也不断降低．这也充分凸显了新推导模型在计算空中重力矢量中的优越性．3.1.2　实际计算结果与地面实测结果比较航空重力矢量两个水平分量的精度评估是其难点问题，为了验证非奇异的严密直接赋值模型计算重力矢量水平分量的精度，利用我国某物理大地测量试验区的地面实测重力数据，计算了试验区10个实测天文/GNSS点的重力矢量水平分量，转换成垂线偏差后，与实测垂线偏差测量结果进行比较，结果见表3．10.13245/j.hust.220911.T003表3计算垂线偏差与地面实测垂线偏差比较重力矢量δλδφ最大值4.842.33最小值-6.79-6.52平均值-1.65-3.00标准差3.352.57均方根误差3.583.86mGal可以看出：计算结果与实测值差值的标准差在3.5 mGal以内，均方根误差在4.0 mGal以内(后续以均方根误差作为航空重力矢量测量数据精度水平的评价指标)，因此采用非奇异的严密直接赋值模型对航空重力矢量测量数据进行外符合精度评估的方案是完全可行的．3.2　空中重力矢量测量数据的精度评估3.2.1　测线数据外符合精度评估2017年，西安测绘研究所在某测区组织了航空重力矢量野外飞行试验，以第1架次的6条重复测线(依次编号为L1~L6)为例，重力矢量垂直分量的内符合较好，精度优于1 mGal；而由于测量结果中的姿态误差较大，导致重力矢量两个水平分量存在较大的系统误差，个别测线的最大不符值超过40 mGal，并且难以通过机场停机点的垂线偏差值对水平分量进行修正．因此，在测区设置了4个地面控制点，通过传统天文观测方法测量得到这4个点的垂线偏差信息，进而对测线上的水平分量进行系统误差修正．采用地面重力点对重力矢量水平分量数据进行校准后(垂直分量不作校准处理)，接着利用测区地面5°×5°，1°×1°，20'×20'，5'×5'和1'×1'的实测重力异常数据，由非奇异的严密直接赋值模型计算了测线上每个采样点的重力矢量三分量，并与测线实测数据进行比较，其差值结果统计见表4．10.13245/j.hust.220911.T004表4测线数据外符合精度统计测线重力矢量最大值最小值平均值标准差均方根误差L1δλ34.40-10.7611.158.7711.79δφ15.26-30.01-11.898.3411.84δr10.81-10.64-0.872.032.12L2δλ36.43-12.6911.609.2012.32δφ6.35-38.75-12.998.6212.59δr9.85-6.02-1.921.942.37L3δλ33.85-11.169.178.5010.69δφ13.16-27.52-10.467.9710.87δr7.96-10.07-2.331.732.39L4δλ35.58-13.4311.358.7611.88δφ8.44-38.29-11.928.5912.03δr10.21-7.74-1.451.882.14L5δλ33.02-12.059.838.0910.51δφ11.38-31.97-9.448.5810.70δr9.61-10.89-1.741.782.16L6δλ35.25-13.0411.399.3012.30δφ8.52-30.67-12.298.8812.42δr10.14-6.27-2.212.002.54mGal可以看出：6条重复测线的外符合精度，东向分量标准差平均值为8.8 mGal，均方根误差平均值为11.6 mGal；北向分量标准差平均值为8.5 mGal，均方根误差平均值为11.7 mGal；垂直分量标准差平均值为1.9 mGal，均方根误差平均值为2.3 mGal．3.2.2　格网数据外符合精度评估采用反距离加权方法，将所有校准后的测线数据进行格网化处理，获得空中3 550 m飞行高度的5′×5′格网数据．接着利用测区地面5°×5°，1°×1°，20'×20'，5'×5'和1'×1'的实测重力异常数据，由非奇异的严密直接赋值模型计算了每个格网点的重力矢量三分量，并与实测格网数据进行比较，其差值结果统计见表5．10.13245/j.hust.220911.T005表5格网数据外符合精度统计重力矢量δλδφδr最大值19.311.67.37最小值-33.7-20.9-10.24平均值-2.53-2.34-2.49标准差7.896.452.05均方根误差8.266.832.70mGal由表5可见：航空重力矢量垂直分量精度优于3 mGal，水平分量精度优于8.5 mGal．垂直分量测量结果与国外同时期航空重力仪的精度水平相当[26]，水平分量还存在一定差距．4 结论a．推导了计算空中重力矢量的斯托克斯非奇异的严密直接赋值模型，解决了经典斯托克斯赋值模型存在的奇异性问题；b．非奇异的严密直接赋值模型的积分核函数离散化误差明显小于经典赋值模型，凸显了新推导模型在计算空中重力矢量时的优越性；c．航空重力矢量测量系统的格网数据的外符合精度，水平分量优于8.5 mGal，垂直分量优于3 mGal．