外部扰动引力逼近的最小二乘配置法研究

翟　曜， 赵东明， 常　岑， 蒲薪吉

信息工程大学地理空间信息学院，河南 郑州，450052



外部扰动引力逼近的最小二乘配置法研究

翟曜， 赵东明， 常岑， 蒲薪吉

信息工程大学地理空间信息学院，河南 郑州，450052

本文给出了利用附加参数最小二乘配置方法计算局部空间扰动场元的公式与计算方法。文中采用不同范围的初始数据，计算了一定高度上的局部空间扰动引力垂直分量的格网值，并将其与模型计算值进行比较分析。在移去恢复低阶模型阶数分别为36阶、90阶、180阶、360阶和720阶的情况下，使用附加参数最小二乘配置方法计算了不同高度上的扰动引力垂直分量值。结果表明，当高度一定时，所求格网范围越小并越接近已知点中心，计算精度越高。其中，移去恢复720阶低阶模型得到的计算精度最高；在已知点附近的空间范围内，计算点的高度越低，得到的扰动引力垂直分量的精度越高。

“移去-恢复”；附加参数最小二乘配置；扰动引力垂直分量

1　引　言

自上世纪60年代最小二乘配置理论被提出以来，它至今一直是物理大地测量研究的热点问题之一。该理论首次将数随机性概念引入至地球重力场的研究中，这不仅是方法论的一次突破，也是地球重力场研究在理论认识上的一次突破。在早期对最小二乘配置理论的研究中，Moritz、Tscherning和Rapp做出的贡献最为突出。Moritz最早提出了最小二乘配置理论，并在1980年出版的《高等物理大地测量》一书中全面地介绍了最小二乘配置理论[1]，其中，关于附加参数最小二乘配置的阐述是本文行文的基础之一；Tscherning与Rapp于1974年提出了全球阶方差模型[2]，凭此可以推导出闭合的全球协方差函数模型，这使得利用最小二乘配置进行大范围、乃至全球的计算成为可能。在后续的研究中为提高计算效率，罗志才等人提出频域最小二乘配置法[3]，其有效避免了计算过程中超大型矩阵求逆的问题；章传银等人推导出一种综合多种类型、不同高度重力场元经验协方差函数的通用表达方法[4]，实现了局部重力场元配置计算的一体化。在针对协方差函数的研究中，陆仲连对局部重力场中的协方差函数进行了深入研究[5]；夏哲仁等人利用回归分析方法，对局部协方差函数进行了逼近[6]，揭示了重力异常与高程的线性相关特性。

本文将“移去-恢复”方法融入附加参数最小二乘配置法中，推导出利用附加参数最小二乘配置方法计算其他扰动场元的公式，并利用其计算了不同范围及高度的扰动引力垂直分量值，为计算空间扰动引力提供了一种新的方法。

2　附加参数最小二乘配置

最小二乘配置的数学模型为：

l=AX+BS+n

(1)

其中，l为观测值向量；X为系统性参数；A为系数矩阵，其表明参数X对观测值l的影响。随机量包括两部分：信号S和误差V1。利用拉格朗日乘数法，可得到最小二乘解为：

(2)

(3)

式中，Cst表示待估信号与观测信号之间的协方差阵，Ctt表示观测信号的协方差阵。

在传统已知重力异常推求其他扰动场元的计算中，一般不考虑系统参数，即A=0，则：

(4)

但当待估信号与观测信号都具有高程属性，即空间异常与高程有较强的相关性时，则必须考虑系统参数的影响。此时，我们将空间异常分为系统部分和信号部分，其中，系统部分与高程相关。对应(1)式，AX即为观测值系统部分，X为参数，且有：

(5)

A为系数，有：

(6)

其中，Hi为已知点高程。

在进行局部地区的研究时，可以把地面近似当作平面，此时r→∞，可得：

(7)

根据协方差传播定律，待估信号扰动引力径向分量δz与观测信号重力异常的协方差函数为：

Cst=C(δzP,ΔgQ)=C(-ΔgP,ΔgQ)=-C(ΔgP,ΔgQ)

(8)

为了进一步推导，我们需要重力异常推估的附加参数最小二乘配置公式，即文献[5]中的公式：

(9)

(10)

此时，将式(7)、(8)、(9)代入式(3)可得未知点扰动引力径向分量信号估值：

(11)

(12)

式(2)、(3)、(10)即为利用附加参数最小二乘配置方法推求扰动引力径向分量的计算公式。对于其它由重力异常推导出的扰动场元Y，式(10)变为：

(13)

其中，Li为重力异常与Y的泛函关系。总结起来附加参数最小二乘配置的计算步骤为:

(3)由式(11)求得未知点的扰动场元Y值。

