顧及起算點誤差的區域參考框架約束方法研究

王盼龍 侯汶材 蔣光偉 王 斌 程傳錄 李 康

1 自然資源部大地測量數據處理中心,西安市友誼東路334號,710054 2 中建新疆建工(集團)有限公司西北分公司,西安市科技西路2528號,710086

區域參考框架是全球參考框架的加密和延伸,選擇合適的平差基準對區域參考框架的建立和維持至關重要。利用GNSS技術建立和維持區域參考框架時,通常選擇區域內部及周邊均勻分布的IGS站作為起算點,并利用國際地球自轉服務(International Earth Rotation Service, IERS)組織發布的國際地球參考框架(international terrestrial reference frame, ITRF)的坐標和線性速度成果,計算指定歷元下的IGS站坐標作為先驗坐標。然而,隨著IGS站的長期運行,受地殼變化、地震和更換天線等因素影響,該方法計算的ITRF坐標常包含粗差[1]。因此,使用合適的基準約束消除或削弱起算點的粗差對保證框架維持成果的可靠性至關重要。

GNSS自由網是秩虧的,為獲取可靠的框架成果,需要附加合適的基準約束條件轉換到指定的參考框架[2]。常用的GNSS網平差方法有參數加權法和最小約束法。參數加權法是給予起算點坐標指定的約束,將整網結果強制約束到由起算點構成的坐標框架,是我國工程領域應用最普遍的方法[3],大多數測繪工程均采用該方法實現全球框架到區域框架的轉換工作。然而,該方法需要給予起算點合適的坐標精度,實踐中常憑經驗給予坐標約束。當起算站點存在粗差時,會造成網形扭曲,嚴重影響待解點的定位結果[3-4]。最小約束法是利用相似變換方法構建起算點網形與自由網網形最優符合的約束條件[4],該方法附加的條件數為秩虧數,能有效保留控制網本身的基準信息。但當起算點坐標存在粗差時,仍然會造成整網存在系統誤差。采用IGS官網發布的站點瞬時框架坐標,對其進行參數約束平差是最常用且能有效消除粗差影響的框架維持方法[5]。然而,該方法依賴于IGS官網發布的瞬時框架坐標成果,會大大影響框架維持的時效性。

在已有研究基礎上,為解決起算點粗差對框架成果的影響,本文提出基于Helmert轉換模型的基準約束方法。首先對自由網解進行最小約束平差,然后利用Helmert轉換法識別并剔除存在粗差的起算點,最后利用最優起算點的先驗坐標構建相似變換約束條件,平差得到GNSS框架網的坐標成果。同時設計對比實驗,采用IGS站ITRF2014坐標分別構建參數加權約束、最小約束、顧及IGS站先驗坐標誤差的Helmert模型約束解算中國大陸構造環境監測網絡(下文簡稱陸態網)框架成果,分析成果的精度和一致性,以驗證本文所提算法的性能。

1 區域參考框架的基準約束方法

1.1 參數加權約束法

現代大地測量中不存在完全未知的參數,在確定GNSS網的起算基準時,通常采用全球IGS站的先驗坐標和精度,平差數學模型如下[2]:

(1)

式中,L為觀測值,E(L)為L的期望值,D(L)為L的方差,μX為參數X的先驗值,ΣX為參數X的先驗方差-協方差陣,QX為參數X的先驗單位協方差矩陣,PX為參數X的先驗單位權陣,ΣXL為X與L的協方差且為正定矩陣。取X的近似坐標為X0=μX,列出誤差方程:

(2)

(3)

(4)

參數加權法能解決GNSS自由網平差的秩虧問題,可獲得網的絕對定位結果。在實際使用中,究竟要給予何種程度的約束,在很大程度上需要憑經驗確定,尤其是當起算點坐標中存在粗差時,會導致整個網形的扭曲,嚴重影響框架成果的可靠性[6]。

1.2 最小約束法

最小約束的含義是在處理觀測量得到點位坐標時,僅包含定義其自身參考框架的必需數據或約束[6]。因此,需要在式(2)上附加秩虧數d(d=μ-t)個約束條件,其中,μ為待解參數個數,t為必要觀測個數。通常在GNSS網中,確定兩個參考框架的轉換關系是利用7個參數構成的相似轉換模型進行描述[7],即

X2=X1+Dθ

(5)

(6)

已知兩個坐標系下3個以上的公共點,利用最小二乘法可解得θ:

