基于抗差估計(jì)的變形監(jiān)測基準(zhǔn)網(wǎng)坐標(biāo)系統(tǒng)轉(zhuǎn)換

呂偉才， 張瑞華， 張則友， 周學(xué)年， 陳 芳

（1.安徽理工大學(xué) 測繪學(xué)院，安徽 淮南 232001；2.淮南礦業(yè)集團(tuán) 李嘴孜煤礦，安徽 淮南 232074）

0 引言

李嘴孜煤礦由原李嘴孜煤礦和孔集煤礦合并而成。2008年，李嘴孜煤礦孔井工業(yè)廣場報(bào)廢，對(duì)主、副井井筒－42.8m以上進(jìn)行了充填，且主井在－250m水平、副井在－140m水平設(shè)置了隔斷平臺(tái)。按該礦生產(chǎn)計(jì)劃，將逐步回采原孔井工業(yè)廣場保護(hù)煤柱。為監(jiān)測井下開采對(duì)主、副井井筒的影響，掌握地表及井筒變形形態(tài)，需在工業(yè)廣場及井筒附近地表布設(shè)工業(yè)廣場監(jiān)測系統(tǒng)。

工業(yè)廣場平面監(jiān)測系統(tǒng)（高程監(jiān)測系統(tǒng)采用幾何水準(zhǔn)方法建立）由10個(gè)點(diǎn)組成的基準(zhǔn)網(wǎng)和由54個(gè)監(jiān)測點(diǎn)組成的監(jiān)測網(wǎng)構(gòu)成，基準(zhǔn)網(wǎng)和監(jiān)測網(wǎng)的基準(zhǔn)時(shí)段測量，采用GPS靜態(tài)相對(duì)定位進(jìn)行。為分析工業(yè)廣場井筒移動(dòng)變形規(guī)律與井下開采的關(guān)系，必須進(jìn)行GPS坐標(biāo)系和礦區(qū)坐標(biāo)系之間的坐標(biāo)系統(tǒng)轉(zhuǎn)換。本文主要討論GPS基準(zhǔn)網(wǎng)在平面上進(jìn)行坐標(biāo)系統(tǒng)轉(zhuǎn)換的數(shù)學(xué)模型及轉(zhuǎn)換后GPS網(wǎng)的質(zhì)量評(píng)價(jià)方法，其他內(nèi)容將另文介紹。

1 GPS基準(zhǔn)網(wǎng)平面坐標(biāo)系統(tǒng)轉(zhuǎn)換模型

1.1 GPS基準(zhǔn)網(wǎng)概況

變形監(jiān)測系統(tǒng)的平面基準(zhǔn)網(wǎng)由10個(gè)GPS點(diǎn)組成，如圖1所示，其中TZHN為三等控制點(diǎn)，PWZ5、PDK5、PLZZ 為 四 等 控 制 點(diǎn)，XZHZ、CPYDM、CAUST為安徽理工大學(xué)建立的C級(jí)GPS點(diǎn)。這些點(diǎn)均具有礦區(qū)坐標(biāo)系統(tǒng)（BJ－54坐標(biāo)系）下的平面坐標(biāo)，為坐標(biāo)系統(tǒng)轉(zhuǎn)換提供了基礎(chǔ)條件。JZ01～JZ03為新建監(jiān)測系統(tǒng)的工作基準(zhǔn)網(wǎng)。

圖1 GPS基準(zhǔn)網(wǎng)示意圖

GPS平面基準(zhǔn)網(wǎng)采用10臺(tái)中海達(dá)雙頻接收機(jī)觀測2個(gè)時(shí)段，每時(shí)段長為2h?；€向量解算采用HDS2003隨機(jī)軟件和廣播星歷進(jìn)行，經(jīng)同步邊、重復(fù)邊、同步環(huán)、異步環(huán)檢核［1］，外業(yè)觀測質(zhì)量滿足D級(jí)GPS網(wǎng)的要求。隨后的數(shù)據(jù)處理，采用具有抗差功能的“GPS（監(jiān)測）網(wǎng)數(shù)據(jù)處理軟件包（GMDPS 2.5）”［2］進(jìn)行。

