IGG抗差估計在高程擬合中的應用研究

李廣來

（華北理工大學礦業工程學院，河北 唐山 063210）

在穩健估計中，是通過選用不同的權函數來排除粗差。權函數是穩健估計的核心問題，常用的權函數有多種，如Huber函數、殘差絕對和最小函數、Tukey函數、Hampel函數、IGGI、IGGⅢ等。其中IGG方案是我國最早的權函數，對測量數據的抗差估計有較好的效果[1]。

一般使用的高程數據是正常高程，而在工程測量中GPS獲取的高程一般是WGS-84坐標下的大地高程，兩者之間存在高程異常。所以要通過擬合的方式計算高程異常，從而換算成正常高進行使用。又由于數據采集過程中會受到周圍建筑物、衛星信號等影響，測量數據難免會混入粗差。為此，筆者選用IGGI、IGGⅢ對含有粗差的GPS靜態數據和水準數據進行高程擬合，探討IGG方案對粗差的抗性。

1 IGG方案

IGG方案是基于測量誤差的有界性提出來的，它對測量抗差估計比較有效[2]。IGG屬于之間有淘汰區的M估計，權因子之間的變化較為平緩，因此K的取值和余差值小的變化影響不大。這種估計方案充分考慮了測量數據實際情況，是一種適合處理測量誤差的抗差方案。IGG方案把數據分為三部分。1）當觀測值無粗差時，用最小二乘即可得到最優解；2）當測值超出一定的范圍時，采用降權估計；3）當個別觀測值明顯異常時，采用淘汰法估計（零權估計）。

1.1 IGGI方案

1.2 IGGⅢ方案

我國最早出現的相關等價權函數，是楊元喜在周江文IGG方案下擴充提出的IGGⅢ方案[1]。

在IGGⅢ方案中關鍵在于權函數的選取，不同權函數對IGGⅢ的影響如圖1所示。

圖1 權函數曲線圖

1.3 選權迭代法

在穩健估計處理粗差過程中，關鍵在于穩健權陣的選擇，一般應用最多的方法是選權迭代法[3]。選權迭代一般就是改變權重然后不斷迭代，在每次迭代平差中，改變的是使含有粗差的觀測值的權逐次減少，最后把粗差的權趨于零來達到剔除粗差的目的[4]。選權迭代也是MATLAB進行仿真實驗的關鍵思路，選權迭代的一般處理過程如下[5]。

1）列誤差方程，令各權因子初值為1，即令w1=w2=…wn=1（w為穩健權陣），p=I（p為觀測權陣）[2]。

3）重復1）、2）構造新的穩健權陣，依次迭代直到前后每個改正數的差值小于規定的閾值，迭代終止。

2 技術設計流程圖

本研究的技術路線大致為數據采集、MATLAB虛擬仿真、數據分析，詳細過程如圖2所示。

圖2 技術設計流程圖

3 高程擬合

由于我國使用的是正常高系統，而GPS高程大部分是在WGS-84坐標下的大地高，無法正常使用，所以要對GPS測得數據進行擬合[6]，高程擬合中各個曲面的對應關系如圖3所示。

圖3 高程擬合示意圖

似大地水準面沿鉛垂線方向至參考橢球面的距離稱為高程異常[7]，一般用δ來表示。

式中：H為大地高，Hr為正常高，δ為高程異常。

目前高程擬合常用的方法有曲面擬合、多面函數法、樣條插值法等，在本研究中選用二次曲面擬合，曲面擬合回歸分析中一種最常用的方法，而且在實際測量數據處理中廣泛應用。對于高程擬合來說，一般采取一次、二次、三次曲面擬合，曲面擬合次數的選取要以實際問題為主，具體問題具體分析[8]。在高程變化大的區域，擬合次數過小則難以表達區域高程的變化；在高程變化小的區域，擬合次數過大則會導致最后數據冗余，而結果和擬合次數小的差別不大[9]。

曲面擬合方程的一般形式為一次擬合、二次擬合、三次擬合。

本例中以二次擬合為例，誤差方程如下：

對于二次擬合要求6個系數，1個已知點可列1個方程，所以至少6個已知點，聯立（10）、（11）式，由系數矩陣B和常數項L，即可求解擬合系數。

4 實例驗證與精度分析

本例選用華北理工大學測區內靜態GPS測量數據，同時在這些點上聯測了二等水準。

