利用熵函數模型解算GPS高程異常

鄧捷利,陳天偉,楊茜茜,鄭旭東

(桂林理工大學 a.廣西空間信息與測繪重點實驗室;b.測繪地理信息學院,廣西 桂林 541006)

GPS定位技術被廣泛應用于測量工作之中,GPS測量獲取的是大地高程,但工程中通用的高程系統為正常高,可采用水準測量方法獲得。大地高與正常高之間存在差值,稱為高程異常,即H=h+ξ, 其中，H表示大地高,h為正常高,ξ為高程異常。若能獲取某點的高程異常就能將該點的GPS高程應用于工程之中,所以確定高程異常在實際工程應用中尤為重要。

獲取高程異常有多種方法,其中一種是在進行GPS測量的同時,再對部分點進行水準測量,就可得知這部分點的高程異常,通過這些已知高程異常的點,就能插值或擬合出測區內其他點的高程異常[1]。不同擬合方法適用條件不同,控制點呈直線排列的情形時(如公路、 鐵路等),使用曲線擬合方法;當控制點呈面狀分布且高程異常變化不明顯才適合使用平面擬合;地形復雜則使用曲面擬合；還有一些常用擬合模型,加權平均值法(多用于山區地形)、曲面樣條法、神經網絡法等[2],本文主要討論通過加權均值法給出權系數初值再使用最大熵法插值。

待定點所使用的已知點是離散的、隨機的,所以插值得到高程異常變量也有著同樣的性質。為使插值模型使用的數據充分離散,讓模型更穩健,使用最大熵法。依據最大熵原理[3],當系統熵值最大,模型誤差將會最小。

1 加權均值法原理

加權均值法是先對已知值進行定權，然后通過加權平均來求得期望值：

(1)

該插值方法的核函數是具有地理意義的,表達了點間距對高程異常插值的影響。核函數具有地理意義的還有克里金法、楊赤中法、分形插值等,這些插值方法的權系數也適合作為最大熵法的權系數初值。

2 最大熵原理

最大熵理論的數學模型為

(2)

其中:S[p]為系統的熵值;pi為該方法求得的權系數，pi≥0;gj表示各階統計矩函數； E(gj)表示統計矩的期望;i=1,2,…n;j=1,2,…,m。

對于待定點而言,聯測點(已知高程異常)的分布是離散的、隨機的,為了能反映出空間屬性,則該最大熵模型的約束條件應該包括待定點高程異常點的數學期望值、方差。

模型系統的熵值:

S[p]=-(p1lnp1+p2lnp2+…+pnlnpn)/ln 2;

(3)

約束條件為

(4)

以上約束條件從上到下分別為0、1、2階矩約束。E1=E(ξ)是高程異常的期望值,因為隨機變量的原點矩等于相應的樣本矩[5],可由加權均值法算得。E2可以通過高程異常樣本值和E1按計算方差估值的公式算得,此式即為二階矩約束條件。本文的目的是在約束條件下使得系統的熵值S[p]趨向于最大,使得模型最大程度地依賴樣本,則模型更穩健。

待定點的高程異常值ξ的估計量為

(5)

3 解算方法

最大熵模型是個附有限制條件的極值問題, 由于該函數是非凸函數, 難以使用拉格朗日乘數法求解, 多用數值分析進行求解。 因為傳統的優化算法存在一些缺點, 對于某些目標函數不易求得最優解。

遺傳算法是一種新興算法,該算法解算快并且可以求出全局最優解。經過多年的研究,該方法因其多種優點得到了廣泛的認同。遺傳算法(GA)借鑒遺傳學中自然選擇進化過程來求解答案,該方法大致可以分為4個步驟:編碼、初始群體的生成、適應度價值評估檢查、操作算子[6]。下文具體論述初始群體生成中的初值確定和適應度價值評估檢查中如何處理約束的問題。

3.1 初值的確定

遺傳算法的尋優性能會受到初值設置的影響,初值設置十分重要。本文熵函數待求值pi是樣本點(已知點)的權系數,而插值的方法多種多樣,權系數的確定方法也復雜多樣,基于上文加權均值法的優點,采用加權均值法所得到的結果作為一階約束的樣本矩，以反距離定權作為初值。

3.2 約束處理

遺傳算法通常解決無約束問題,而的最大熵模型則包括了多個限制條件,因為沒有限制條件的最大熵是沒有意義的[7]。工程中常使用罰函數法來處理約束問題。該方法是通過一個因子把目標函數與約束條件鏈接成新的無約束的目標函數,構造方法一般有兩種:一種是采用加法形式val(x)=f(x)+c(x); 另一種是乘法形式val(x)=f(x)·c(x), 其中c(x)是懲罰項。當解在非可行域,就會進行懲戒,目標函數就會增大；當解可行,懲罰值就低,目標函數就變小,通過迭代使得目標函數最小,從而求出最優解。