以基本位于測區(qū)中部的JZ01點(diǎn)為已知點(diǎn)，在WGS－84空間直角坐標(biāo)系下進(jìn)行GPS網(wǎng)的空間無約束平差和精度評(píng)定［3］，然后將平差坐標(biāo)及其協(xié)因數(shù)陣轉(zhuǎn)換為大地坐標(biāo)系中的成果，最后投影到高斯平面上［4］。該網(wǎng)投影到高斯平面上后，最弱點(diǎn)點(diǎn)位中誤差為±8.0mm，平均為±5.2mm。可見投影后，該GPS網(wǎng)的精度很高，這就為實(shí)施該GPS網(wǎng)的坐標(biāo)系統(tǒng)轉(zhuǎn)換提供了基礎(chǔ)。

1.2 轉(zhuǎn)換參數(shù)的最小二乘估計(jì)模型

設(shè)GPS平面基準(zhǔn)網(wǎng)點(diǎn)經(jīng)過空間無約束平差后，投影到WGS－84高斯平面直角坐標(biāo)系中的坐標(biāo)為（x，y）G，坐標(biāo)聯(lián)測點(diǎn)在BJ－54坐標(biāo)系中的高斯平面直角坐標(biāo)為（x，y）T，坐標(biāo)轉(zhuǎn)換聯(lián)測點(diǎn)個(gè)數(shù)為k，其中m（m＞2）個(gè)用于求解轉(zhuǎn)換參數(shù)，這m個(gè)坐標(biāo)聯(lián)測點(diǎn)稱為坐標(biāo)轉(zhuǎn)換基準(zhǔn)點(diǎn)。在平面上2個(gè)坐標(biāo)系之間的轉(zhuǎn)換模型［5］為：

其中，（x0，y0）為平移參數(shù)；λ為尺度因子；θ為旋轉(zhuǎn)角。在這m個(gè)轉(zhuǎn)換基準(zhǔn)點(diǎn)上，將（1）式寫成誤差方程形式為：

其中，i＝1，2，…，m；（vxi，vyi）表示（xi，yi）T的殘差；wxi＝（xG－xT）i，wyi＝（yG－yT）i為第i個(gè)基準(zhǔn)點(diǎn)的自由項(xiàng)。在這m個(gè)基準(zhǔn)點(diǎn)上，將（2）式寫成矩陣形式，有

其中，A為轉(zhuǎn)換系數(shù)矩陣；W為自由項(xiàng)為待求轉(zhuǎn)換參數(shù)；P為權(quán)陣；Q為m個(gè)點(diǎn)的互協(xié)因數(shù)陣（由GPS網(wǎng)的空間無約束平差的坐標(biāo)協(xié)因數(shù)陣投影到高斯平面上而得）。由最小二乘估計(jì)原理可得：

將求得的＾T代入（1）式可得任一GPS點(diǎn)在新BJ－54系坐標(biāo)系中的高斯平面直角坐標(biāo)（xi，yi）T。

1.3 轉(zhuǎn)換參數(shù)的抗差估計(jì)模型

淮南礦區(qū)的三角網(wǎng)是多年前建立的，由于受開采和其他因素的影響，圖1所示的礦區(qū)三角點(diǎn)（作為坐標(biāo)系統(tǒng)轉(zhuǎn)換的基準(zhǔn)點(diǎn)）中可能發(fā)生了移動(dòng)變形，為了控制某些基準(zhǔn)點(diǎn)中存在的變形對(duì)求轉(zhuǎn)換參數(shù)的影響，可采用抗差估計(jì)［5－7］。根據(jù)最小二乘估計(jì)原理，可得轉(zhuǎn)換參數(shù)初始平差結(jié)果：

進(jìn)行抗差估計(jì)時(shí)，采用的等價(jià)權(quán)函數(shù)為

其中，i，j＝1，2，…，2m；Dj＝｜vj｜／；k0為分位參數(shù)，k1為淘汰點(diǎn)，一般取k0＝1.5，k1＝2.5［3，6－7］；ˉpij為相關(guān)等價(jià)權(quán)元素，由其構(gòu)成相關(guān)等價(jià)權(quán)ˉP；vj為按（5）式計(jì)算的改正數(shù)；＾σ為按（6）式計(jì)算的單位權(quán)中誤差。

