CN110705022A - 一种稀疏球面径向基函数局部重力场建模方法 - Google Patents

Abstract

本发明公开了一种稀疏球面径向基函数局部重力场建模方法，步骤如下：S1：搜集多源重力观测数据；S2：将多源重力观测数据中的长和/或短波部分移去；S3：构建残余位的球面径向基函数模型；S4：根据重力异常与重力扰动位之间的泛函关系和残余位的球面径向基函数模型，确定出各观测量的观测方程；S5：对观测方程中的待估参数进行求解；S6：对正则化超参数进行优选，并将优选的正则化超参数对应的待估参数作为求解后的待估参数；S7：将求解后的待估参数代入残余位的球面径向基函数模型中，恢复长和/或短波部分，得到扰动重力位模型，并导出其他重力场泛函模型。本发明得到的稀疏模型不仅更简单，且精度更高，同时还可以实现自动的核函数选取。

Description

一种稀疏球面径向基函数局部重力场建模方法

技术领域

本发明涉及局部重力场建模技术领域，尤其涉及一种稀疏球面径向基函数局部重力场建模方法。

背景技术

相比传统的采用球谐函数模型进行重力场建模的方法，球面径向基核函数法是一种具有局部支撑特性的灵活方法，特别适用于采用局部重力观测数据进行局部重力场建模问题，因此该方法在进入新千年后得到了国际范围内的广泛关注。

球面径向基核函数法的灵活性来源于其存在较多的超参数，超参数调节的合理与否对最终的建模效果具有重要影响。这些超参数中最为重要的是核函数的位置和数量，即核函数的放置问题，当然二者是相互联系的，同等(层内)放置密度情况下，层数翻倍(一层对应着一个埋藏深度)，对应着待估计参数个数的翻倍。然而这些超参数的调节是一个相当棘手的问题，至今没有很好的方法和准则处理这一问题。比如核函数的个数不能太多也不能太少，太多则会引起问题的病态且所建立的模型过于复杂从而不易应用，太少则无法精确地表示建模对象。虽然常用的L2范数正则化方法可以较有效的解决病态问题，但无法解决模型过于复杂的问题。虽然有些文献提出了一些基于数据驱动的超参数自适应选择方法，但这些方法的自适应性其实很低，仍然需要进一步改进。

某一组超参数对对应着一个新的模型，因此超参数的选择问题实际上是一个模型选择问题。然而传统的信息量准则等模型选择理论具有诸多局限性，比如这些模型选择理论需要逐个地评估模型性能，而备选模型的个数是很多的，这就涉及非常大的计算量，一般是一个NP难问题。

发明内容

发明目的：针对球面径向基函数局部重力场建模中径向基核函数个数及位置选择的问题，本发明提出一种稀疏球面径向基函数局部重力场建模方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：

一种稀疏球面径向基函数局部重力场建模方法，所述建模方法具体包括如下步骤：

S1：搜集建模区域内的多源重力观测数据；

S2：将所述多源重力观测数据中的长波部分和/或短波部分移去，获取得到残余重力观测数据；

S3：通过所述残余重力观测数据，构建残余位的球面径向基函数模型，用以表示所述残余重力观测数据中不同位置处的残余位的大小，所述残余位的球面径向基函数模型具体为：

其中：T_res(x)为x位置处的残余位的大小，N为核函数的个数，β_j为第j个核函数的系数，x为残余位的位置，z_j为第j个核函数的位置；

S4：根据重力异常与重力扰动位之间的泛函关系和所述残余位的球面径向基函数模型，确定出各观测量的观测方程；

S5：通过稀疏正则化对所述观测方程中的参数向量进行求解；

S6：对所述待估参数求解过程中产生的正则化超参数进行优选，确定出优选的正则化超参数，并将所述优选的正则化超参数所对应的待估参数作为求解后的待估参数；

S7：将所述求解后的待估参数代入残余位的球面径向基函数模型中，并恢复所述残余重力观测数据中的长波部分和/或短波部分，得到扰动重力位模型，并导出其他重力场泛函模型。

