地球重力场模型是对地球重力场地面和空间观测数据实施数学上的解析逼近的结果，在大地测量学、地球物理学、地震学、地质学、海洋科学、国防科学、空间科学中的应用非常广泛^[1-5]。利用重力场模型和相应的球谐展开式可以计算多种扰动重力场元(重力异常、扰动重力、垂线偏差、大地水准面高等)^[6-8]。为了获得高精度、高分辨率的地球重力场模型，国内外众多机构投入了巨大的人力、物力与财力^[9-10]，目前地球重力场模型已经发展到2159完全阶次^[11]。2008年美国国家地理空间情报局(US National Geospatial-Intelligence Agency)发布了目前应用范围最广、高精度超高阶地球重力场模型EGM2008^[12-13]。2014年德国地学研究中心(Geo Forschungs Zentrum)和法国空间大地测量组(Groupe Recherches Geodesie Spatiale)联合发布了EIGEN(European Improved Gravity Model of the Earth by New Techniques)系列模型中的最新模型EIGEN-6C4^[14-16]。这两个模型是目前普遍采用的超高阶重力场模型。

卫星重力数据用来解算地球重力场模型的中低阶部分，地面与海洋重力数据用来解算地球重力场模型的高阶部分^[17-19]。目前重力场模型位系数的求解方法主要有两种：数值积分方法(numerical quadrature，NQ)又称调和分析方法，还有最小二乘方法(least squares，LS)。联合卫星重力资料确定中低阶位系数的方法也有多种，包括基于原始观测数据的处理和基于卫星重力场模型位系数的处理，其中利用最小二乘方法解算位系数的法矩阵叠加方法更具有优势，能够更加详细地顾及地面数据和卫星观测数据之间的权重关系^[20]。采用最小二乘方法构建超高阶重力场模型过程中，观测方程中的未知数(模型位系数)和观测值数量大，相应的法方程规模巨大，大型法方程的构建与直接求解难以实现^[21-23]。目前国内外利用最小二乘方法计算超高阶重力场模型，均使用超算机或者采用集成计算机网络，基于并行算法进行的，计算成本高，耗时长^[24-25]。基于最小二乘方法在求解超高阶模型问题上的优势和面临的难题，本文探索该海量计算过程中难题的解决方案，力争实现超高阶重力场模型的小存储、高效率的解算。

1 块对角最小二乘的基本原理

重力异常与位系数之间关系的基本数学模型为^[6]

式中，观测量为Δg_ij^c；待求量为位系数{C_nm^*, S_nm}；(i, j)表示格网中点重力异常数据所在的位置，代表不同的纬度和经度。该方程式是使用最小二乘法求解位系数的基础。

最小二乘方法可得位模型系数的解及相应的协方差阵^[20]

根据式(1)，矩阵A的各个元素由式(3)和式(4)表示

基于以上公式，利用全球5′×5′格网重力异常求解2159阶位系数(不包括C₀₀^*)，其中格网重力异常观测量个数为2160×4320个，未知数的个数即重力场模型位系数个数C_nm^*和S_nm共2158×2162个，因此观测矩阵A大小为(2160×4320)×(2158×2162)，法矩阵N=A^TPA大小为(2158×2162)×(2158×2162)，仅存储系数矩阵A需要开辟317 TB的存储空间，如此庞大的计算量目前是不可能依靠普通计算机完成的。

当重力数据分布于整个旋转椭球面、数据在经度方向分辨率一致、数据的权与经度无关且关于赤道对称时，权矩阵P是个对角阵，因此法矩阵N=A^TPA可表示为

为方便研究，本文近似权矩阵P为单位阵，根据球谐函数的正交特性知，由以上观测方程得到的法方程矩阵N是一稀疏矩阵，即

文献[24]给出了对位系数C_nm^*和S_nm进行重新排序的方法，根据此方法将位系数重新排列，对于N_max阶重力场模型，可利用最小二乘方法恢复且只能恢复至N_max-1阶位系数^[20]，对应的系数矩阵A可列分为4×(N_max-1)块，即

式中，A_t是C_nm^*和S_nm按一定顺序排列后，系数矩阵列分的块矩阵；未知数C₂₀^*，C₄₀^*，C₆₀^*…对应的系数列组合为第一个块A₁；C₃₀^*，C₅₀^*…系数列组合为A₂，以此类推，可以发现，块A_t的列数随着t的增加而逐渐减小，且与块A_t对应的未知参数是次m相等、阶n为非奇即偶的位系数。

以2160阶次重力场模型前2159阶位系数的最小二乘方法求解为例，N_max=2160，系数矩阵列分块个数以及生成的法矩阵块个数均为4×2159，法方程的求解可简化为对4×2159个法矩阵分别求解。这些法矩阵块中最大的块为int(N_max/2)×int(N_max/2)，即1080×1080，因此最大的系数列块矩阵A₁的大小为(2160×4320)×1080，存储A₁需要75.1 G的空间，因此对于普通性能的计算机，块对角化后的超高阶模型求解仍面临存储难题。

