地形熵聚類下利用多面函數擬合高程異常

孫佳龍，郭淑艷，龍冰心，秦思遠

(1. 淮海工學院測繪與海洋信息學院，江蘇 連云港 222005； 2. 海島(礁)測繪技術國家測繪地理信息局重點實驗室，山東 青島 266590； 3. 海岸帶地理環境監測國家測繪地理信息局重點實驗室，廣東 深圳 518060； 4. 江蘇省海洋資源開發研究院，江蘇 連云港 222001； 5. 北京大地宏圖勘測科技有限公司，北京 100029)

高程異常是似大地水準面與參考橢球面之間的高程差值。在一個區域內精確確定高程異常，就可以通過GPS觀測的大地高準確地求出正常高，從而減少對水準測量工作的依賴，提高測量效率[1-3]。在一個區域內確定高程異常，通常是利用已知點的高程異常擬合高程異常曲面，從而插值得出待定點的高程異常。高程異常曲面的擬合方法主要有多項式擬合法、二次曲面擬合法和多面函數擬合法等[4-6]。多面函數法是從幾何觀點出發，解決根據數據點所形成的平差數學曲面問題，該方法主要應用于擬合重力異常、大地水準面差距、垂線偏差等。而在多面函數擬合時，首先需要確定中心點，中心點選擇的不同，在很大程度上決定了最后的擬合精度[7]。人為選擇中心點，會導致擬合的精度隨經驗的變化而不同，從而導致求解的高程異常值因人而異。聚類分析方法可以利用點與點之間的空間距離進行分類，從而在分類結果中有效地選擇出中心點，但利用該方法需首先確定分類的個數，而分類個數有時又難以確定[8-9]。本文以地形熵作為一個參考量，首先確定分類個數，然后以點與點之間的空間距離與地形熵的比值作為指標，從眾多高程異常點中選擇多面函數的中心點，并利用多面函數進行擬合，同時與其他方法進行比較，以期得出一些有益的結論。

1 多面函數

多面函數的理論依據為：任意數學表面與任意不規則的圓滑表面，總可以利用一系列有規則的數學表面總和，以任意精度逼近。因此，由多面函數表示的殘差高程異常為[4]

(1)

式中，aj為待定系數；Q(x,y,xj,yj)為x和y的二次函數，其中心點在(xj,yj)處；n為中心點的個數；ζ可由二次式的和確定，故稱多面函數。

用一組簡單的解析函數的線性組合對這個曲面進行逼近，通常選擇多元二次函數作為解析函數，其形式為[6]

(2)

將式(2)寫成誤差方程，并用矩陣表示為

V=QA-ζ

(3)

待定系數A可根據已知GPS/水準數據和地球重力場模型得到的殘差高程異常，按最小二乘法計算[10-11]

(4)

由此可根據系數A和式(1)求解其他點的殘差高程異常。

2 地形熵聚類多面函數擬合方法

當用多面函數表示殘差高程異常時，由于選擇的中心點個數n要小于已知的高程異常值的個數，因此，需要對所有的高程異常值進行選取。人為選取高程異常值既耗時又費力，且容易出錯。利用聚類分析方法選擇中心點時，需要首先確定分類個數，而對于某一區域的高程異常而言，分類個數有時難以確定。作為一種表示地形復雜程度的指標，地形熵可以較準確地表示局部區域高程值變化的劇烈程度[12-15]。地形熵定義為

(5)

式中，h(i,j)為高程，在本文中，以高程異常值代替高程。局部地形熵反映了地形所含有信息量的大小，因此，局部地形熵可以描述地形的性質。局部區域高程值變化越急劇，起伏變化越大，地形越獨特，則計算出的局部地形熵越小，否則熵就越大。高程異常值變化越劇烈的位置，該點應參與到曲面擬合中，即應作為中心點，擬合多面函數。基于地形熵聚類的多面函數擬合算法過程如下：

(1) 將相鄰兩個高程異常數據點組成新塊；

(2) 分別計算每塊區域的地形熵；

(3) 如某塊熵的相鄰兩塊區域的熵均大于該塊的熵，則該塊熵為極小熵；

(4) 極小熵的個數n即為聚類分析的分類數，利用聚類分析算法對高程異常點進行聚類；

