重庆市都市区高程异常模型的建立及RTK测量坐标转换系统开发

何德平，王昌翰

(重庆市勘测院，重庆 400020)

1 引言

如何实现WGS－84坐标到城市独立坐标系之间的高精度转换是进行RTK测量的关键，在高斯投影后平面采用二维四参数转换，难点是高程转换。当区域不大时，高程采用比较简单的平面拟合法。大范围的WGS－84高与正常高的转换关系通常是利用重力数据、水准测量数据、GPS测量数据经过一定的算法建立高精度的似大地水准面。

重庆市都市区面积约6 000 km2，平均海拔高程为320 m，最高海拔高为 950 m，属丘陵区、山地地形。实践表明，一般商业软件采用的平面高程拟合算法不适合重庆山地地形特点和大范围RTK测量。本文总结了在已有城市控制点资料，包括平面数据、高程数据、GPS数据的基础上，考虑山地地形起伏大的特点，经过对多种高程拟合法对比分析，选择合适算法，建立满足重庆市都市区城市测量与工程建设需要的高精度高程异常模型，并开发RTK测量坐标转换系统。

2 主要研究内容及总体思路

主要研究内容:

(1)重庆市都市区已有城市测量控制点资料，包括平面数据、高程数据、GPS数据兼容性分析;

(2)地形改正研究;

(3)GPS高程拟合算法对比研究;

(4)基于Windows及WindowsCE RTK测量坐标转换系统软件开发。

总体思路是首先对都市区范围内各个时期所测量的已有资料进行兼容性分析，经筛选后，提出满足一定密度及残差要求的坐标转换控制点。经高斯投影后，平面坐标采用四参数转换法，转换参数采用C级GPS点计算。考虑到重庆是山地地形，精确的似大地水准面求取是难点，因此对平面法、二次曲面拟合法、基于地形改正和二次曲面拟合的移去－恢复算法进行对比研究，选用适合重庆特点的高程拟合方法。实现的难点是基于地形改正和二次曲面拟合的移去－恢复算法。

为便于地形改正计算，将重庆市都市区范围分成8 km×8 km的格网，各块分别计算地形改正并进行二次曲面拟合。由于控制点分布及测量误差的影响，各区块接边处转换结果存在一定的差异，因此在实际计算时各分块之间均重叠 3 km～5 km的范围，重叠部分的高程异常取中数。

根据各分块数据，按200 m×200 m格网间距计算都市区高程异常改正模型，按DEM GRD格式存储。由该高程异常模型经双线性内插而求得某一点的高程异常值。

3 相关算法

3.1 地形改正［6］

平面高程拟合法在平坦地区精度较高，在丘陵特别是山地的拟合结果与直接高程测量结果相差较大，说明大地水准面高程异常与地形相关。地形改正按式(1)计算。

式中:G为引力常数;ρ为地球质量密度;h为是流动单元的平均高程;hp为计算点高程;γ为参考椭球面上的正常重力;l为计算点到流动单元的距离;dxdy为流动单元的面积。

3.2 “移去－拟合－恢复”法原理［5］

根据物理大地测量学的理论，大地水准面差距N包含3个分量:

式中，Ngm是长波项，Ns是短波项，δN是由于地表的起伏引起的剩余分量，即地形改正项。在局部地区，引起大地水准面差距N无规律变化的主要原因是地表起伏，若在N中去掉地形改正项δN，则对于局部区域而言，应该是一个近似规则的函数。在已知的GPS水准点上，首先去掉δN，再通过高程拟合得到待定点高程，然后再加上该点处的地形改正σni即得到该点的正常高。

3.3 高程拟合算法

(1)平面拟合法

在小区域且较为平坦的范围内，通常采用平面逼近似大地水准面。某点的高程异常ξ与该点的平面坐标(x，y)有如下关系式:

式中:a1、a2、a3为模型参数。

若公共点数目大于3，则相应的误差方程为:

写成矩阵形式为:

(2)二次曲面拟合法

公共点上的高程异常与平面坐标之间，存在:

式中:a0、a1、a2、a3、a4、a5为待定参数。因此，区域内至少需要6个公共点。当公共点多于6个时，则可列出相应的误差方程。

用矩阵表示:

以上两种算法在列出误差方程后，依据最小二乘原理VTPV=最小，求向量A的解，进而可求出向量ξ。

4 RTK测量坐标转换系统开发

RTK测量坐标转换系统分为运行于Windows和Windows CE两种版本。Windows CE版本运行于GPS的RTK手簿上，可实现野外实时转换并将转换成果通过GPRS传回服务器，便于统一管理。

4.1 数据准备

在高斯投影后平面转换采用二维四参数法，利用240点C级GPS控制点作为公共点。在对现存各个时期所测量的四等及以上高程控制点数据进行高程兼容性分析的基础上，剔出误差大的点，最后选择950个控制点(其平面坐标为一级导线精度)，覆盖重庆市都市区面积约6 000 km2范围，高程点平均密度为 3 km～4 km。

4.2 地形改正的计算步骤［5］

计算过程:

(1)以8 km×8 km分块，并根据各分块周边控制点分布情况，外扩 3 km～4 km为计算区域。

