基于objectARX的地下管線坐標快速轉換標注程序研究

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

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

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

1 坐標實時轉換標注程序原理

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

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

1.1 高斯正算

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

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

1.2 高斯反算

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

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

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

1.3 不同平面坐標系之間的坐標轉換

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

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

1.4 程序流程圖與代碼

程序流程如圖1所示。程序運行后，會彈出如圖2所示的窗口，將事先采集好的3個及以上公共點添加到griddataview表格中，點擊計算，得到七參數；然后點擊標注，將回到圖形工作空間，鼠標進入捕捉狀態，當鼠標拾取到管線點時，則以管線點為基點作一條引線，將管線點的平面直角坐標和通過七參數轉換得到的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； ‘ 中央子午線轉為秒值 如=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；‘計算經偉度值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； ‘ 轉為度

L= L / 3600； ‘ 轉為度

圖2 程序運行界面

2 程序應用實例

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

圖3 程序標注效果圖

3 結 語

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