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

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

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

1 坐標(biāo)實(shí)時(shí)轉(zhuǎn)換標(biāo)注程序原理

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

表1 WGS84坐標(biāo)系與我國(guó)已建坐標(biāo)系的基準(zhǔn)參數(shù)

1.1 高斯正算

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

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

1.2 高斯反算

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

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

本程序使用迭代法對(duì)其求解。

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

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

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

1.4 程序流程圖與代碼

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

圖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； ‘ 求底點(diǎn)緯度值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ì)算經(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 程序運(yùn)行界面

2 程序應(yīng)用實(shí)例

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

圖3 程序標(biāo)注效果圖

3 結(jié) 語(yǔ)

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