1.一种稀疏球面径向基函数局部重力场建模方法,其特征在于,所述建模方法具体包括如下步骤:S1:搜集建模区域内的多源重力观测数据;
S2:将所述多源重力观测数据中的长波部分和/或短波部分移去,获取得到残余重力观测数据;
S3:通过所述残余重力观测数据,构建残余位的球面径向基函数模型,用以表示所述残余重力观测数据中不同位置处的残余位的大小,所述残余位的球面径向基函数模型具体为:其中:Tres(x)为x位置处的残余位的大小,N为核函数的个数,βj为第j个核函数的系数,x为残余位的位置,zj为第j个核函数的位置;
S4:根据重力异常与重力扰动位之间的泛函关系和所述残余位的球面径向基函数模型,确定出各观测量的观测方程;
S5:通过稀疏正则化对所述观测方程中的参数向量进行求解;
S6:对所述待估参数求解过程中产生的正则化超参数进行优选,确定出优选的正则化超参数,并将所述优选的正则化超参数所对应的待估参数作为求解后的待估参数;
S7:将所述求解后的待估参数代入残余位的球面径向基函数模型中,并恢复所述残余重力观测数据中的长波部分和/或短波部分,得到扰动重力位模型,并导出其他重力场泛函模型。
2.根据权利要求1所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S2中,根据参考重力场模型将所述多源重力观测数据中的长波部分移去,根据地形数据将所述多源重力观测数据中的短波部分移去。
3.根据权利要求1或2所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S3中,所述第j个核函数的位置zj的大小通过多层方法进行选定,具体为:任意一个所述观测量的位置均对应有两个核函数的位置,所述两个核函数的位置对应有两个不同的深度,且所述两个核函数的位置中的纬度和经度与观测量的位置纬度和经度相同,所述选定的两个核函数位置,具体为:其中:z1为第一个选定的核函数的位置,z2为第二个选定的核函数的位置,Rj为第j个核函数的位置所对应的地心向径, 为第j个残余位的位置所对应的纬度,λj为第j个残余位的位置所对应的经度,R1为第一个选定的核函数对应的地心向径,R2为第二个选定的核函数对应的地心向径。
4.根据权利要求3所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S4中,确定出各观测量的观测方程,具体如下:S4.1:根据所述重力异常与重力扰动位之间的泛函关系,由所述残余位的球面径向基函数模型获取残余重力观测数据中重力异常观测的表达式,具体为:其中:Δg(xi)为残余重力观测数据中第i个重力异常观测,N为核函数的个数,Bij为观测矩阵B中的第i行、第j列元素,βj为第j个核函数的系数;
S4.2:通过所述残余重力观测数据中重力异常观测的表达式,确定出各观测量的观测方程,具体为:y=Bβ+e
其中:y为Δg(xi)堆栈而成的观测向量,B为观测矩阵,β为βj堆栈而成的参数向量,e为观测误差,Δg(xi)为残余重力观测数据中第i个重力异常观测,βj为第j个核函数的系数。
5.根据权利要求4所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S4.1中,所述观测矩阵B中的第i行、第j列元素Bij,具体为:其中:
Bij为观测矩阵B中的第i行、第j列元素,ri为第i个观测位置xi所对应的地心向径,Rj为第j个核函数的位置zj所对应的地心向径,zj为第j个核函数的位置,xi为第i个重力异常观测所在的位置。
6.根据权利要求4所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S5中,通过稀疏正则化对所述观测方程中的参数向量进行求解,具体如下:S5.1:通过各观测量的观测方程,确定协方差矩阵的大小,具体为:Q=cov[e]
其中:Q为协方差矩阵,e为观测向量的观测误差;
S5.2:根据协方差矩阵的大小,对所述观测方程中的参数向量进行基追踪原子分解,具体为:其中:为参数向量的估计值,y为Δg(xi)堆栈而成的观测向量,Δg(xi)为残余重力观测数据中第i个重力异常观测,B为观测矩阵,β为βj堆栈而成的参数向量,βj为第j个核函数的系数,Q为协方差矩阵,μ为正则化超参数;
S5.3:通过FISTA方法对所述参数向量的估计值进行求解,确定出求解后的参数向量。
7.根据权利要求6所述的一种稀疏球面径向基函数局部重力场建模方法,其特征在于,在所述步骤S5.3中,通过FISTA方法对所述参数向量的估计值进行求解,具体如下:S5.3.1:将所述核函数的系数进行初始化,具体为:
其中:θ1和s1均为第一次迭代运算引入的辅助变量,β0为迭代前核函数系数的初步估计值,B为观测矩阵,Q为协方差矩阵,μ为正则化超参数,y为Δg(xi)堆栈而成的观测向量,Δg(xi)为残余重力观测数据中第i个重力异常观测,I为单位矩阵;
S5.3.2:根据所述迭代前核函数系数的初步估计值,进行迭代计算,确定出第k次迭代时核函数系数的估计值,具体为:其中:βk为第k次迭代时核函数系数的估计值, 为软阈值算
子,θk+1和sk+1均为第k+1次迭代运算引入的辅助变量,βk-1为第k-1次迭代时核函数系数的估计值,sk为第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(xi)堆栈而成的观测向量,Δg(xi)为残余重力观测数据中第i个重力异常观测;
第二步:对所述合页函数进行求解,具体为:
其中:
(|[a]j|-μt)+为合页函数,[a]j为a的第j个分量,μ为正则化超参数,λmax为最大特征值,B为观测矩阵,Q为协方差矩阵,y为Δg(xi)堆栈而成的观测向量,Δg(xi)为残余重力观测数据中第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个正则化超参数各自对应的交叉检核误差中选出最小值,所述最小交叉检核误差对应的正则化超参数即为优选的正则化超参数。