陳 毅，袁顯貴，楊一洋，彭 成

（1.東莞市萬江區(qū)測繪隊，廣東 東莞 523000；2.廣東方緯科技有限公司，廣東 廣州 510000；3.湖南省地質(zhì)礦產(chǎn)勘查開發(fā)局四〇三隊，湖南 常德 415000；4.杭州中房信息科技有限公司，浙江 杭州 310000）

通常，在地下管線竣工測量中直接獲取的都是平面直角坐標；而在提交測量成果時，驗收單位經(jīng)常要求測量單位將管線點的平面直角坐標和WGS84坐標同時標注在管線竣工測量圖上。目前常見的處理方法是先通過坐標轉(zhuǎn)換軟件將管線點的平面直角坐標轉(zhuǎn)換成WGS84坐標，再手工將管線點的平面直角坐標和WGS84坐標標注在竣工測量圖上。該方法不僅繁瑣，而且很容易出錯。為了解決這一問題，本文利用.NET與objectARX的AutoCAD二次開發(fā)技術(shù)，開發(fā)了坐標快速轉(zhuǎn)換標注程序，通過鼠標獲取管線點的平面直角坐標，再快速地將平面直角坐標轉(zhuǎn)換為WGS84坐標，然后按規(guī)范要求通過引線將兩種坐標標注在管線竣工測量圖上。

1 坐標實時轉(zhuǎn)換標注程序原理

由于WGS84坐標與我國已建的坐標系橢球基準不同[1]（具體參數(shù)見表1），所以地方坐標系平面直角坐標轉(zhuǎn)換到WGS84坐標，包括高斯正反算以及不同坐標系之間的坐標轉(zhuǎn)換。

表1 WGS84坐標系與我國已建坐標系的基準參數(shù)

1.1 高斯正算

高斯正算是指將大地坐標(B,L,H)轉(zhuǎn)換到平面直角坐標系(x,y,z)的過程[2]。高斯正算公式為：

式中，m= (L -L0)cosB ；為赤道到緯度B的子午線弧長[3-4]。

1.2 高斯反算

高斯反算是指將平面直角坐標系(x,y,z)轉(zhuǎn)換到大地坐標(B,L,H)的過程[5–6]。高斯反算公式為：

式中， L0為中央子午線經(jīng)度；Bf為對應于x的底點緯度，計算公式為：

本程序使用迭代法對其求解。

1.3 不同平面坐標系之間的坐標轉(zhuǎn)換

不同坐標系之間的坐標轉(zhuǎn)換，即利用公共點求解轉(zhuǎn)換參數(shù)的過程。常見方法包括三參數(shù)法、四參數(shù)法和七參數(shù)法。求解三參數(shù)至少需要一個公共點，四參數(shù)至少需要兩個公共點，七參數(shù)至少需要3個公共點。本文采用布爾莎模型求解七參數(shù)，并進行不同平面坐標系的轉(zhuǎn)換[7]。其具體公式為：

式中，ΔX0，ΔY0，ΔZ0為坐標平移量；R(ω)為3個方向的旋轉(zhuǎn)矩陣；m為縮放因子。

1.4 程序流程圖與代碼

程序流程如圖1所示。程序運行后，會彈出如圖2所示的窗口，將事先采集好的3個及以上公共點添加到griddataview表格中，點擊計算，得到七參數(shù)；然后點擊標注，將回到圖形工作空間，鼠標進入捕捉狀態(tài)，當鼠標拾取到管線點時，則以管線點為基點作一條引線，將管線點的平面直角坐標和通過七參數(shù)轉(zhuǎn)換得到的WGS84坐標同時標注出來。

圖1 程序流程圖

程序的核心代碼為：

decimal y1= y-500000.0M；

if (this.IsBigNumber== true)

{

y1= y1-(this.L0 / (int)this.Strip) * 1000000M；

}

y= y1；

decimal l0= this.L0 * 3600； ‘ 中央子午線轉(zhuǎn)為秒值 如=105×3 600

decimal e4= pow(e2,2)；

decimal e6= pow(e2,3)；

decimal e8= pow(e2,4)；

decimal e10= pow(e2,5)；

decimal e_12= pow(e2,6)；