进一步地讲，在所述步骤S2中，根据参考重力场模型将所述多源重力观测数据中的长波部分移去，根据地形数据将所述多源重力观测数据中的短波部分移去。

进一步地讲，在所述步骤S3中，所述第j个核函数的位置z_j的大小通过多层方法进行选定，具体为：

任意一个所述观测量的位置均对应有两个核函数的位置，所述两个核函数的位置对应有两个不同的深度，且所述两个核函数的位置中的纬度和经度与观测量的位置纬度和经度相同，所述选定的两个核函数位置，具体为：

其中：z₁为第一个选定的核函数的位置，z₂为第二个选定的核函数的位置，R_j为第j个核函数的位置所对应的地心向径，为第j个残余位的位置所对应的纬度，λ_j为第j个残余位的位置所对应的经度，R₁为第一个选定的核函数对应的地心向径，R₂为第二个选定的核函数对应的地心向径。

进一步地讲，在所述步骤S4中，确定出各观测量的观测方程，具体如下：

S4.1：根据所述重力异常与重力扰动位之间的泛函关系，由所述残余位的球面径向基函数模型获取残余重力观测数据中重力异常观测的表达式，具体为：

其中：Δg(x_i)为残余重力观测数据中第i个重力异常观测，N为核函数的个数，B_ij为观测矩阵B中的第i行、第j列元素，β_j为第j个核函数的系数；

S4.2：通过所述残余重力观测数据中重力异常观测的表达式，确定出各观测量的观测方程，具体为：

y＝Bβ+e

其中：y为Δg(x_i)堆栈而成的观测向量，B为观测矩阵，β为β_j堆栈而成的参数向量，e为观测误差，Δg(x_i)为残余重力观测数据中第i个重力异常观测，β_j为第j个核函数的系数。

进一步地讲，在所述步骤S4.1中，所述观测矩阵B中的第i行、第j列元素B_ij，具体为：

其中：

B_ij为观测矩阵B中的第i行、第j列元素，r_i为第i个观测位置x_i所对应的地心向径，R_j为第j个核函数的位置z_j所对应的地心向径，z_j为第j个核函数的位置，x_i为第i个重力异常观测所在的位置。

进一步地讲，在所述步骤S5中，通过稀疏正则化对所述观测方程中的参数向量进行求解，具体如下：

S5.1：通过各观测量的观测方程，确定协方差矩阵的大小，具体为：

Q＝cov[e]

其中：Q为协方差矩阵，e为观测向量的观测误差；

S5.2：根据协方差矩阵的大小，对所述观测方程中的参数向量进行基追踪原子分解，具体为：

其中：

为参数向量的估计值，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，B为观测矩阵，β为β_j堆栈而成的参数向量，β_j为第j个核函数的系数，Q为协方差矩阵，μ为正则化超参数；

S5.3：通过FISTA方法对所述参数向量的估计值进行求解，确定出求解后的参数向量。

进一步地讲，在所述步骤S5.3中，通过FISTA方法对所述参数向量的估计值进行求解，具体如下：

S5.3.1：将所述核函数的系数进行初始化，具体为：

其中：θ₁和s₁均为第一次迭代运算引入的辅助变量，β₀为迭代前核函数系数的初步估计值，B为观测矩阵，Q为协方差矩阵，μ为正则化超参数，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，I为单位矩阵；

S5.3.2：根据所述迭代前核函数系数的初步估计值，进行迭代计算，确定出第k次迭代时核函数系数的估计值，具体为：

其中：β_k为第k次迭代时核函数系数的估计值，

为软阈值算子，θ_k+1和s_k+1均为第k+1次迭代运算引入的辅助变量，β_k-1为第k-1次迭代时核函数系数的估计值，s_k为第k次迭代运算引入的辅助变量；

S5.3.3：根据所述第k次迭代时核函数系数的估计值、第k-1次迭代时核函数系数的估计值以及预设阈值的大小，当三者关系满足如下公式时，结束迭代计算，所述公式具体为：

||β_k-β_k-1||_∞≤ε

