劉曉剛，孫中苗，管 斌，范昊鵬,4

1. 地理信息工程國家重點實驗室，陜西 西安 710054； 2. 西安測繪研究所，陜西 西安 710054； 3. 大地測量與地球動力學國家重點實驗室，湖北 武漢 430077； 4. 信息工程大學地理空間信息學院，河南 鄭州 450052

航空重力測量技術，因其可以在沙漠、沼澤、冰川、原始森林、陸海交界等一些難以開展地面重力測量的區(qū)域進行作業(yè)，快速經(jīng)濟地獲取精度良好、分布均勻、大面積的地球重力場中高頻信息，從而成為地球重力場研究的最為熱門的領域之一[1-2]。

航空重力測量技術主要包括標量測量技術和矢量測量技術。航空重力標量測量技術已經(jīng)非常成熟，而航空重力矢量測量技術正在國家高分辨率對地觀測系統(tǒng)重大專項的支持下開展研究。我國的航空重力矢量測量儀初樣機已經(jīng)研制成功，其集成試驗在經(jīng)歷了飛機選型論證及加改裝、試驗區(qū)勘選及測量、系統(tǒng)安裝集成及測試標定等準備工作后，2015年6—10月，西安測繪研究所組織國防科技大學、第一測繪導航基地等單位在內(nèi)蒙古地區(qū)進行了我國首次航空重力矢量測量儀集成飛行試驗，獲得了第一手的航空和地面重力矢量測量數(shù)據(jù)。2017年6月，西安測繪研究所組織上述單位在山西地區(qū)進行了我國第二次航空重力矢量測量儀集成飛行試驗。因此，這些原始測量數(shù)據(jù)的預處理及其專業(yè)處理的成果，是評判我國首套航空重力矢量測量儀系統(tǒng)的研制是否滿足要求的重要標準，而向下延拓則是其專業(yè)處理中非常重要的一個環(huán)節(jié)。

航空重力測量數(shù)據(jù)向下延拓結(jié)果的好壞，直接影響到其進一步應用，如不同類型重力測量數(shù)據(jù)的融合、全球或區(qū)域地球重力場模型的構(gòu)建、全球或區(qū)域(似)大地水準面的精化、水下重力匹配輔助導航中重力基準圖的生成等。將空中重力測量數(shù)據(jù)向下延拓到地球表面采用的是逆Poisson積分方程，屬于典型的不適定問題[3-6]。為了減弱方程的病態(tài)性對于延拓結(jié)果的影響，國內(nèi)外學者提出了很多解決思路，如最小二乘配置法[7-11]、直接代表法[12]、迭代法[13]、解析法[14-16]和正則化法[17-26]等。

目前使用的大多數(shù)向下延拓方法都是在Poisson積分方程的基礎上發(fā)展而來的，而局部區(qū)域空中格網(wǎng)重力測量數(shù)據(jù)的向下延拓，需要先將Poisson積分進行離散化處理，從而達到離散求和的目的。本文在對傳統(tǒng)向下延拓模型進行數(shù)值分析的基礎上，提出了Poisson積分迭代法這一延拓思路，并與傳統(tǒng)最小二乘法、改進最小二乘法、Tikhonov正則化法等延拓模型進行了精度比較。

1 傳統(tǒng)最小二乘法延拓模型

根據(jù)Poisson積分公式，將地面重力異常數(shù)據(jù)向上延拓，有[7]

(1)

由于在實際情況下空中和地面重力異常均是離散值，式(1)可寫成如下形式的觀測方程

gh=AX

(2)

式中，gh為空中重力異常Δgh所組成的向量；X為地面區(qū)域的重力異常Δg0所組成的向量；A為相應的系數(shù)陣。在航空重力測量中，空中重力異常Δgh是已知觀測值，地面重力異常為待估參數(shù)，這是Poisson積分方程的逆問題。

式(2)的傳統(tǒng)離散化形式為