(2)以计算点P周边4 km×4 km范围分成小矩形块 100 m×100 m作为流动单元，如图1所示。

图1 地形改正计算示意图

(3)利用该区域DEM计算流动单元的平均高程h，即h=(h1+h2+h3+h4)/4，hi表示流动单元4个角点的高程。

(4)按照式(1)计算P点的地形改正。共中各参数如下:

①采用都市区1∶5万DEM(格网间距25 m)

②G=6.67259e－11;ρ=2.67e3

③γ=978032.68［1+0.0053024sin2(φ)

－0.0000058sin2(2φ)］mGal

地形改正计算量大，为了避免实时计算地形改正速度较慢的问题，事先计算整个都市区范围的高程异常改正模型，并存储格式为格网间距为 200 m×200 m DEM GRD格式。任意某点高程异常值利用高程异常改正模型采用双线性内插法求得。

4.3 双线性内插［2］

由正方形4个角点高计算内部某点高采用双线性内插法。

式中:a、b、c、d 为正方形四个格网点，l是格网边长。

4.4 系统的主要功能有

(1)重庆市独立坐标系转WGS－84坐标。

(2)WGS－84坐标转独立坐标系。高程转换提供平面转换、二次曲面转换、基于地形改正和二次曲面拟合的移去恢复算法三种转换方式。

(3)高程异常改正模型加密。

(4)控制点检校。系统具有对已知控制点检校测量功能，导入已知控制成果，得出检校精度。

(5)GPS RTK手簿与服务器通过GPRS保持通信，将实时转换成果传回服务器，以便于单位对控制点数据及时统一管理。

(6)按GPS RTK测量规范要求，提供图根导线点、碎部点、一级导线点、二级导线点的成果报表格式。

5 RTK测量坐标转换系统运行效果及数据测试

5.1 系统界面

系统界面如图2所示。

图2 系统界面

5.2 计算结果

(1)地形改正

计算结果表明，地形改正随地形起伏变化比较明显，与重庆平均高程面 300 m高相差不大的地区，改正数在1 cm以下，500 m以上的区域，改正数在 3 km～6 km之间。重庆市都市区地形改正计算结果如图3所示。

图3 重庆都市区地形改正晕渲图

(2)重庆市都市区高程异常模型计算

基于地形改正和二次曲面拟合的移去恢复算法的重庆市都市区高程异常模型如图4所示。

图4 重庆市都市区高程异常模型晕渲图

5.3 三种坐标转换方法精度比较

为了检验WGS－84坐标到重庆市独立坐标系平面及高程的转换精度，在都市区布设均匀分布的12个控制点，通过将转换结果与实际测量结果比较，评定数据转换的精度。

对12个控制点GPS观测 8 h，分别解算WGS－84坐标和重庆市独立坐标系坐标，高程按二等水准精度接测。

高程转换利用系统提供的平面拟合、二次曲面拟合、基于地形改正和二次曲面拟合的移去－恢复算法三种方法分别计算正常高，转换坐标与实际测量坐标的平面位置较差、高程较差及中误差如表1所示。

WGS－84系坐标到重庆市独立坐标系转换精度统计表 表1

结果表明，三种高程转换方法中，基于地形改正和二次曲面拟合的移去恢复算法更适合重庆地形特点，精度更高。

6 结论

(1)二次曲面高程拟合比平面拟合法更能适合重庆山地城市的地形特点，拟合精度更高。

(2)基于地形改正和二次曲面拟合的移去恢复算法精度高，适合重庆地形特点，特别是离平均高程面大的地区高程拟合精度明显提高。

(3)利用重庆市都市区已有的950个高程较正点，采用基于地形改正和二次曲面拟合的移去恢复算法得到的都市区高程异常模型高程精度总体在 5 cm以下，大部分区域都在 3 cm以下。

(4)开发的RTK测量坐标转换系统能将WGS－84坐标转换为重庆市独立坐标系坐标，具有自动化程度高，成果规范的特点，平面及高程精度能满足城市测量和工程建设的要求。

［1］李建成等.地球重力场逼近理论与中国2000似大地水准面的确定［M］.武汉:武汉测绘科技大学出版社，2003

［2］李志林，朱庆.数字高程模型［M］.武汉:武汉测绘科技大学出版社，2000

［3］刘大杰，施一民，过静珺.全球定位系统(GPS)的原理与数据处理［M］.上海:同济大学出版社，2003

［4］CJJ 8－99.城市测量规范［S］.

［5］王爱生，欧吉坤，赵长胜.“移去－拟合－恢复”算法进行高程转换和地形改正计算公式探讨［J］.测绘通报，2005(4)

［6］宁津生，罗自才，杨洁吉等.深圳市1 km高分辨率厘米级高精度大地水准面的测定［J］.测绘学报，2003，32(2)

［7］赵建虎，刘经南，张红梅.顾及非格网数据考虑地形改正的GPS水准高程拟合［J］.武汉测绘科技大学学报，1999，24(4)

［8］赵景权，周景龙.GPS高程转换中考虑地形改正的几何算法［J］.东北测绘，2001(1)