其中：β_k为第k次迭代时核函数系数的估计值，β_k-1为第k-1次迭代时核函数系数的估计值，ε为预设阈值；

S5.3.4：将所述第k次迭代时核函数系数的估计值作为参数向量的估计值的大小。

进一步地讲，在所述步骤S5.3.2中，所述软阈值算子的求解过程，具体如下：

第一步：给定一个向量a，对所述向量进行软阈值操作，获取得到一个同维数的向量，所述向量的第j个元素，具体为：

其中：

a表示软阈值算子中的向量，[a]_j为a的第j个分量，(|[a]_j|-μt)₊为合页函数，sgn([a]_j)为符号函数，θ_k为第k次迭代运算引入的辅助变量，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测；

第二步：对所述合页函数进行求解，具体为：

其中：

(|[a]_j|-μt)₊为合页函数，[a]_j为a的第j个分量，μ为正则化超参数，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测；

第三步：对所述符号函数进行求解，具体为：

其中：sgn([a]_j)为符号函数，[a]_j为a的第j个分量。

进一步地讲，在所述步骤S6中，对所述待估参数求解过程中产生的正则化超参数进行优选，具体如下：

S6.1：从预先设置的M个正则化超参数中选取一个正则化超参数；

S6.2：将所有所述多源重力观测数据通过无放回随机采样方法均分为10组，并从所述10组数据中，依次选出一组数据作为检核组的数据，并根据剩余9组数据重复步骤S1-步骤S5，建立所述残余位的球面径向基函数模型；

S6.3：通过所述残余位的球面径向基函数模型对检核组中的数据进行预测，获取所述检核组的预测值，根据所述检核组的实际观测值，将所述实际观测值与预测值相减得到预测误差，同时计算所述预测误差的平方和，并将所述10组数据分别作为检核组时所对应的预测误差的平方和相加，获取得到所述选取的正则化超参数所对应的交叉检核误差；

S6.4：除所述选取的正则化超参数外，从所述预先设置的M个正则化超参数中重新选取一个正则化超参数，并重复步骤S6.2-步骤S6.3，直至所述M个正则化超参数全部被选取；

S6.6：从所述M个正则化超参数各自对应的交叉检核误差中选出最小值，所述最小交叉检核误差对应的正则化超参数即为优选的正则化超参数。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

(1)本发明的稀疏球面径向基函数局部重力场建模方法通过选取一组冗余的径向基核函数，并采用稀疏正则化方法对这些核函数的系数进行求解，从而不仅得到的稀疏模型更简单，且精度更高，同时还可以实现自动的核函数选取，提高球面径向基核函数局部重力场建模技术的自动化程度；

(2)本发明的稀疏球面径向基函数局部重力场建模方法适用于多源、不规则分布的重力场观测数据，且无需对重力场观测数据进行向下延拓、格网化等预处理，同时该建模方法还能够选择足够多的核函数，且核函数的个数甚至可以多于观测量的个数，从而具有更好的表示能力。

附图说明

图1是本发明的数据区域和建模区域的分布示意图；

图2是本发明的稀疏球面径向基函数局部重力场建模方法的流程示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。

实施例1

参考图2，本实施例提供了一种稀疏球面径向基函数局部重力场建模方法，具体包括如下步骤：

步骤S1：搜集建模区域内的多源重力观测数据。在本实施例中，该多源重力观测数据包括重力异常数据、扰动重力数据、垂线偏差数据、重力梯度数据、高程差数据。同样的，在建模区域以外或附近区域中有对建模有作用的数据，也需要进行搜集。参考图1，可以看出数据区域的范围可以大于建模区域。

步骤S2：将多源重力观测数据中的长波部分和/或短波部分进行删除，获取得到残余重力观测数据。具体地讲，根据参考重力场模型删除多源重力观测数据中的长波部分。当多源重力观测数据中的短波部分也需要删除时，根据地形数据删除短波部分。其中，多源重力观测数据中的短波部分是否删除取决于使用者的自身选择。

