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

翟　曜， 趙東明， 常　岑， 蒲薪吉

信息工程大學地理空間信息學院，河南 鄭州，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