θ=(DTPXD)-1DTPX(X2-X1)

(7)

設B=(DTPXD)-1DTPX,式(7)可轉化為:

θ=B(X2-X1)

(8)

為消除式(2)的秩虧問題,附加約束如下:

(9)

(10)

因此,附加最小約束的法方程為:

(11)

式中,x0為轉換點的坐標差。合并式(10)和(11),得:

(12)

由式(12)可知,相似變換約束實質上是將待定點的控制網轉換到已知點所在的框架下。因為最小約束法附加的條件數等于法方程的秩虧數,因此不干擾控制網本身的基準信息。該方法可為秩虧自由網提供充分必要的起算條件,平差結果不改變由觀測信息所確定的內在網形,可最大限度地保留觀測精度[8]。

2 Helmert轉換基準約束方法

由于受地殼變化、地震、更換天線等因素影響,基準點可能存在粗差,這種粗差在解算前可能并不清楚。起算點的坐標粗差會扭曲網平差結果,引起觀測網的變形[9]。本文提出的Helmert轉換法可以有效剔除存在粗差的起算點,并實現平差后網形與自由網網形的最佳符合。Helmert轉換法的具體解算流程如圖1所示。其中,Helmert相似變換模型為[7]:

(13)

式中,TX、TY、TZ為3個平移參數,RX、RY、RZ為3個旋轉參數,TS為尺度參數。利用IGS站先驗坐標構建的最小約束條件,平差得到整網結果,即(X0,Y0,Z0)、(X1,Y1,Z1)均為已知量,建立7個轉換參數的誤差方程:

(14)

(15)

3 實驗與分析

實驗選取中國境內和周邊15個IGS站和249個陸態網GNSS基準站,觀測時間為2020-08-17～23,采用高精度GNSS數據處理軟件GAMIT/GLOBK 10.70分區解算獲取單天的自由網解,然后合并子區獲取整個區域的單天自由網解,最終合并單天解獲取周自由網解SINEX文件。SINEX文件中IGS基準站的先驗坐標由IERS發布的ITRF2014框架的坐標和線性速度成果歸算獲取。利用GAMIT進行基線解算時,主要的控制參數設置見表1,基線解算結果見表2。

表1 主要控制參數

表2 基線解算精度

區域參考框架的可靠性至關重要,將其與ITRF參考框架聯系起來是實現區域參考框架最有效的手段。本文采用不同的基準約束方案解算陸態網在ITRF2014框架下的坐標,利用Python編程實現3種平差方法,并構建4種基準約束條件進行整網平差。

方案1:以15個IGS站為起算點,N、E、U方向上的坐標先驗精度分別為1 mm、1 mm、2 mm。

方案2:以15個IGS站構建七參數相似變換,7個參數先驗精度分別為10-5m、10-5m、10-5m、10-3ppb、10-6″、10-6″、10-6″。

方案3:采用本文提出的Helmert轉換法,將15個IGS站作為起算點構建相似變換,解算整網成果。利用Helmert坐標轉換模型,計算IGS站的成果與先驗坐標的偏差值,將誤差超限最大的IGS站作為待解點,利用剩余IGS站重新構建七參數相似變換,迭代計算,直至不存在誤差超限的IGS站。其中,7個參數的先驗精度分別為10-5m、10-5m、10-5m、10-3ppb、10-6″、10-6″、10-6″。

方案4:以方案3中滿足限差的IGS站為起算點,N、E、U方向上的坐標先驗精度分別為1 mm、1 mm、2 mm。

為了驗證解算結果的精度和可靠性,下載IGS官網發布的周解igs20P2119.ssc文件,提取中國境內及周邊15個IGS站的坐標及解算精度作為起算基準進行整網平差。將解算的成果作為參考值,對比不同方案解算成果的精度和可靠性:

(16)

分別采用4種基準約束平差方案獲取陸態網的ITRF2014框架成果,計算4種方案結果與參考值的誤差值(誤差的絕對值),圖2為平面方向誤差值分布結果,圖3為大地高方向誤差值分布結果。

圖2 平面解算誤差分布Fig.2 Distribution of plane error

圖3 大地高解算誤差分布Fig.3 Distribution of ellipsoidal error

