目前,获取海域重力场信息的主要技术手段包括卫星重力探测、卫星测高、航空和海面船测重力测量，在陆海交界区域，还包括陆地地面重力测量^[1,2]。如何有效地整合处理由各种不同观测手段获取的具有不同频谱属性、分辨率、空间分布和误差特性的多型数据，是进一步拓展海域重力场信息保障应用必须面对的现实问题。

关于多源重力数据融合问题，国内外学者已经进行了广泛而深入的研究，并提出了许多有效的处理方法，主要有统计法和解析法两大类型，最小二乘配置(简称配置法)是统计法的典型代表，由于该方法可以联合处理不同类型的重力数据，因此在多源重力数据融合处理中得到了广泛应用^[3,4,5,6,7]。但由于配置法构建经验协方差函数必须以足够分辨率的观测数据为基础，而在实际应用中这样的条件是很难满足的，要想获得较高逼近度的协方差函数模型特别是三维空间协方差函数模型并非易事，加上协方差矩阵求逆还存在不稳定性问题，配置法融合处理效果也因此受到了较大的制约^[7,8]。在解析法方面，文献^[9,10]利用残差重力异常修正重力位模型系数的思想，提出了融合卫星、航空、地面(海面)重力数据的迭代计算方案；文献^[11]专门针对同一基准面上的同类多源数据，提出了分步平差、拟合、推估和内插相结合的解析融合处理方法。点质量方法同样属于解析法的范畴，它源于解算大地测量边值问题的Bjerhammar换置理论^[1,2]，该理论以某一虚拟球作为边界面，将离散分布于地球表面上的观测值向下解析延拓到一个完全处于地球内部的虚拟球面上，把原有的自由边界面的边值问题严格转换为固定的球面边值问题，从而可以最简单的Stokes积分公式给出问题的解。点质量方法因其计算模型结构简单、使用方便而备受关注，并在局部重力场特别是在外部重力场逼近计算中得到了较好的应用^[2]。但必须指出的是，求解虚拟点质量所涉及的向下延拓迭代过程或矩阵求逆常受到边界面和边界值粗糙化或奇异性的影响，使数值解算过程失稳，因为向下延拓总是欠适定的^[1]。最近一个时期，有多位学者将点质量方法应用于航空重力测量数据的向下延拓计算，同样发现点质量的解算过程存在明显的不稳定性^[12,13]。针对这种情况，本文将Tikhonov正则化方法引入点质量模型解算，并将其拓展应用于海域多源重力数据的融合处理。基于EGM2008位模型模拟产生航空重力和海面船测重力数据进行了融合处理仿真试验，数值计算与分析结果表明，对点质量计算模型进行正则化改造是可行、有效的。 1 点质量法计算模型

“虚拟质量”方法的理论支撑是著名的Runge-Krarup定理^[14]，即在地球外部空间定义的任一正则谐函数φ，总可用定义在地球内部任意给定的一个球外部的正则谐函数ψ来一致逼近。所谓虚拟点质量模型，就是以分布于多个Bjerhammar球面上的虚拟扰动质点系来等效地代替地球实际引力场源，只要求这些扰动质点在地球表面产生的位及其导出量以要求的精度逼近场元观测量。这种以扰动质点表征外部扰动位解的方法实际上是用质点位的线性组合来逼近地球外部扰动位，公式形式为^[2,15]:

式中，B为边界算子；∑代表地球表面；Ω为∑所界的闭点集关于三维Euler空间R³的补集；GδM代表牛顿引力常数与扰动点质量的乘积；F(Q)代表地表已知观测量。对应于重力异常、高程异常和垂线偏差观测量的边值条件方程可分别写为:

式(2)~式(5)中，(r_i,φ_i,λ_i)代表第i个观测量的球坐标；(R_j,φ_j,λ_j)代表第j个质点的球坐标；R_j=R－D_j为Bjerhammar球半径，R为地球平均半径；D_j代表第j个质点的埋藏深度；γ为地球平均重力；m为质点的个数；l_ij和ψ_ij分别代表第i个观测量与第j个质点之间的空间距离和球心角。

在实际应用中，已知观测量可以是重力异常、高程异常、垂线偏差等各类重力场参量中的一种，也可以是它们的任意组合，将它们的观测方程组成线性方程组进行联合求解，这正是点质量模型能够用于融合多源重力数据的内在特性。用L 表示已知的观测值向量，X表示未知的点质量参数向量，A表示联系已知观测值和未知点质量的系数矩阵，则有:

当已知观测量个数n等于待求点质量个数m时，方程(6)的点质量解为:

当n>m时，式(6)的最小二乘解为:

式中，P为观测值向量L的权矩阵。当L由不同类型的观测量组成时，可按定权公式确定它们的权比。假设重力异常和高程异常的测量精度分别为m_Δg和m_ζ，它们的权因子为p_Δg=σ₀²/m_Δg²，p_ζ=σ₀²/m_ζ²，取单位权中误差σ₀²=m_Δg²，则有p_Δg=1，p_ζ=m_Δg²/m_ζ²。