2 块对角化最小二乘方法的改进 2.1 法矩阵约化求解

由上述分析可知，利用分块矩阵求解法矩阵，可简化为对一个块的求解

利用式(8)求解法矩阵过程中，需要存储系数列块矩阵A_t，最大的A_t矩阵为A₁，存储该矩阵需要75.1 GB空间，较难实现，针对此问题，试分析矩阵A_t的特点，寻找简化方法。A_t可表示为式(9)、式(10)两种形式之一

由上述公式可看出，矩阵A_t具有极强的规律性，为简化矩阵计算，设

则可得

则法矩阵N_t的计算如下

综合式(11)、式(12)，可得N_t的表达式为

因为格网重力异常分布是按照全球等角格网分布的，即经度间隔是一致的，可得

则

当N_max=2160时，式(16)的近似误差在10^-10量级，式(8)的解算可简化为T_mθi^TT_mθi的求解，在计算过程中只需要存储每一个纬度的T_mθi值。

对于求解2160阶次的超高阶重力场模型而言，使用5′×5′格网重力异常，共计2160个纬度值，因此所有T_mθi值的存储量为2160×1080，即需要开辟17.8 MB的空间。综上所述，经过块对角化以及本文提出的法矩阵快速计算方案，求解一法矩阵的过程，大大节约了存储空间，提升效果见表 1。

由表 1可知，经过待求位系数的重新排序的块对角化以及法矩阵约化求解改进后，使得系数矩阵和法矩阵的必要存储空间有了大幅度的缩减，使得计算变得简单高效。

2.2 Legendre函数计算存储策略

根据重力场模型求解的数学模型及系数矩阵的表达式可知，核心的计算工作量在Legendre函数的求解。求解纬度θ_i、阶次(n, m)的缔合Legendre函数值P_nm(cos θ_i)，需要从P₀₀(cos θ_i)起算，利用已成熟的递推方法，得到所有阶次直至P_{Nmax, Nmax}(cos θ_i)的值，常见递推方式见图 1。因此，对于系数矩阵A的求解而言，需要一次求解各个纬度值所有阶次的Legendre函数值，但鉴于A矩阵过于庞大不可存取，上文对A进行了列分块。

分块后的计算工作等价于处理4×(N_max-1)个A_t矩阵的求解，然而单独计算矩阵A_t需要计算全部纬度的P_nm(cos θ_i)值，因此按照传统的递推求解方式，将所有Legedre函数值就算结果计算一次写入文件，反复读取4×(N_max-1)次，或者将Legendre函数值循环计算4×(N_max-1)次。

对于2160阶次的重力场模型求解，当一次性存取所有的Legendre函数值，需要存储变量个数为2160×(2159×2164)/2，对应的存储空间37.6 GB，这么大的数据量很难使用内存存储，如果写入文件进行读取的话会大大降低计算效率。如果每求一个A_t矩阵便遍历求解一次legendre值，那么需要循环计算次数为2160×4，这也将大大降低计算效率。

通过对系数矩阵A_t的分析，系数矩阵A_t计算过程中只使用到某一固定m的P_nm(cos θ)值，因此可以考虑采用列推法。如果列推法可行，只需存储每一列前两个值P_mm(cos θ)和P_{m+1, m}(cos θ)，当使用某一m次P_nm(cos θ)时，可直接推导该列而不用计算其他次的Legendre函数值。然而，由于Legendre函数值随着阶次的升高，其值的大小不断减小直到小于double型变量存储范围发生溢出，被存储为0，而被记为0的值在列推法中直接影响了其他较大Legendre值的求解。如图 2所示，计算纬度75°的Legendre函数值时，约600次以上的值均为0，显然是不正确的。

列推法无法应用，这里采用跨阶次递推方法进行Legendre求解，见图 1(c)，其不会发生列推法中的溢出问题。已知P_n，0(cos θ)值，利用跨阶次递推，可求得P_n，2(cos θ)的值，同理，可依次求解m=4、6、8…时P_nm(cos θ)的值，已知P_n，1(cos θ)的值，利用跨阶次递推，可求得P_n，3(cos θ)的值，同理，可依次求解m=5、7、9…时的P_nm(cos θ)值。

因此，求解A_t矩阵中需要使用m次的P_nm(cos θ)值时，可利用前一矩阵A_t-1计算过程中存储的所有纬度的P_{n, m-1}(cos θ)值，利用跨阶次递推方法，得到所有纬度的P_nm(cos θ)值并存储，原P_{n, m-1}(cos θ)值可释放。这样，Legendre的计算只需要进行一次且存储问题也得到了解决。采用该方法需要存储m=0和m=1的所有纬度和“阶数n”的P_nm(cos θ)值，共存储2160×2159×2个数值，所需存储量71.2 MB，该数据量在普通计算机上存储已经绰绰有余。

2.3 Legendre函数对称性