（7）式定義的等價(jià)權(quán)函數(shù)，根據(jù)Dj的大小，對(duì)基準(zhǔn)點(diǎn)采用不同的處理方法：當(dāng)Dj＜1.5時(shí)，認(rèn)為基準(zhǔn)點(diǎn)中不存在位移，按初始權(quán)參與轉(zhuǎn)換參數(shù)的求解；當(dāng)Dj在1.5～2.5范圍內(nèi)時(shí)，認(rèn)為基準(zhǔn)點(diǎn)中可能存在位移，對(duì)該基準(zhǔn)點(diǎn)采取降權(quán)處理；當(dāng)Dj＞2.5時(shí)，認(rèn)為該基準(zhǔn)點(diǎn)存在顯著的位移，此時(shí)該基準(zhǔn)點(diǎn)不參與轉(zhuǎn)換參數(shù)的解算。因而，（7）式定義的等價(jià)權(quán)函數(shù)對(duì)存在位移的基準(zhǔn)點(diǎn)具有抗差作用。

以等價(jià)權(quán)ˉP代替初始權(quán)P進(jìn)行迭代最小二乘估計(jì)，求得第k次迭代時(shí)的轉(zhuǎn)換參數(shù)＾T（k）、殘差V（k）及單位權(quán)中誤差＾σ（k）：

其中，tp為被剔除掉的基準(zhǔn)點(diǎn)分量個(gè)數(shù)。

由 （8）式、（9）式按 （7）式計(jì)算第k＋1次迭代時(shí)的等價(jià)權(quán)ˉP（k＋1）。

重復(fù)上述過程，直到前后2次迭代計(jì)算的轉(zhuǎn)換參數(shù)之差滿足預(yù)期的收斂精度時(shí)，停止迭代計(jì)算，此時(shí)即獲得了轉(zhuǎn)換參數(shù)的抗差估值。所求得的轉(zhuǎn)換參數(shù)是否顯著，還需要通過統(tǒng)計(jì)檢驗(yàn)的方法進(jìn)行確定。

2 GPS網(wǎng)坐標(biāo)系統(tǒng)轉(zhuǎn)換質(zhì)量評(píng)價(jià)

2.1 GPS網(wǎng)坐標(biāo)系統(tǒng)轉(zhuǎn)換后的質(zhì)量評(píng)價(jià)模型

轉(zhuǎn)換后GPS平面基準(zhǔn)網(wǎng)的質(zhì)量可以用點(diǎn)位中誤差、邊長中誤差、邊長相對(duì)中誤差和坐標(biāo)方位角中誤差等指標(biāo)來衡量，這些指標(biāo)均可由轉(zhuǎn)換后GPS點(diǎn)坐標(biāo)的協(xié)因數(shù)陣導(dǎo)出。

（1）GPS點(diǎn)轉(zhuǎn)換后坐標(biāo)的協(xié)因數(shù)陣。由坐標(biāo)轉(zhuǎn)換模型，可得任一GPS點(diǎn)（xi，yi）G轉(zhuǎn)換后的坐標(biāo)（xi，yi）T為：

其中，（xi，yi）T是 GPS 點(diǎn) 轉(zhuǎn) 換 后 的 坐 標(biāo)；（xi，yi）G是GPS點(diǎn)空間平差后投影到高斯平面上的坐標(biāo)；（x0，y0，λ，θ）為坐標(biāo)轉(zhuǎn)換參數(shù)。對(duì)（11）式全微分得：

由于尺度比λ和旋轉(zhuǎn)角θ都是很小的量，故（12）式可寫成：

由協(xié)因數(shù)陣傳播定律可得轉(zhuǎn)換后GPS點(diǎn)的協(xié)因數(shù)陣為：

其中，Q為抗差估計(jì)后轉(zhuǎn)換參數(shù)的協(xié)因數(shù)陣；QxG為GPS點(diǎn)投影到高斯平面上的坐標(biāo)協(xié)因數(shù)陣；矩陣A由各GPS點(diǎn)的投影坐標(biāo)按（15）式組成，即

A的維數(shù)為（2n，4），其中n為GPS點(diǎn)總數(shù)。由（14）式可得某一點(diǎn)的協(xié)因數(shù)陣或某2點(diǎn)的互協(xié)因數(shù)陣。

（2）GPS網(wǎng)轉(zhuǎn)換后的點(diǎn)位中誤差。雖然由（14）式可得任一GPS點(diǎn)轉(zhuǎn)換后的協(xié)因數(shù)陣，但由于GPS網(wǎng)在空間無約束平差及坐標(biāo)轉(zhuǎn)換過程中的單位權(quán)方差因子不同，因此應(yīng)將（14）式轉(zhuǎn)化為協(xié)方差陣形式，即

則可得該點(diǎn)轉(zhuǎn)換后的點(diǎn)位中誤差為：