由式(7)或式(8)求得点质量后，即可根据实际需求，参照式(2)至式(5)计算地球表面及其外部空间任意点上的不同种类的重力场参数。设待求的重力场参量用L′表示，其对应的系数矩阵为C_P，则对应于式(7)和式(8)点质量解的重力场参数估值分别为：

2 计算模型不适定性分析与正则化改造

如前所述，点质量方法的核心思想源于重力场源的等效性原理和解析延拓理论，由于虚拟点质量解算过程涉及的向下延拓反问题在数学上总是欠适定的，因此若不采用特殊的方法求解，将得不到合理的估值。为分析造成解不稳定的原因，对系数矩阵A作如下奇异值分解^[16]：

式中，U和V分别代表A的左右奇异向量，U=[u₁,u₂,…，u_n]，V=[v₁,v₂,…，v_n]，且满足u_iu_i^Τ= I ，v_i v_i^Τ= I ；Λ 是对角矩阵，其对角线上的元素为 A的递减奇异值λ_i，即 Λ =diag[λ₁,…,λ_n]，λ₁≥λ₂≥…≥λ_n≥0。代入式(7)和式(8)，可得点质量解的统一谱分解式为：

式中，A⁺代表A^－1或(A^ΤPA)^－1A^ΤP。令重力场观测向量L=+e，其中代表观测向量真值，e为观测量的随机误差向量，则式(12)可写成：

式中，代表点质量估值的真值。如果观测方程是不适定的，则系数矩阵A的奇异值λ_i将单调递减趋于零。由式(13)知，观测噪声对真解的影响主要反映在等式右端第二项，当奇异值λ_i特别小时，即当i→∞，λ_i→0时，e 将在高频段被成倍放大，很小的观测噪声也会引起待估参数较大的偏差，其影响量甚至会掩盖掉反映真解的右端第一项。这正是解算重力场向下延拓反问题出现不稳定性的内因所在。要想获得稳定的点质量解，必须对不适定观测方程进行正则化处理，以抑制高频段观测误差对待估参数的影响。

关于正则化方法应用于不适定问题的求解，国内外已经开展了多年研究^[17]，在测量数据处理领域，特别是在航空重力测量向下延拓方面，我国学者对正则化方法的应用也做了许多有益的探讨^{[8,12,13,18,19]}。在众多正则化方法中，以Tikhonov正则化法研究最为深入，应用最为广泛。其他主要还有阻尼奇异值分解法、岭估计法、半参数模型法等。不同的正则化方法采用不同的正则化准则，本文重点研究Tikhonov正则化法的应用问题。

采用移去-恢复技术对重力场观测向量进行位模型处理，可将观测方程(6)改写为：

式中，l 代表移去位模型影响后的剩余观测值向量。 对式(14)引入Tikhonov正则化算法，其正则化准则和正则化解可分别写为^[17,18]：

式中，α为正则化参数； H 表示正则化矩阵，通常可取为单位阵。对系数矩阵 A 作奇异值分解后，可得到与式(12)相对应的正则化谱分解式为：

式(17)对应的估值均方误差为^[13]：

式中，σ₀² 代表单位权方差。由式(18)知，Tikhonov正则化估值的误差包括两部分：第一项是由测量误差引起的估值误差，随着正则化参数α的增大而减小；第二项是由正则化引起的估值偏差，随着α的增大而增大。当α=0时，第二项偏差为0，此时Tikhonov正则化估计回归传统的最小二乘估计。由此可见，相对于最小二乘估计，正则化估计是通过适当地牺牲了估计值的有偏性来换取方差的减小。选择合适的α可起到平衡式(18)右端两项误差的作用，使得它们的和为最小，这正是通过正则化算法可获得准确、稳定解的原因。因此，选择合适的正则化参数α就成为有效运用Tikhonov正则化算法的关键。根据文献^[8,16]的研究结果，认为基于广义交叉检核(general cross-validation,GCV)准则的正则化参数确定方法是可行、有效的。本文因此继续沿用了该算法。 3 数值计算与分析

为了验证正则化点质量方法在多源重力数据融合处理中的有效性，本文专门设计了基于EGM2008重力场位模型的航空重力和海面船测重力测量两种数据源的融合处理仿真试验。 3.1 试验数据与方法步骤

考虑到航空重力测量具有快速高效但精度较低，以及船载重力测量具有精度高但效率较低的特点，仿真试验将航空重力测量数据设计为全覆盖，船载重力测量数据设计为稀疏控制。使用2 160阶EGM2008重力场位模型分别计算某海域1 km、3 km、5 km高空处的重力异常 Δg₁、Δg₃、Δg₅，计算范围取为3 °×3°，格网间距为5′×5′。同时计算相同范围内格网间距为45′×45′的稀疏海面重力异常Δg₀作为控制信息。为了评估多源数据融合处理的解算效果，同时计算了相同区域内5′×5′格网的海面重力异常 Δg′ ₀。 为减弱远区效应影响，采用移去-恢复技术，取EGM2008模型的前360阶作为参考场，可分别求得相对应的残差重力异常 δΔg₀、δΔg₁、δΔg₃、δΔg₅和δΔg′ ₀。