通过对系数矩阵A_t进行分析可知，当n-m为偶数时，A_t矩阵上半部分与下半部分是对称的，当n-m为奇数时，A_t矩阵上半部分与下半部分数值相等，但符号相反。究其原因，是由于全球格网数据分布是关于赤道对称的，满足θ_i=π-θ_N-1-i，且缔合Legendre函数满足以下特性

因此在计算矩阵A_t及法矩阵时，只需要计算上半部分即可，最终求得的法矩阵N乘以2即为最终的法矩阵，该方法能够节省一半的存储空间及一半的计算时间。

3 试验及分析 3.1 精度改进

采用上述改进的计算方法和存储方式，利用EGM2008模型2159阶位系数模拟计算全球5′×5′格网重力异常点值数据，实现2159阶次超高阶重力场模型的最小二乘方法实现，并与EGM2008模型前2159阶位系数进行比较，作为模型恢复方法精度的衡量标准，并与轮胎调和分析方法恢复重力场模型的精度进行对比分析^[26]，结果见图 3。图 3分别统计了两种方法恢复的模型系数与EGM2008系数差的绝对值的自然对数。

使用EGM2008作为标准模型，计算两种方法恢复的2159阶模型系数的误差阶RMS并进行比较，比较结果见图 4，其计算公式为^[3]

由图 4可看出，块对角最小二乘方法恢复模型系数的误差阶RMS明显小于使用数值积分方法恢复的模型系数，最小二乘方法恢复的模型整体精度比数值积分方法高出至少5个数量级。

3.2 效率改进

为验证本文所提改进方法在构建全球重力场模型的优势，使用3.2 GHz主频、2 GB内存PC计算机，基于vs2008试验平台，进行模拟计算。

试验统计了分别使用3°×3°、1°×1°、30′×30′及5′×5′的格网重力异常数据，利用传统块对角最小二乘法恢复59阶、179阶、359阶和2159阶次模型的计算时间，并利用改进后的最小二乘方法进行同样的计算，统计结果见表 2。

如表 2所示，采用块对角最小二乘方法计算模型位系数，效率相对低下，计算359阶位系数需要近2 h，对于2159阶重力场模型，由于其系数矩阵A_t存储量庞大，无法实现系数的求解。相比传统块对角最小二乘法，采用改进后的方法计算效率有了明显提升，与恢复359阶模型使用2 h相比，改进方法只需要20 s。更重要的是，改进后的方法实现了2159超高阶重力场模型的求解，耗时不足2 h。

综上，改进后的最小二乘方法不仅实现了普通计算机上超高阶模型的位系数求解，并且效率得到了大的提升。

3.3 数据的精度估计

为了说明相比数值积分方法，最小二乘方法在精度评估方面的优势，本文对EGM2008模型前359阶位系数模拟的30′×30′格网重力异常加入白噪声，利用最小二乘方法恢复模型位系数过程中，记录V^TPV的估计值，利用式(20)可获取观测值的单位权中误差

式中，n为观测量个数；t为未知数个数。

通过对上述30′×30′重力异常数据分别加入5×10^-5 m/s²和10×10^-5 m/s²白噪声后，进行模型位系数的恢复，最终得到的统计结果见表 3。

根据表 3结果可知，最小二乘方法能够在一定程度上反映出数据与模型的吻合程度，即一定程度上评估数据的精度情况，这是最小二乘相对于数值积分方法不可比拟的优势。

4 结论

数值积分方法实质上就是法矩阵为对角阵的最小二乘方法，即认为模型系数之间均是不相关的，协方差均为0，而块对角最小二乘则估计了同次不同阶系数之间的相关性，其更符合实际情况，也能够恢复更高精度的重力场模型。数值积分方法不能对重力异常数据进行精度评估，块对角最小二乘方法可以。相比数值积分方法，最小二乘方法无法求解N_max阶系数，例如5′×5′分辨率的格网重力异常只能求解2159阶次的系数，第2160阶的系数必需使用数值积分方法才能求得。另外，数值积分方法使用的重力异常数据通常为点值，即包含全频段信号，所反演的系数才是实际模型位系数，而最小二乘方法所使用的格网重力异常数据应只包含所求阶次的能量，否则其信号在位系数之间会发生泄漏。

本文基于Pavlis提出的块对角最小二乘方法，通过对法矩阵求解方程进行约化、对Legendre函数的计算和存储方式进行设计，结合缔合Legendre函数关于赤道的对称性，解决了大型矩阵存储及计算效率低下的难题，实现了超高阶重力场模型最小二乘方法的小存储、高效率的解算。结果显示，精度相比数值积分方法提高近5个量级，计算效率提高近300倍，并精确估计了重力异常数据的精度。

本文最小二乘改进方法的局限性：采用上述高效的超高阶模型位系数最小二乘方法，在式(8)中，有权矩阵P为单位权阵这一假设。对其进行拓展，当位于同一经度的格网重力异常精度相同，或位于同一纬度的格网重力异常精度相同，这两种情况下也可以使用上述改进的快速方法。当出现每一片区域都有自己不同的精度这一情况下，上述方法是不适用的，这也是该改进方法最大的局限性。