本次實驗共測了17個GPS點，其中選取15個點進行高程擬合實驗，另外選取5個點進行最后擬合系數的驗證，說明IGGI、IGGⅢ擬合回歸方程的正確性。本例中，選擇曲面二次擬合，對測區進行高程擬合實驗。測區選點的原則：1）本研究選擇的是二次曲面擬合模型，該模型有6個未知點，一個已知點可建立一個誤差方程，所以至少需要6個聯測點的坐標才可求解參數。2）因為地球不是一個規則的曲面，地面高低起伏，不同地區的高程異常也不一樣，因此選點盡可能地要包圍整個測區，如圖4所示。靜態GPS和水準原始數據如表1所示。

圖4 測區點位選擇圖

表1 靜態GPS和水準原始數據

4.1 IGGⅠ方案和IGGⅢ的橫向對比

對于IGGI、IGGⅢ方案，都令K0=1.5、K1=2.5，比較IGGI、IGGⅢ方案對粗差的抗性程度。IGGI、IGGⅢ收斂速度對比如圖5所示。

從圖5中的折線變化趨勢可以看出，IGGI方案一共迭代了8次，從第7次后中誤差趨于穩定。IGGⅢ方案一共迭代了6次，從第4次后中誤差趨于穩定。可以看出，IGGⅢ方案比IGGI方案迭代次數更少、收斂速度更快，而且最后穩定的中誤差基本相同。表明IGGⅢ方案比IGGI方案在速度上更有優勢，在后面的數據對比中，選用IGGⅢ方案。

圖5 IGGI、IGGⅢ收斂速度對比圖

4.2 IGGⅢ的縱向對比

對于IGGⅢ來說，通過改變不同的K值，測試在不同K值下的抗差估計效果。K值對于IGGⅢ來說相當于閾值，閾值越小，對粗差的門限越小，粗差剔除也會相對越快，但是也有可能剔除掉偶然誤差。IGGⅢ中不同K值結果對比如表2所示。

表2 IGGⅢ中不同K值結果對比

從上表可以看出，隨著K0、K1取值的不同，即閾值從小到大，迭代次數依次增加，中誤差也逐漸減小。雖然閾值越小迭代速度變快，但是可能會剔除錯誤，導致最后擬合效果變差。所以要具體情況具體分析，選擇合適的閾值。

4.3 最小二乘和IGGⅢ方案的對比

驗證最小二乘和IGGⅢ方案對抗差的不同效果，選取5個控制點作為檢驗點，檢驗點數據如表3所示。

表3 檢查點數據

通過用擬合出二次項的系數，來計算新的高程異常值和原來的高程異常之間的差異。最小二乘和IGGⅢ結果對比如表4所示，而最小二乘擬合和IGGⅢ收斂速度對比如圖6所示。

圖6 最小二乘擬合和IGGⅢ收斂速度對比圖

表4 最小二乘和IGGⅢ結果對比

本研究中利用MATLAB 2016a對IGGⅢ擬合的平面和最小二乘擬合的平面在三維空間中進行顯示，并進行曲面對比，如圖7、圖8、圖9所示。即可以從曲面的趨勢走向來直觀展示兩者之間的差異。

圖7 最小二乘高程擬合曲面圖

圖8 IGGⅢ高程擬合曲面圖

圖9 最小二乘和IGGⅢ高程擬合對比圖

從圖7、圖8、圖9中可以看出，對于含有粗差的數據，用典型的最小二乘法擬合，計算時誤差平均分配，從而導致擬合出新的高程異常和原來的高程異常差值較大；而用IGGⅢ方案擬合的高程，對誤差進行了降權處理，擬合出的高程異常和原來的差值較小。

5 結論

筆者采用標準化殘差的方法構建IGGI和IGGⅢ權函數，第一種在IGGI和IGGⅢ中都取K0=1.5、K1=2.5，第二種在IGGⅢ權函數中取K0=1、K1=3，分別對GPS靜態測量和水準聯測的數據進行高程擬合。通過實驗分析可以看出，IGGI和IGGⅢ都對含有粗差的數據有較好的抗性，IGGⅢ比IGGI抗差效果更好；在IGGⅢ中對于K值的選取，應根據具體問題具體分析，對數據要求高時減小閾值，反之加大閾值。在日常的基礎測量工作中，粗差不可避免地會混入其中，而高程擬合又在實際測量工作中十分普遍，采用IGGⅢ方案可以很好地處理粗差。所以IGGⅢ不僅在高程擬合中，而且在所有包含粗差的測量數據處理中，都有很強的實用價值。