步骤S3：通过残余重力观测数据，构建残余位的球面径向基函数模型，用以表示所述残余重力观测数据中不同位置处的残余位的大小。具体地讲，残余位的球面径向基函数模型，具体为：

其中：T_res(x)为x位置处的残余位的大小，N为核函数的个数，β_j为第j个核函数的系数，x为残余位的位置，z_j为第j个核函数的位置。

在本实施例中，第j个核函数的位置z_j的大小通过多层方法进行选定，具体为：

任意一个观测量位置x均对应有两个核函数的位置，该对应的两个核函数的位置对应有两个不同的深度，即选定的两个核函数的位置与平均地球表面的距离不同。且两个核函数位置中的纬度和经度为相应观测量的位置x的纬度和经度相同。从而第j个残余位的位置x的大小为：

时，两个核函数的位置，具体为：

其中：z₁为第一个选定的核函数的位置，z₂为第二个选定的核函数的位置，R_j为第j个核函数的位置所对应的地心向径，

为第j个残余位的位置所对应的纬度，λ_j为第j个残余位的位置所对应的经度，R₁为第一个选定的核函数对应的地心向径，R₂为第二个选定的核函数对应的地心向径，r_j为第j个残余位的位置所对应的地心向径。

在本实施例中，具体地讲，第一个选定的核函数对应的深度R₁、第二个选定的核函数对应的深度R₂，具体为：

其中：R₁为第一个选定的核函数对应的地心向径，R₂为第二个选定的核函数对应的地心向径，R为地球平均半径。

步骤S4：根据重力异常与重力扰动位之间的泛函关系和残余位的球面径向基函数模型，确定出各观测量的观测方程，具体如下：

步骤S4.1：根据重力异常与重力扰动位之间的泛函关系，由步骤S3中的残余位的球面径向基函数模型，可以获取得到残余重力观测数据中重力异常观测的表达式，具体为：

其中：Δg(x_i)为残余重力观测数据中第i个重力异常观测，N为核函数的个数，B_ij为观测矩阵B中的第i行、第j列元素，β_j为第j个核函数的系数。

在本实施例中，观测矩阵B中的第i行、第j列元素B_ij，具体为：

其中：

B_ij为观测矩阵B中的第i行、第j列元素，r_i为第i个观测位置x_i所对应的地心向径，R_j为第j个核函数的位置z_j所对应的地心向径，z_j为第j个核函数的位置，x_i为第i个重力异常观测所在的位置。

步骤S4.2：通过残余重力观测数据中重力异常观测的表达式，确定出各观测量的观测方程，具体为：

y＝Bβ+e

其中：y为Δg(x_i)堆栈而成的观测向量，B为观测矩阵，β为β_j堆栈而成的参数向量，e为观测向量的观测误差，Δg(x_i)为残余重力观测数据中第i个重力异常观测，β_j为第j个核函数的系数。

步骤S5：通过稀疏正则化对观测方程中的参数向量进行求解，具体如下：

步骤S5.1：通过各观测量的观测方程，确定协方差矩阵的大小，具体为：

Q＝cov[e]

其中：Q为协方差矩阵，e为观测向量的观测误差。

步骤S5.2：根据协方差矩阵的大小，通过Lasso建模方法，对步骤S4.2中观测方程中的参数向量β进行基追踪原子分解，具体为：

其中：

为参数向量的估计值，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，B为观测矩阵，β为β_j堆栈而成的参数向量，β_j为第j个核函数的系数，Q为协方差矩阵，μ为正则化超参数。

步骤S5.3：通过FISTA方法对参数向量的估计值进行求解，确定出求解后的参数向量。其中通过FISTA方法对参数向量的估计值进行求解，具体如下：

步骤S5.3.1：将核函数的系数进行初始化，具体为：

其中：θ₁和s₁均为第一次迭代运算引入的辅助变量，β₀为迭代前核函数系数的初步估计值，B为观测矩阵，Q为协方差矩阵，μ为正则化超参数，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，I为单位矩阵。

步骤S5.3.2：根据迭代前核函数系数的初步估计值β₀，进行迭代计算，确定出第k次迭代时核函数系数的估计值，具体为：