3　“移去-恢复”技术

“移去-恢复”理论最早由Forsberg & Tsche-rning提出[7]，其为重力场的相关计算提供了分频段计算方式，即在计算时移去输入数据中包含的低频信号，然后在计算结果中恢复相应部分信号。在附加参数最小二乘配置中引入“移去-恢复”技术，可以有效地控制中、长波段的误差。其具体步骤如图1所示：

图1　“移去-恢复”技术计算流程图

4　实验分析

实验的数据来源于EGM2008地球重力场模型生成的一片水平范围为2°×2°、分辨率为5′×5′的空间重力异常格网数据。数据共包含24×24个格网点，所有点高程随机分布在0～6000m之间。

首先，在移去恢复低阶模型阶数为180阶和720阶的情况下，计算中心与已知点水平位置中心重合、高程为3000m、1.9165°×1.9165°范围内5′×5′的23×23个格网点的扰动引力垂直分量；然后，以EGM2008地球重力场模型生成该区域3000m的扰动引垂直分量作为所求值的真值，以之与配置值进行比较，结果如图2所示：

(a) 　“移去-恢复”低阶模型阶数为180阶时附加参数配置结果

(b) 　“移去-恢复”低阶模型阶数为720阶时附加参数配置结果

(c)　模型生成的结果

(d)　图(a)与(c)的差值

(e)　图(b)与(c)的差值图2　附加参数最小二乘配置值与模型真值之间的比较结果(单位：mGal)

从图2可以看出，配置结果与模型计算值的差值绝大多数在-5～5mGal之间，差值较大的点有很大一部分出现在所求区域的边缘，而中心区域的差值较小。其统计结果如表1所示：

表1不同计算范围内的附加参数最小二乘配置值与模型真值差值的统计结果

计算结果范围(块)23×2321×2119×1917×1715×1513×1311×119×9180阶Std(mGal)2.99932.45792.40022.24662.30112.25582.14752.0711720阶Std(mGal)1.47331.39891.36101.37021.38111.36771.36941.2227

进行横向比较后可以看出，计算结果范围越小、越靠近中心，计算精度越高；若在计算结果中去除最外层的点,可明显增加计算精度。为了进一步证明这一结论，增大初始值范围至2.5°×2.5°，用同样的方法计算“移去-恢复”低阶模型阶数为180阶和720阶时，2.4165°×2.4165°范围内的29×29个点值。其统计结果如表2所示：

表2增大初始值范围后附加参数最小二乘配置值与模型真值在不同范围内差值统计结果

计算结果范围(块)29×2927×2725×2523×2321×2119×1917×1715×1513×1311×119×9180阶Std(mGal)3.06172.46172.39292.45592.42572.39642.24132.27942.22532.11362.0456720阶Std(mGal)1.47211.38871.37961.39621.39201.36051.35831.36621.36601.36211.2132

表2中表现出的趋势与表1相同，但其计算 个点的精度明显高于表1中的结果、边缘仍出现了精度较差的点。结合以上的实验结果可以得到图3与图4。从中可以看出，随着计算范围缩小，计算个点的精度逐渐提高，且最外围的点较其他点精度明显较低。其主要原因是：边缘点在水平范围内四周的强相关点较少，导致其精度较低。因此，对计算结果中靠近初始值水平范围边缘的点采取剔除的处理方式。在后续的计算中发现此结论在移去恢复不同阶数、计算不同高度时仍成立。因此，在接下来的计算分析中以均已剔除边缘点为前提。

图3　“移去-恢复”低阶模型阶数为180阶时不同计算范围的附加参数最小二乘配置值与模型真值的差值

与此同时，我们可以发现图2(d)与2(a)、2(c)中的配置结果在空间上有一定的相关性，而在图2(e)中则很难发现这种相关性。这主要是因为图2(d)的结果“移去-恢复”低阶模型的阶数较低，其与模型值相减后的结果受模型的长波项影响较大，因此显现出模型长波项的局部趋势；而图2(e)的结果“移去-恢复”低阶模型的阶数较高，其与模型值相减后的结果受模型的长波项影响较小，其随机性更高，因此这种相关性很难显现出来。这也说明，当移去模型阶数越高，剩余的模型长波项局部趋势越小，随机性更高。

为了研究“移去-恢复”低阶模型的具体阶数，仍以原始的 个点为初始数据，在“移去-恢复”低阶模型阶数分别为36阶、90阶、180阶、360阶和720阶的情况下，计算不同高度的21×21个平面格网点值，以EGM2008地球重力场模型生成对应高度扰动引力垂直分量值作为真值，其差值结果统计见表3与图5：