(3)

式中，矩陣Aij的表達式為

(4)

對系數(shù)矩陣A作如下奇異值分解

A=UΛVT

(5)

式中，U、V分別由矩陣A的左、右特征向量組成；U=[u1u2…un]，V=[v1v2…vn]，U∈Rn×n,V∈Rn×n，且UTU=UUT=I，VTV=VVT=I；Λ∈Rn×n為對角陣，其對角元素是矩陣A的奇異值λ1、λ2、…、λn，且λi按遞減的次序排列。

考慮到式(5)，式(2)的最小二乘解及其譜分解形式為[19]

(6)

式中，λi為A的奇異值；ui、vi分別為其相應的左右奇異向量；rank(A)表示矩陣A的秩；而l=gh為誤差方程的自由項。

式(6)即為傳統(tǒng)的最小二乘法延拓模型。

2 改進最小二乘法延拓模型

采用式(4)所示的Poisson積分離散化形式在某些情況下不能得到正確的結(jié)果，甚至誤差很大。主要原因是，與實際結(jié)果相比，處于計算點向徑方向的流動點對計算點的貢獻較大，如果以該點的核函數(shù)值作為整個格網(wǎng)的核函數(shù)值，將會帶來較大誤差。因此，有學者給出了修改向徑方向流動點加權(quán)系數(shù)的離散化形式。矩陣A的對角線元素為[27-30]

(7)

矩陣A的非對角線元素為

(8)

式中，ψ0表示選取積分半徑；NC表示積分半徑內(nèi)的測量點個數(shù)。

該離散化公式的推導思路是先將積分半徑內(nèi)的其他格網(wǎng)對計算點的延拓貢獻值求出，再將整個積分范圍內(nèi)的理論貢獻值與之相減，即得到向徑方向的流動點對計算點的貢獻。

將式(6)中的矩陣A用式(7)和式(8)來替換，就可以得到改進的最小二乘法延拓模型。

3 Tikhonov正則化法延拓模型

航空重力測量數(shù)據(jù)向下延拓時，觀測方程是病態(tài)的，其系數(shù)矩陣A的奇異值單調(diào)地趨向零，由式(6)難以獲得穩(wěn)定的解。因此，需對病態(tài)方程進行正則化處理，以抑制高頻測量噪聲對延拓結(jié)果的影響。在眾多正則化方法中，以Tikhonov正則化法的應用最為廣泛，所采用的準則如下

(9)

根據(jù)式(9)的約束條件，可得Tikhonov正則化解為

(10)

Tikhonov正則化解是有偏的，其偏差為

(11)

Tikhonov正則化解的平均均方誤差為

MSE(Xα)=D(Xα)+bias(Xα)bias(Xα)T=

(12)

將式(5)代入式(10)、(12)，可得Tikhonov正則化解的譜分解式為[19]

(13)

通過實際計算表明，采用式(13)來獲得Tikhonov正則化解的速度要優(yōu)于式(10)，這是因為譜分解式避免了許多大型矩陣的乘法和求逆運算，節(jié)約了大量計算時間。

(14)

上式表明，Tikhonov正則化解的誤差由兩部分組成：第一項為觀測誤差所引起的估值誤差，隨α的增大單調(diào)減??；第二項為正則化所引起的估值誤差，隨α的增大單調(diào)增大，當α=0時，該項誤差為0。選擇合適的正則化參數(shù)α可起到平衡這兩項誤差的作用，使得它們的和為最小，這正是通過正則化算法可獲得精確、穩(wěn)定解的原因。

4 改進Poisson積分迭代法延拓模型