其中：β_k为第k次迭代时核函数系数的估计值，

为软阈值算子，θ_k+1和s_k+1均为第k+1次迭代运算引入的辅助变量，β_k-1为第k-1次迭代时核函数系数的估计值，s_k为第k次迭代运算引入的辅助变量。

在本实施例中，软阈值算子

的求解过程，具体如下：

第一步：给定一个向量a，对该向量进行软阈值操作，获取得到一个同维数的向量，该向量的第j个元素，具体为：

其中：

a表示软阈值算子中的向量，[a]_j为a的第j个分量，(|[a]_j|-μt)₊为合页函数，sgn([a]_j)为符号函数，θ_k为第k次迭代运算引入的辅助变量，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测。

第二步：对合页函数(|[a]_j|-μt)₊进行求解，具体为：

其中：

(|[a]_j|-μt)₊为合页函数，[a]_j为a的第j个分量，μ为正则化超参数，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测。

第三步：对符号函数sgn([a]_j)进行求解，具体为：

其中：sgn([a]_j)为符号函数，[a]_j为a的第j个分量。

步骤S5.3.3：根据第k次迭代时核函数系数的估计值β_k、第k-1次迭代时核函数系数的估计值β_k-1以及预设阈值的大小，当三者关系满足如下公式时，结束迭代计算，所述公式具体为：

||β_k-β_k-1||_∞≤ε

其中：β_k为第k次迭代时核函数系数的估计值，β_k-1为第k-1次迭代时核函数系数的估计值，ε为预设阈值。

在本实施例中，预设阈值ε的大小选择为10^-6。

步骤S5.3.4：当迭代计算结束后，根据确定出的第k次迭代时核函数系数的估计值β_k，确定出参数向量的估计值

的大小，即参数向量的估计值

的大小与第k次迭代时核函数系数的估计值β_k的大小相等。

步骤S6：对待估参数求解过程中产生的正则化超参数通过10折交叉检核法进行优选，确定出优选的正则化超参数，并将优选的正则化超参数所对应的待估参数作为求解后的待估参数，具体如下：

步骤S6.1：工作人员根据实际使用状况及数据预先设置M个正则化超参数μ，并从预先设置的M个正则化超参数中选取一个正则化超参数。

步骤S6.2：将所有多源重力观测数据通过无放回随机采样方法均分为10组，并从均分的10组数据中，依次选出一组数据作为检核组的数据，并根据剩余9组数据重复步骤S1-步骤S5，建立所述残余位的球面径向基函数模型。

步骤S6.3：通过步骤S6.2建立的残余位的球面径向基函数模型对检核组中的数据进行预测，获取得到检核组的预测值。

同时根据检核组的实际观测值，将实际观测值与预测值相减得到预测误差，同时计算得到预测误差的平方和，并将10组数据分别作为检核组时所对应的预测误差的平方和进行相加，获取得到交叉检核误差，其中该交叉检核误差即为选取的正则化超参数所对应的交叉检核误差。

步骤S6.4：除已选取过的正则化超参数外，从预先设置的M个正则化超参数中重新选取一个正则化超参数，并重复步骤S6.2-步骤S6.3，直至M个正则化超参数全部被选取。

步骤S6.6：从M个正则化超参数各自对应的交叉检核误差中选出最小值，其中最小交叉检核误差对应的正则化超参数即为优选的正则化超参数。

步骤S7：将确定出的求解后的参数向量β代入步骤S3中的残余位的球面径向基函数模型中，并根据代入了参数向量β后的残余位的球面径向基函数模型恢复残余重力观测数据中的的长波部分和/或短波部分。具体地讲，将代入了参数向量β后的残余位的球面径向基函数模型和参考重力场长波模型相加，即可恢复残余重力观测数据中的长波部分，将将代入了参数向量β后的残余位的球面径向基函数模型和地形数据短波模型相加，即可恢复残余重力观测数据中的短波部分。