表3“移去-恢复”不同阶数在不同高度上的附加参数最小二乘配置值与模型真值的差值std统计结果(单位：mGal)

高度(m)阶数100030005000700012000200002160阶-36阶2.30893.28634.30795.19436.93948.92462160阶-90阶1.92383.07634.17805.10846.91898.98352160阶-180阶1.25982.45793.56624.46956.15307.94112160阶-360阶1.01222.40123.50074.39666.11917.89952160阶-720阶0.67861.39202.28632.68774.09095.9825

图5　“移去-恢复”不同阶数在不同高度上的附加参数最小二乘配置值与模型真值的差值

从表3和图5可以看出，随着高度的增加，“移去-恢复”低阶模型的计算精度逐渐降低。在6000m以下其标准差尚能保持在5mGal以下，但当所求点高度远超过已知点分布的高度范围时，其精度已无法得到保证。同时，可以看到“移去-恢复”低阶模型阶数为36阶与90阶时得到的计算结果精度相当；“移去-恢复”低阶模型阶数为180阶与360阶时得到的计算结果精度相当，且精度优于“移去-恢复”36阶与90阶低阶模型得到的结果；而“移去-恢复”低阶模型阶数为720阶时得到的结果相对于其他结果是最优的。经过计算，当“移去-恢复”低阶模型阶数继续增加后，精度提升十分有限；并且由于文中实验数据来源于模型生成，因此当“移去-恢复”低阶模型阶数过高时，实验结果将无法体现出附加参数最小二乘配置的效果。所以“移去-恢复”低阶模型阶数以720阶为宜。

5　结　论

通过实验得知：

(1)高度一定时，所求格网范围越小并越接近已知点中心，附加参数最小二乘配置的计算精度越高。当计算点水平位置位于初始值分布范围边缘或更远的地方时，其计算结果的精度难以得到保证。

(2)随着计算高度的增加，附加参数最小二乘配置的计算精度将逐渐降低；当所求点的高度远大于已知点分布的高度范围时，其计算结果的精度将无法得到保证。

(3)在“移去-恢复”低阶模型阶数的选择上，相对于36阶、90阶、180阶和360阶，“移去-恢复”720阶低阶模型得到的效果更好。

总之，基于“移去-恢复”的附加参数最小二乘配置方法可有效计算初始值附近区域的空间扰动引力垂直分量值。

[1]赫尔墨特·莫里兹.高等物理大地测量学[M].北京:测绘出版社,1984.

[2]C.C.Tscherning and R.H.Rapp. Closed Covariance Expressions for Gravity Anomalies, Geoid Deflections of the Vertical Implied by Anomaly Degree Variance Model [R]. Report No.208, Ohio: The Ohio State University.1974.

[3]罗志才,宁津生,晁定波.卫星重力梯度向下延拓的频域最小二乘配置法[J].测绘科学技术学报, 1996(3):161-165.

[4]章传银,丁剑,晁定波.局部重力场最小二乘配置通用表示技术[J].武汉大学学报·信息科学版,2007(5):431-434.

[5]陆仲连.关于在局部扰动重力场中最小二乘配置法的应用[J].测绘科学技术学报,1984(00):3-16.

[6]夏哲仁,林丽.局部重力异常协方差函数逼近[J].测绘学报,1995(1):23-27.

[7]Forsberg, R. and C.C.Tscherning: The Use of Height Data in Gravity Field Approximation by Collocation[J]. J.Geophys.Res, 1981,86(B9):7843-7854.

Research on the Least Squares Collocation Method for External Disturbance Gravity

Zhai Yao， Zhao Dongming， Chang Cen， Pu Xinji

Institute of Geospatial Information, Information Engineering University, Zhengzhou 450052, China

The formula and calculation method of local space perturbation field elements are put forward using least squares collocation with extra parameters. Initial data of different ranges are used to calculate local spatial grid values of vertical component of disturbing gravity at certain height and are compared with model results. Under the conditions that the remove-restore degree of the low degree model are 36 degree, 90 degree, 180 degree, 360 degree and 720 degree, the vertical component values of disturbing gravity at different heights are calculated respectively using least squares collocation with extra parameters. The results show that a higher calculation accuracy will be achieved as the grid range decreases and comes nearer to the center of unknown point. To be specific, the highest accuracy is obtained by the remove-restore low degree model of 720 degree, and around the known point, the lower the calculated point is, the higher accuracy will be achieved.

remove-restore; the least squares collocation with parameters; vertical component of disturbing gravity

2015-11-12。

地理信息工程国家重点实验室开放研究基金资助项目(SKLGIE 2013-Z-1-1)。

翟曜(1991—)，男，硕士研究生，主要从事物理大地测量研究。

P228.1

A