4 實例分析

一般來說, 地形平緩區域的建模精度要優于地形復雜區域, 故實例采用地形復雜的數據。 本文實驗數據來源于黎劍[8]對區域GPS高程異常擬合及建模方法的研究, 對象為云南某地區, 該地區為山區地形,共27個GPS水準聯測點(大地高和正常高已知), 剔除3個點, 將16個點作為樣本點, 8個點作為待定點。 點分布如圖1所示。

大地高和正常高相減可得到8個待定的高程異常, 記錄在表1的第2列。 參照文獻[8]中使用二次曲面擬合以及多面函數擬合得到的值記錄在表1的3、 4列。 通過第1節所描述的方法即可得到某個待定點周圍已知點的權系數,通過周圍已知點權系數以及高程異常即可獲取該待定點高程異常的插值。在搜索方式采用鄰域搜索圓,搜索點數為4個的條件下使用加權均值法,得到的插值記錄在表1第5列。

圖1 24個點點位分布圖Fig.1 Distribution view of the 24 points

表1 高程異常以及4種方法的插值結果

Table 1 Elevation anomalies and interpolation results of four methods m

點號真值二次曲面多面函數加權均值最大熵法431.97231.88631.79831.94131.927732.31232.28632.29232.22732.2691032.64932.43232.51932.55332.5331232.08832.20432.15132.12332.1051632.27832.2732.29832.22832.2571832.02932.10932.03632.08532.082132.47232.43732.49432.45232.4582332.17232.20232.20432.16632.153

在Matlab平臺利用遺傳算法函數進行編程,對最大熵模型進行解算,對該模型的解算本質就是附有限制條件的極值問題,因為函數本身性質,只能通過數值分析的方法進行求解。解算大致分為3個步驟:首先,將式(3)作為目標函數,將該函數編制函數文件sangobj.m;然后,將式(4)的約束條件編為罰函數文件sangcstr.m,其中邊界約束pi>0在主程序中列出;最后,設置初值,如本文第3節描述,使用加權均值法得到的權系數作為種群初值。由以上步驟解算得各pi值為權系數。通過以上最大熵法獲得的權系數和對應已知點的高程異常則可得到待定點的插值,結果見表1第6列。通過真值和擬合值之差可以得到殘差,使用殘差進行精度分析，結果見表2。

表2 各模型精度分析

Table 2 Precision analysis of each model cm

二次曲面多面函數加權均值法最大熵法最大殘差21.717.49.611.6最小殘差0.820.6-1.7中誤差9.88.25.34.9

加權均值法與前兩種方法對比得出:由于加權均值法中反距離加權法的核函數具有地理意義,更能反映地形變化對高程異常的影響,并使用搜索半徑進行分段擬合,所以加權均值法在地形較為復雜的地區比起前兩種方法擬合精度要更高一些。加權均值法在7號點擬合表現差于前兩種方法,是因為7號點位于圖幅邊緣,該點并未被已知點包圍,已知點分布不均勻,所以導致插值精度降低。16號點因為該點附近的樣本點高程異常都比較相近,可推測出16號附近的高程異常變化不大,所以各種擬合方法差別不大。21號點各擬合方法精度也比較接近,與16號點相似,該點周圍高程異常差變化不明顯,且該點周圍的4個已知點接近對稱分布,對該點的擬合如同對相近的數值求均值,所以各種擬合方法差別很小。

對比加權均值法和最大熵法得出:使用最大熵模型結合遺傳算法得出的插值數據精度大多高于加權均值法,雖然在最大殘差、最小殘差上最大熵法并未表現出優越性,但最大熵法中誤差較低,比起加權均值法,其整體更接近真值,只有在少部分點表現出較差的結果,如23號點改進后比改進前要差,因為該點附近的已知點分布不均勻,多位于該點上方,該點的搜索范圍下方缺少已知點信息，由于對樣本依賴程度較大,插值結果偏向下方已知點的信息，精度反而變低。

5 結束語

當已知點分布較為均勻時,與核函數沒有地理意義的二次曲面擬合法以及多面函數法等相比,采用加權平均值法,即使山區地形大地水準面變化較大,也能較好地解決高程異常擬合問題。

最大熵法最大程度地排除了主觀因素的影響,讓模型最大程度地依賴樣本數據,使得模型更加穩健。與加權平均值法相比，最大熵法擬合出的數據整體上更為靠近真值,證明在地形復雜的區域使用最大熵法進行擬合有著良好的表現。