通过残余重力观测数据中的的长波部分和/或短波部分的回复，即可得到扰动重力位模型，从而可以得到其他重力场泛函模型，譬如：重力异常、重力扰动、垂线偏差、大地水准面、重力梯度。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。

Claims (9)

1.一种稀疏球面径向基函数局部重力场建模方法，其特征在于，所述建模方法具体包括如下步骤：

S1：搜集建模区域内的多源重力观测数据；

S2：将所述多源重力观测数据中的长波部分和/或短波部分移去，获取得到残余重力观测数据；

S3：通过所述残余重力观测数据，构建残余位的球面径向基函数模型，用以表示所述残余重力观测数据中不同位置处的残余位的大小，所述残余位的球面径向基函数模型具体为：

其中：T_res(x)为x位置处的残余位的大小，N为核函数的个数，β_j为第j个核函数的系数，x为残余位的位置，z_j为第j个核函数的位置；

S4：根据重力异常与重力扰动位之间的泛函关系和所述残余位的球面径向基函数模型，确定出各观测量的观测方程；

S5：通过稀疏正则化对所述观测方程中的参数向量进行求解；

S6：对所述待估参数求解过程中产生的正则化超参数进行优选，确定出优选的正则化超参数，并将所述优选的正则化超参数所对应的待估参数作为求解后的待估参数；

S7：将所述求解后的待估参数代入残余位的球面径向基函数模型中，并恢复所述残余重力观测数据中的长波部分和/或短波部分，得到扰动重力位模型，并导出其他重力场泛函模型。

2.根据权利要求1所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S2中，根据参考重力场模型将所述多源重力观测数据中的长波部分移去，根据地形数据将所述多源重力观测数据中的短波部分移去。

3.根据权利要求1或2所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S3中，所述第j个核函数的位置z_j的大小通过多层方法进行选定，具体为：

任意一个所述观测量的位置均对应有两个核函数的位置，所述两个核函数的位置对应有两个不同的深度，且所述两个核函数的位置中的纬度和经度与观测量的位置纬度和经度相同，所述选定的两个核函数位置，具体为：

其中：z₁为第一个选定的核函数的位置，z₂为第二个选定的核函数的位置，R_j为第j个核函数的位置所对应的地心向径，

为第j个残余位的位置所对应的纬度，λ_j为第j个残余位的位置所对应的经度，R₁为第一个选定的核函数对应的地心向径，R₂为第二个选定的核函数对应的地心向径。

4.根据权利要求3所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S4中，确定出各观测量的观测方程，具体如下：

S4.1：根据所述重力异常与重力扰动位之间的泛函关系，由所述残余位的球面径向基函数模型获取残余重力观测数据中重力异常观测的表达式，具体为：

其中：Δg(x_i)为残余重力观测数据中第i个重力异常观测，N为核函数的个数，B_ij为观测矩阵B中的第i行、第j列元素，β_j为第j个核函数的系数；

S4.2：通过所述残余重力观测数据中重力异常观测的表达式，确定出各观测量的观测方程，具体为：

y＝Bβ+e

其中：y为Δg(x_i)堆栈而成的观测向量，B为观测矩阵，β为β_j堆栈而成的参数向量，e为观测误差，Δg(x_i)为残余重力观测数据中第i个重力异常观测，β_j为第j个核函数的系数。

5.根据权利要求4所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S4.1中，所述观测矩阵B中的第i行、第j列元素B_ij，具体为：

其中：

B_ij为观测矩阵B中的第i行、第j列元素，r_i为第i个观测位置x_i所对应的地心向径，R_j为第j个核函数的位置z_j所对应的地心向径，z_j为第j个核函数的位置，x_i为第i个重力异常观测所在的位置。

6.根据权利要求4所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S5中，通过稀疏正则化对所述观测方程中的参数向量进行求解，具体如下：

S5.1：通过各观测量的观测方程，确定协方差矩阵的大小，具体为：

Q＝cov[e]

其中：Q为协方差矩阵，e为观测向量的观测误差；

S5.2：根据协方差矩阵的大小，对所述观测方程中的参数向量进行基追踪原子分解，具体为：