(5) 聚類分析中各類別的距離與各塊區域地形熵的比值作為新的指標量，進行排序；

(6) 選擇指標量最大的n個點作為多面函數的中心點，其余點為檢核點。

3 實例解算與分析

某地區D級GPS控制網布設了18個觀測點，對其進行了GPS觀測和三等水準測量，利用GPS獲得的大地高和水準測量得到的正常高，得到18個GPS/水準點的高程異常。各點的相對位置及高程異常值大小如圖1所示。

圖1 高程異常點分布情況

從圖1可以看出，這些點分布不均勻，在中間區域分布較少，在高程異常值較小的區域，點位個數也偏少。對這18個點的高程異常進行分塊計算地形熵，如圖2所示。

從圖2可以看出，第2、6、9、12、15和17塊區域為極小熵，與極小熵相鄰較小的塊熵為第3、5、10、11和16，由此可以得到重合的點號分別為3、6、10、12、16和17號點，共6個點。將聚類分析的分類個數取為6，進行聚類分析，結果如圖3所示。

圖2 各塊區域地形熵

圖3 聚類分析結果

從分類結果可以看出，距離比較近的點被分為一類，如7和8號點。由于9和11號點距離其他點均較遠，因此也被分為一類。通過聚類分析，可以較好地將各點根據空間距離的遠近分成不同的類。通過聚類分析，將距離最大的6個點，即9、10、11、12、16、18作為聚類分析選擇的中心點。

無論是單獨采用地形熵還是聚類分析，兩者選擇的點是有差異的。本文將空間距離與地形熵的比值作為新的指標量，既考慮了點與點之間的平面距離，又反映了高程異常值的變化情況，通過該指標，地形熵聚類法重新選擇了中心點。利用不同方法選擇的中心點進行了多面函數的擬合，求得檢核點的高程異常，與檢核點的實際高程異常求差，結果如圖4所示。

圖4 多種方法選擇中心點后的擬合結果

從圖4可以看出，地形熵多面函數擬合的結果與實際高程異常值相比，在不同點差異較大，差異最小時接近0，而差異最大時超過0.06 m，說明地形熵多面函數在局部區域可以很好地選擇中心點，實現有效擬合，但在個別區域差值較大。而聚類分析多面函數比地形熵多面函數變化幅度小，個別點高程異常差大于0.04 m，說明聚類分析多面函數在個別點與實際高程異常的差異也較大。而地形熵聚類的多面函數變化的幅度相對最小，最大差異值為0.03 m，為3種方法中最小的。3種方法的標準差和外符合精度見表1。

表1 擬合結果殘差比較

從表1可以看出，10、12和16號點在3種方法中都被選擇，說明這3個點具有很強的代表性。而3種方法選擇的其他中心點卻有所不同，說明采用不同的指標量得到的結果具有差異性。計算不同中心點的多面函數所擬合的高程異常，然后求其與實際高程異常的差值見表1，由地形熵多面函數擬合的結果在3種方法中無論是標準差還是外符合精度都是最低的，說明單獨采用地形熵選擇的中心點還不能很好地反映整個區域的高程異常變化，特別是隨坐標的變化。而聚類分析得到精度較高，說明在高程異常變化不大的情況下，以空間距離的大小作為指標量，可以很好地獲取中心點，從而提高擬合精度。當將地形熵和聚類分析方法結合后，得到的標準差和外符合精度都是最小的，說明同時考慮坐標和高程異常的變化可以在一定程度上提高曲面擬合精度，在外符合精度上分別比聚類分析和地形熵提高了28%和48%，說明了該方法的有效性。

4 結 語

利用地形熵可以較好地表征區域高程異常變化情況，從而為聚類分析提供分類個數。通過計算空間距離與地形熵的比值可以較好地反映整個區域的平面和高程異常變化整體情況。通過3種方法選擇了多面函數的中心點，利用多面函數對整個區域的高程異常進行擬合，擬合結果顯示距離與地形熵的比值是一個可以綜合反映高程異常隨空間位置變化的指標量，利用其選擇的中心點并進行多面函數擬合的精度也最高，然而這個指標量對于一些特殊的區域是否適合還需進一步研究。