由圖2可知,方案1和2平面坐標誤差較大,方案3和4平面誤差較小。方案1設定所有IGS站1 mm的平面坐標精度約束,但CHAN、DAEJ、STK2和URUM站平面的先驗坐標與周解的誤差超過3 cm,導致網形扭曲,且距離誤差點越近,站點的誤差值越大。方案2采用最小約束法,保證了解算前后網形的一致性,但起算點粗差同樣會代入整網,影響平差結果。方案3解算誤差分布均勻,且均優于1 cm,說明該方法可有效削弱起算點粗差對解算結果的影響。方案4與方案3解算結果相近,說明利用方案3選取的起算點的平面坐標不再包含粗差,滿足1 mm先驗坐標精度。

由圖3可知,方案1和4大地高誤差較大,方案2和3大地高誤差較小。方案1大地高誤差分布不均勻,IGS站大地高誤差較大,陸態網站大地高誤差較小。方案2大地高誤差分布均勻,IGS站和陸態網站大地高誤差均較小,說明最小約束法能有效削弱大地高方向的誤差。方案3大地高誤差分布均勻且最小,說明Helmert轉換法能削弱粗差對平差結果的影響,是最優的區域框架基準約束方法。方案4中部分IGS站大地高誤差較大,且中國東部的陸態網站大地高誤差較大。

為進一步分析不同方案定位結果的精度和可靠性,分別統計4種基準約束平差方案的站點誤差,表3(單位cm)和表4(單位cm)分別為15個IGS站和249個陸態網站點坐標的誤差統計結果(最大值為絕對值的最大值)。

表3 IGS基準站坐標誤差

表4 陸態網基準站坐標誤差

由表3和表4可知,方案3中IGS站和陸態網站的解算誤差均優于1 cm,各方向RMS值均優于5 mm,是4種方案中的最優方案。方案1中IGS站的誤差較大,原因為部分站點存在未知的粗差,因此給定的先驗精度不合理會造成網形扭曲。方案2相對于方案1大地高方向誤差較小,說明最小約束法能在一定程度上降低網形扭曲造成的誤差。方案3利用Helmert轉換法可有效識別并剔除存在粗差的站點,同時利用最小約束法確保解算前后網形的一致性,實現整網平差成果的高精度和高可靠性。方案4相對于方案3大地高方向誤差較大,其原因為剔除的IGS站大地高方向限差為3 cm,當選定的IGS站采用2 mm精度約束時,高程方向的粗差同樣會造成網形扭曲。

由表2可知,基線解算結果滿足mm級精度,保持網形的情況下引入合適的基準約束,即可獲取高精度的框架成果。因此,解算前后單位權方差的一致性與網形變化強相關。為分析不同解算方案對網形的影響,分別對解算前后單位權方差進行卡方檢驗(χ2)。卡方檢驗的相關系數可表征兩個變量之間的相關程度,相關系數絕對值越大,相關性越強;反之,相關系數越接近于0,相關度越弱。相關強度可以通過表5進行判斷[10]。

表5 相關系數

由表6可知,方案3的χ2值最小,說明解算前后單位權方差變化小,即網形一致性最好,實現了解算前后站點間相對關系不變。方案1和方案4的χ2值超過0.4,說明解算前后單位權方差變化大,給定的基準約束造成網形變化較大。采用方案3解算時,利用Helmert模型計算出未超限的11個IGS站的先驗坐標到平差結果的轉換參數,結果見表7。

表6 χ2檢驗結果

表7 Helmert轉換法的轉換參數

當起算站點先驗坐標存在未知誤差時,采用方案3的平差算法可以有效識別先驗坐標存在較大粗差的站點。同時,最小約束平差方法可保證整網間站點相對關系的可靠性,利用轉換參數吸收起算點的小粗差,獲取高精度和高可靠性的區域框架成果。

4 結 語

為解決起算點IGS站的坐標粗差對區域框架維持的影響,本文提出基于Helmert轉換法的基準約束方法。該算法基于高精度的GNSS自由網成果,利用Helmert轉換模型保持解算前后網形的一致性,然后迭代計算起算點的坐標誤差,實現粗差點的識別與剔除。實驗結果表明,本文提出的Helmert轉換法與傳統的參數加權法、最小約束法相比,定位精度優于傳統策略,且能有效削弱起算點坐標粗差對網形結構的影響,提高框架成果的可靠性。

致謝：本文使用的GNSS數據為中國地震局提供的陸態網基準站觀測數據,數據處理采用GAMIT/GLOBK軟件,在此一并表示感謝。