如何正確地確定正則化參數(shù)，則是整個正則化算法的關鍵，參數(shù)選擇的優(yōu)劣直接影響了最終延拓結(jié)果的精度。在模擬試驗中，最佳正則化參數(shù)的確定，是在延拓面數(shù)據(jù)已知的情況下根據(jù)最小延拓誤差來估算的。而對于實際的航空重力測量數(shù)據(jù)向下延拓，延拓面并不一定存在已知的測量數(shù)據(jù)。因此，在這種情況下如何確定正則化參數(shù)，存在較大的難度。為了減弱延拓模型的病態(tài)性，提高延拓結(jié)果的精度，本文提出Poisson積分迭代法這一延拓思路，進而建立相應的延拓模型。

由于R/r也是調(diào)和函數(shù)，因此有

(15)

在式(5)兩端同時乘以Δg0(R,θ,λ)，有

(16)

將式(1)減去式(16)，通過變換可以得到

(17)

對式(17)進行改化，可以得到

(18)

由于式(18)兩端都含有Δg0(R,θ,λ)項，因此，該式可以進一步改化成迭代形式。作為第一次近似，命

(19)

再用下式進行迭代計算

sinθ′dθ′dλ′

(20)

式(20)則為本文建立的Poisson積分迭代法向下延拓模型，該式可以進一步修改成矩陣計算的形式

(21)

式中，矩陣A可以用式(4)來計算(本文稱之為新建模型Ⅰ，即Poisson積分迭代法向下延拓模型)，也可以用式(7)和式(8)來計算(新建模型Ⅱ，即改進Poisson積分迭代法向下延拓模型)。當前后兩次計算值之間的差異小于ε時，就可以終止迭代。

5 數(shù)值試驗與結(jié)果分析

5.1 試驗數(shù)據(jù)說明

采用中國某地區(qū)實測航空重力測量數(shù)據(jù)，分辨率是5′×5′，飛行高度為3400 m。其等值線圖如圖1所示。

相應的地面重力測量數(shù)據(jù)，其等值線圖如圖2 所示。

航空和地面重力測量數(shù)據(jù)的統(tǒng)計結(jié)果如表1所示。

表1 數(shù)據(jù)統(tǒng)計結(jié)果

圖1 航空重力測量數(shù)據(jù)等值線圖(單位：mGal)Fig.1 Isoline chart of airborne gravimetry data (Unit：mGal)

圖2 地面重力測量數(shù)據(jù)等值線圖(單位：mGal)Fig.2 Isoline chart of ground gravimetry data (Unit：mGal)

根據(jù)數(shù)字高程模型SRTM30，采用Global Mapper軟件獲得了該地區(qū)5′×5′的地形數(shù)據(jù)，其等值線圖如圖3所示。

圖3 地形數(shù)據(jù)等值線圖(單位：m)Fig.3 Isoline chart of topography data(Unit：m)

地面重力測量數(shù)據(jù)對應的地形數(shù)據(jù)統(tǒng)計結(jié)果如表2所示。

表2 試驗區(qū)地形數(shù)據(jù)結(jié)果統(tǒng)計

以地面點的平均地形高度1 476.115 m作為向下延拓的基準面，則延拓高度是1 923.885 m。

根據(jù)圖1、圖2與圖3的對比可以看出，航空和地面重力數(shù)據(jù)的變化趨勢，與地形的起伏密切相關。地形起伏較為劇烈的地區(qū)，重力數(shù)據(jù)也有相應的變化。

5.2 對比試驗結(jié)果

首先，利用傳統(tǒng)的最小二乘法延拓模型，將航空重力測量數(shù)據(jù)向下延拓到地面平均高程面上，其等值線圖如圖4所示。將延拓結(jié)果與地面重力測量數(shù)據(jù)進行比較，精度統(tǒng)計結(jié)果列于表3中。

圖4 傳統(tǒng)最小二乘法延拓模型的延拓結(jié)果等值線圖(單位：mGal)Fig.4 Isoline chart of the continuation results of traditional least squares model(Unit：mGal)