decimal A1= 1 + 3 * e2 / 4 + 45 * e4 / 64 + 175 * e6 /256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 /1048576；

decimal B1= 3 * e2 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 +2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288；

decimal C1= 15 * e4 / 256 + 105 * e6 / 1024 + 2205 *e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608；

decimal D1= 35 * e6 / 3072 + 105 * e8 / 4096 +10395 * e10 / 262144 + 55055 * e_12 / 1048576；

decimal E1= 315 * e8 / 131072 + 3465 * e10 /524288 + 99099 * e_12 / 8388608；

decimal F1= 693 * e10 / 1310720 + 9009 * e_12 / 5242880；

decimal G1= 1001 * e_12 / 8388608； ‘ 求底點緯度值Bf

decimal B0= x / (a * (1-e2) * A1)；

decimal Bf= 0.0M；

decimal FB= 0.0M；

decimal FB1= 0.0M；

decimal a0= a * (1-e2)；

decimal delta= Math.Abs(Bf-B0)；

while (delta ＞ 4.8E-11M) ‘ 0.000000000048M

{

Bf= B0；

FB= a0 * (A1 * Bf-B1 * sin(2 * Bf) + C1 * sin(4 * Bf) -D1* sin(6 * Bf) + E1 * sin(8 * Bf)-F1 * sin(10 * Bf) + G1 * sin(12 * Bf))；

FB1= a0 * (A1-2 * B1 * cos(2 * Bf) + 4 * C1 *cos(4 * Bf)-6 * D1 * cos(6 * Bf) + 8 * E1 * cos(8 * Bf)-10 *F1 * cos(10 * Bf) + 12 * G1 * cos(12 * Bf))；

B0= Bf + (x-FB) / FB1；

delta= Math.Abs(Bf-B0)；

}

decimal sinBf= sin(Bf)；

decimal sinBf2= sinBf * sinBf；

decimal W= sqrt(1-e2 * sinBf2)；

decimal W3= W * W * W；

decimal N= a / W；

decimal t= tan(Bf)；

decimal t2= t * t；

decimal t4= t2 * t2；

decimal cosBf= cos(Bf)；

decimal cosBf2= cosBf * cosBf；

decimal yy= e12 * cosBf2； //η2

decimal Mf= a0 / W3；

decimal y_N= y / N；

decimal y_N2= y_N * y_N；decimal y_N4= y_N2 * y_N2；‘計算經(jīng)偉度值B、L

decimal t_B= Bf*p-(p * t / (2 * Mf) * y * y_N) * (1-(5 + 3* t2 + yy-9 * yy * t2) * y_N2 + (61 + 90 * t2 + 45 * t4) * y_N4 / 360)；

decimal t_L= (p / cosBf) * y_N * (1-(1 + 2 * t2 + yy) * y_N2 / 6 + (5 + 28 * t2 + 24 * t4 + 6 * yy + 8 * yy * t2) * y_N4 / 120)；

L= t_L + l0；

B= t_B / 3600； ‘ 轉(zhuǎn)為度

L= L / 3600； ‘ 轉(zhuǎn)為度

圖2 程序運行界面

2 程序應用實例

本文以東莞市萬江區(qū)測繪隊承接的110kV莆鳳莆泰線架空改造線路竣工測量工程為例，該工程線路全長9.6 km，采用全站儀與GPS對管線進行竣工測量，管線點共計948個，施測坐標為珠區(qū)坐標系。運行本程序后，彈出七參數(shù)解算窗口，輸入3個公共點坐標，將自動解算七參數(shù)；再點擊標注，便能按照規(guī)范要求自動將管線點的珠區(qū)坐標和WGS84坐標標注出來。標注效果如圖3所示。

圖3 程序標注效果圖

3 結(jié) 語

本文通過自主研發(fā)的大地坐標快速標注程序，實現(xiàn)了地方坐標系向WGS84坐標系的快速轉(zhuǎn)換，并按照規(guī)范要求自動將地方坐標與大地坐標標注在管線竣工測量圖上。與傳統(tǒng)作業(yè)方法相比，該方法節(jié)省了80%以上的作業(yè)時間，工作效率得到了大幅提升，且坐標轉(zhuǎn)換精度與傳統(tǒng)作業(yè)方法基本一致，具有非常高的應用價值。