（3）GPS網(wǎng)轉(zhuǎn)換后基線向量中誤差。GPS網(wǎng)轉(zhuǎn)換后的基線向量可用轉(zhuǎn)換后基線端點(diǎn)坐標(biāo)之差來表示，則任一基線可表示為：

其中，（x，y）為轉(zhuǎn)換后的GPS點(diǎn)的坐標(biāo)。由協(xié)方差傳播定律可得轉(zhuǎn)換后任一基線向量的方差為：

其中，a＝Δxij／Sij；b＝Δyij／Sij。則 GPS網(wǎng)轉(zhuǎn)換后任一基線向量的中誤差為：

相對(duì)中誤差為｜MSij｜／Sij。

（4）GPS網(wǎng)轉(zhuǎn)換后坐標(biāo)方位角中誤差。轉(zhuǎn)換后GPS網(wǎng)的任一基線向量的坐標(biāo)方位角為：

根據(jù)協(xié)方差傳播定律可得轉(zhuǎn)換后任一基線向量的坐標(biāo)方位角之方差為：

點(diǎn)位中誤差線、基線向量相對(duì)中誤差和坐標(biāo)方位角中誤差3項(xiàng)指標(biāo)，從不同角度反映了GPS平面基準(zhǔn)網(wǎng)轉(zhuǎn)換到BJ－54坐標(biāo)系下后的精度。根據(jù)這些指標(biāo)，可判定轉(zhuǎn)換后是否保留GPS網(wǎng)高精度的特點(diǎn)（點(diǎn)位誤差、基線向量中誤差）及GPS網(wǎng)是否發(fā)生了扭曲變形（坐標(biāo)方位角中誤差）。

2.2 GPS網(wǎng)坐標(biāo)系統(tǒng)轉(zhuǎn)換后的質(zhì)量

在圖1所示的GPS平面基準(zhǔn)網(wǎng)中，有7個(gè)坐標(biāo)聯(lián)測點(diǎn)。為保證轉(zhuǎn)換后GPS網(wǎng)的質(zhì)量，設(shè)計(jì)了4種轉(zhuǎn)換方案，通過分析比較，最后采用的方案為：TZHN、XZHZ、CAUST、CPYDM 作為坐標(biāo)系統(tǒng)轉(zhuǎn)換基準(zhǔn)點(diǎn)，PWZ5、PDK5、PLZZ等3點(diǎn)用于檢核轉(zhuǎn)換模型的精度。采用該方案進(jìn)行坐標(biāo)系統(tǒng)轉(zhuǎn)換，分別采用最小二乘估計(jì)（LS估計(jì)）和抗差估計(jì)（RS估計(jì)），GPS平面基準(zhǔn)網(wǎng)轉(zhuǎn)換到BJ－54坐標(biāo)系下后的點(diǎn)位中誤差、邊長相對(duì)中誤差和坐標(biāo)方位角中誤差的統(tǒng)計(jì)信息，見表1～表3所列。

采用最小二乘估計(jì)（LS估計(jì)）和抗差估計(jì)（RS估計(jì)）時(shí)，對(duì)于點(diǎn)位中誤差，最高精度分別為1.65、1.14cm，最低精度分別為4.36、1.71cm，平均精度分別為2.70、1.34cm；對(duì)于邊長相對(duì)中誤差，最高精度分別為2.33×10－6、1.41×10－6，最低精度分別為16.67×10－6、5.00×10－6，平均精度分別為7.69×10－6、3.57×10－6；對(duì)于坐標(biāo)方位角中誤差，最高精度分別為0.35"、0.06"，最低精度分別為7.63"、0.11"，平均精度分別為2.42"、0.09"。結(jié)合表1～表3可以看出，對(duì)于工業(yè)廣場平面監(jiān)測系統(tǒng)GPS平面基準(zhǔn)網(wǎng)的坐標(biāo)系統(tǒng)轉(zhuǎn)換，采用抗差估計(jì)（RS估計(jì)）的精度明顯優(yōu)于采用最小二乘估計(jì)（LS估計(jì)）的精度，更好地保留了GPS技術(shù)高精度的特點(diǎn)，為監(jiān)測井下開采對(duì)工業(yè)廣場主副井平面變形提供了可靠的基準(zhǔn)。

表1 點(diǎn)位中誤差統(tǒng)計(jì)表

表2 邊長相對(duì)中誤差統(tǒng)計(jì)表

表3 坐標(biāo)方位角中誤差統(tǒng)計(jì)表

3 結(jié)束語

