陈 毅，袁显贵，杨一洋，彭 成

（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%以上的作业时间，工作效率得到了大幅提升，且坐标转换精度与传统作业方法基本一致，具有非常高的应用价值。