为了模拟航空和海面重力测量数据的实际特性，对仿真的空中位模型重力异常分别加入含系统偏差的随机噪声，即观测误差e₁ (系统偏差为1 mGal(即10^－5 m/s²)，标准差为±2 mGal，均方差为±2.23 mGal)和观测误差e₂(系统偏差为2 mGal，标准差为±3 mGal，均方差为±3.61 mGal)；对仿真的海面位模型重力异常加入观测误差e₀₀(系统偏差为0，标准差为±1 mGal，均方差为±1 mGal)。具体试验方法和步骤如下。

1 )模拟产生观测误差，并将其加入空中和海面位模型重力异常 Δg，得到含误差的仿真航空和海面重力异常Δg_e；

2) 应用移去-恢复技术，从含误差的航空和海面重力异常Δg_e中移去参考场重力异常，得到航空和海面残差重力异常δΔg_e；

3) 引入Tikhonov正则化算法，并基于GCV准则确定相对应的正则化参数α。图 1展示了5 km高度对应于误差e₂情形下的GCV函数变化曲线(航高8 000 m)；

4) 基于正则化点质量法模型解算待定的点质量参数，点质量埋藏深度统一取为D_j=20 km；

5) 利用第4)步获得的点质量估值，依据式(9)和式(10)计算海面未测点上的残差重力异常 δΔg′ _0e，并在其基础上恢复参考场重力异常，得到海面未测点上的重力异常Δg′ _0e；

6) 将基于正则化点质量法的融合处理结果Δg′ _0e与原始的海面重力异常Δg′ ₀进行比较，以检验融合处理结果的精度。

图 1 观测误差为e2时的GCV函数变化曲线图Fig. 1 GCV Function for Tikhonov Regularization with Error e2

图选项

采用传统点质量法和正则化配置法，分别针对单独使用航空重力、联合使用航空和海面重力两种数据情况，以及无误差 (e₀)、误差e₁和e₀₀以及误差e₂和e₀₀三种误差情形，逐一进行了多源数据融合处理计算，并将解算结果与模拟真值进行了比较。表 1列出了 三种情形下使用不同方法和不同数据组合获得的比对统计结果。四组结果对应的四种解算模式具体定义如下：① 基于传统点质量法，单独使用航空重力数据进行向下延拓计算；② 基于正则化点质量法，并采用GCV准则选择正则化参数，单独使用航空重力数据进行向下延拓计算；③ 基于传统点质量法，联合使用航空和海面重力数据进行融合处理，其中，权因子由误差e₁和e₀₀、误差e₂和e₀₀的大小确定；④ 基于正则化点质量法，并采用GCV准则选择正则化参数，联合使用航空和海面重力数据进行融合处理，同样，权因子由误差e₁和e₀₀、误差e₂和e₀₀的大小确定。

1) 传统点质量观测方程存在比较严重的病态性是不争的事实，且病态程度随航空重力测量高度的加大而明显增大(见表 1中系数矩阵的条件数一栏)。但对于无观测误差作用下的向下延拓过程，不管是使用传统模型的模式①和③，还是使用正则化模型的模式②和④，其对应的点质量解算结果都是可靠、有效的。这说明点质量法几乎不存在由模型简化即系数矩阵构制过程引入的模型误差，正则化方法对无观测误差系统没有作用。

2) 在观测噪声作用下，传统点质量解的稳定性和精度都明显降低，且随着观测误差量值和测量高度的增大而快速降低，这与前面所作的理论分析也是相吻合的。引入正则化算法后，点质量解的稳定性和精度都得到明显改善，观测误差量值越大，测量高度越高，其相对改善幅度愈加显著。即使受到量值较大的观测误差e₂的作用，模式②在不同高度上也分别获得了标准差为2.63 mGal、3.28 mGal和3.98 mGal的比对精度，充分验证正则化算法能够有效抑制病态观测方程系数矩阵对观测噪声的放大效应。

3) 相对于单独使用航空重力数据的模式②，联合使用海面稀疏重力和航空重力数据进行点质量解算的模式④的计算结果精度在标准差和均方根差值指标上都有一定程度的改善，但改善幅度比较有限。这个结果可能跟点质量观测方程中的核函数收敛较快有关，因为其在一定程度上限制了海面稀疏重力数据的控制作用范围。因此，要想提高多源重力数据的融合处理成果质量，必须有相应密度的高精度控制测量数据作保证。 4 结 语

本文的数值计算与比对结果表明，对传统点质量模型进行正则化改造，可完全消除或显著减弱观测方程系数矩阵的病态性，有效抑制病态矩阵小奇异值引起噪声放大对点质量解造成的污染，提高解算结果的精度和稳定性。但需要指出的是，这里只对航空重力和海面船测重力数据进行了融合处理仿真试验，所得结论未必具有普适性。为了更加全面地掌握点质量计算模型的技术特性，有必要继续开展联合重力异常、高程异常和垂线偏差等多种观测量求解点质量，进而恢复地球重力场的试验和研究，以期获得对正则化点质量计算模型更全面、更准确的评价。