礦區(qū)三角網(wǎng)，一般是在多年前甚至十幾年、幾十年前建立的，三角點(diǎn)難免受到開采的影響而產(chǎn)生變形，這就產(chǎn)生了如何選擇穩(wěn)定的三角點(diǎn)以保證轉(zhuǎn)換后GPS網(wǎng)的精度問題。本文以建立李嘴孜煤礦工業(yè)廣場監(jiān)測系統(tǒng)的GPS平面基準(zhǔn)網(wǎng)為例，采用抗差估計(jì)理論，實(shí)現(xiàn)平差系統(tǒng)自動(dòng)識(shí)別穩(wěn)定基準(zhǔn)點(diǎn)，建立了坐標(biāo)系統(tǒng)轉(zhuǎn)換參數(shù)求取的抗差估計(jì)模型和轉(zhuǎn)換后GPS網(wǎng)的質(zhì)量評(píng)價(jià)模型。目前，有的礦區(qū)已建立了單基站或網(wǎng)絡(luò)CORS系統(tǒng)，以解決礦山的日常測量和放樣問題，但CORS系統(tǒng)測量的精度，除受CORS系統(tǒng)本身的影響外，所采用的坐標(biāo)系統(tǒng)和高程系統(tǒng)轉(zhuǎn)換模型是否精確，是影響CORS系統(tǒng)測量精度的主要因素之一。采用本文的方法，也可以建立CORS系統(tǒng)測量的精確坐標(biāo)系統(tǒng)轉(zhuǎn)換模型。

目前，有的礦業(yè)集團(tuán)公司的常規(guī)控制網(wǎng)受開采影響破壞嚴(yán)重，致使采用GPS技術(shù)進(jìn)行礦山測量工作時(shí)，若建立地表移動(dòng)觀測站［8］，難以在測區(qū)周圍找到合適數(shù)量（6～8個(gè)）的穩(wěn)定三角點(diǎn)，或難以找到分布合理的三角點(diǎn)，使建立的GPS網(wǎng)坐標(biāo)系統(tǒng)轉(zhuǎn)換模型的適用性較差，或者轉(zhuǎn)換后GPS網(wǎng)的精度較低。因此，對(duì)于分布范圍廣泛的大型礦業(yè)公司，主管部門應(yīng)根據(jù)各礦井的開拓計(jì)劃，建立覆蓋本公司所有礦井范圍的GPS控制網(wǎng)，獲取坐標(biāo)系統(tǒng)轉(zhuǎn)換模型和高程基準(zhǔn)面，更好地服務(wù)于礦山的安全生產(chǎn)。

［1］安徽理工大學(xué)導(dǎo)航定位技術(shù)應(yīng)用研究所.孔井工業(yè)廣場地表移動(dòng)及主、副、風(fēng)井井筒沉陷觀測工程基準(zhǔn)網(wǎng)測量技術(shù)總結(jié)報(bào)告［R］.淮南：安徽理工大學(xué)，2010.

［2］余學(xué)祥，張華海，高井祥，等.GPS監(jiān)測網(wǎng)數(shù)據(jù)處理軟件包的研制［J］.礦山測量，2003（4）：11－14.

［3］余學(xué)祥，徐紹銓，呂偉才.GPS變形監(jiān)測數(shù)據(jù)處理自動(dòng)化：似單差法的理論與方法［M］.徐州：中國礦業(yè)大學(xué)出版社，2004：82－103.

［4］余學(xué)祥，呂偉才.空間直角坐標(biāo)的協(xié)因數(shù)陣轉(zhuǎn)換到高斯平面上的計(jì)算公式［J］.測繪信息與工程，1997（4）：18－21.

［5］Lü W C，Cheng S G，Yang H S，et al.The application of GPS technology to build a mine－subsidence observation station［J］.Journal of China University of Mining and Technology，2008，18（3）：377－380.

［6］高井祥，張華海，余學(xué)祥.礦區(qū)GPS網(wǎng)坐標(biāo)轉(zhuǎn)換的抗差模型［J］.中國礦業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，1999，28（2）：99－103.

［7］呂偉才，秦永洋，孫興平，等.顧橋礦區(qū)GPS基準(zhǔn)網(wǎng)坐標(biāo)轉(zhuǎn)換的高崩潰污染率抗差估計(jì)［J］.山東科技大學(xué)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，28（6）：9－15.

［8］呂偉才，秦永洋，孫興平，等.Kalman濾波在地表移動(dòng)觀測站沉降監(jiān)測中的應(yīng)用研究［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，34（9）：1370－1374.