其次，利用改進的最小二乘法延拓模型，將航空重力測量數(shù)據(jù)向下延拓到地面平均高程面上，其等值線圖如圖5所示。將延拓結(jié)果與地面重力測量數(shù)據(jù)進行比較，精度統(tǒng)計結(jié)果列于表3中。

再次，在使用Tikhonov正則化法延拓模型之前，需要確定其正則化參數(shù)。根據(jù)地面已知重力測量數(shù)據(jù)，計算了延拓誤差與正則化參數(shù)的關系，如圖6所示，并確定了最優(yōu)正則化參數(shù)α=0.031。采用該正則化參數(shù)，將航空重力測量數(shù)據(jù)向下延拓，延拓結(jié)果等值線圖如圖7所示，精度統(tǒng)計結(jié)果列入表3中。

圖5 改進最小二乘法延拓模型的延拓結(jié)果等值線圖(單位：mGal)Fig.5 Isoline chart of the continuation results of improved least squares model(Unit：mGal)

表3 延拓結(jié)果統(tǒng)計

圖6 延拓誤差與正則化參數(shù)關系圖Fig.6 Relationship of continuation errorsand regularization factors

最后，利用本文建立的Poisson積分迭代法延拓模型將航空重力測量數(shù)據(jù)向下延拓，延拓結(jié)果的等值線圖如圖8、圖9所示，精度統(tǒng)計結(jié)果列于表3中。

圖7 Tikhonov正則化法延拓模型的延拓結(jié)果等值線圖(單位：mGal)Fig.7 Isoline chart of the continuation results of Tikhonov regularization model(Unit：mGal)

圖8 Poisson積分迭代法延拓模型I的延拓結(jié)果等值線圖(單位：mGal)Fig.8 Isoline chart of the continuation results of Poisson integral iteration model I(Unit：mGal)

圖9 Poisson積分迭代法延拓模型Ⅱ的延拓結(jié)果等值線圖(單位：mGal)Fig.9 Isoline chart of the continuation results of Poisson integral iteration model Ⅱ(Unit：mGal)

根據(jù)表3的統(tǒng)計結(jié)果，相比較于傳統(tǒng)的最小二乘法延拓模型，改進的最小二乘法延拓模型的結(jié)果精度提高了約15 mGal，采用Tikhonov正則化法對延拓模型的病態(tài)性進行修正以后，延拓結(jié)果精度進一步提高了約0.34 mGal，本文新建立的Poisson積分迭代法延拓模型精度要優(yōu)于改進的最小二乘法延拓模型，精度提高了約0.21 mGal，但稍遜于Tikhonov正則化法延拓模型，精度略低于0.13 mGal。

通過Poisson積分迭代法延拓模型Ⅰ和Ⅱ的結(jié)果可以看出，不論采用傳統(tǒng)還是改進的Poisson積分離散化公式，延拓結(jié)果幾乎一致，這說明本文建立的Poisson積分迭代法延拓模型不僅可以修正傳統(tǒng)Poisson積分離散化誤差給延拓模型帶來的病態(tài)性影響，而且對延拓模型本身的病態(tài)性有一定的抑制作用(精度優(yōu)于改進的最小二乘法延拓模型)。

6 結(jié)束語

航空重力測量技術，因其可以快速經(jīng)濟有效地獲取分布均勻、精度良好、大面積的重力場中高頻信息，在我國陸地和陸海交界區(qū)域的重力測量中得到了很好的應用。在國內(nèi)外學者研究成果的基礎上，本文利用實測航空和地面重力數(shù)據(jù)，對傳統(tǒng)最小二乘法、改進最小二乘法、Tikhonov正則化法等延拓模型進行了數(shù)值分析，最后建立了Poisson積分迭代法和改進Poisson積分迭代法延拓模型。本文的研究成果，可應用于我國航空重力標量和矢量測量數(shù)據(jù)的處理中，從而為獲得更高精度、更高分辨率的地球重力場產(chǎn)品提供一定的技術支撐。