其中：

为参数向量的估计值，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，B为观测矩阵，β为β_j堆栈而成的参数向量，β_j为第j个核函数的系数，Q为协方差矩阵，μ为正则化超参数；

S5.3：通过FISTA方法对所述参数向量的估计值进行求解，确定出求解后的参数向量。

7.根据权利要求6所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S5.3中，通过FISTA方法对所述参数向量的估计值进行求解，具体如下：

S5.3.1：将所述核函数的系数进行初始化，具体为：

其中：θ₁和s₁均为第一次迭代运算引入的辅助变量，β₀为迭代前核函数系数的初步估计值，B为观测矩阵，Q为协方差矩阵，μ为正则化超参数，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测，I为单位矩阵；

S5.3.2：根据所述迭代前核函数系数的初步估计值，进行迭代计算，确定出第k次迭代时核函数系数的估计值，具体为：

其中：β_k为第k次迭代时核函数系数的估计值，

为软阈值算子，θ_k+1和s_k+1均为第k+1次迭代运算引入的辅助变量，β_k-1为第k-1次迭代时核函数系数的估计值，s_k为第k次迭代运算引入的辅助变量；

S5.3.3：根据所述第k次迭代时核函数系数的估计值、第k-1次迭代时核函数系数的估计值以及预设阈值的大小，当三者关系满足如下公式时，结束迭代计算，所述公式具体为：

||β_k-β_k-1||_∞≤ε

其中：β_k为第k次迭代时核函数系数的估计值，β_k-1为第k-1次迭代时核函数系数的估计值，ε为预设阈值；

S5.3.4：将所述第k次迭代时核函数系数的估计值作为参数向量的估计值的大小。

8.根据权利要求7所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S5.3.2中，所述软阈值算子的求解过程，具体如下：

第一步：给定一个向量a，对所述向量进行软阈值操作，获取得到一个同维数的向量，所述向量的第j个元素，具体为：

其中：

a表示软阈值算子中的向量，[a]_j为a的第j个分量，(|[a]_j|-μt)₊为合页函数，sgn([a]_j)为符号函数，θ_k为第k次迭代运算引入的辅助变量，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测；

第二步：对所述合页函数进行求解，具体为：

其中：

(|[a]_j|-μt)₊为合页函数，[a]_j为a的第j个分量，μ为正则化超参数，λ_max为最大特征值，B为观测矩阵，Q为协方差矩阵，y为Δg(x_i)堆栈而成的观测向量，Δg(x_i)为残余重力观测数据中第i个重力异常观测；

第三步：对所述符号函数进行求解，具体为：

其中：sgn([a]_j)为符号函数，[a]_j为a的第j个分量。

9.根据权利要求6所述的一种稀疏球面径向基函数局部重力场建模方法，其特征在于，在所述步骤S6中，对所述待估参数求解过程中产生的正则化超参数进行优选，具体如下：

S6.1：从预先设置的M个正则化超参数中选取一个正则化超参数；

S6.2：将所有所述多源重力观测数据通过无放回随机采样方法均分为10组，并从所述10组数据中，依次选出一组数据作为检核组的数据，并根据剩余9组数据重复步骤S1-步骤S5，建立所述残余位的球面径向基函数模型；

S6.3：通过所述残余位的球面径向基函数模型对检核组中的数据进行预测，获取所述检核组的预测值，根据所述检核组的实际观测值，将所述实际观测值与预测值相减得到预测误差，同时计算所述预测误差的平方和，并将所述10组数据分别作为检核组时所对应的预测误差的平方和相加，获取得到所述选取的正则化超参数所对应的交叉检核误差；

S6.4：除所述选取的正则化超参数外，从所述预先设置的M个正则化超参数中重新选取一个正则化超参数，并重复步骤S6.2-步骤S6.3，直至所述M个正则化超参数全部被选取；

S6.6：从所述M个正则化超参数各自对应的交叉检核误差中选出最小值，所述最小交叉检核误差对应的正则化超参数即为优选的正则化超参数。