一种基于DEM的CGCS2000城市平面坐标系统确定方法

高昭良

(1.武汉大学测绘遥感信息工程国家重点实验室，湖北武汉 430079;2.福州市勘测院，福建福州 350003)

1 引言

从2008年7月1日起正式启用中国大地坐标系统2000(CGCS 2000)作为国家法定的坐标系，作为我国新一代的平面基准［1］。该系统有利于GPS快速、精确地获取高精度城市坐标和高程成果，有利于GIS与GPS应用的结合，进一步提升数字城市的综合服务能力。

目前全国各地城市基本还是沿用原有的1954北京坐标系、1980西安坐标系统及地方坐标系统，因此，各城市都需要进行采用CGCS2000建立城市平面坐标系统。陈俊勇院士等探讨国家小比例尺地图的“西安80坐标系”更换为“大地参考系1980”(GRS80)，给出这些变化量的大小区间并讨论它对地图表示可能产生的影响［2］。李江卫、解斌等人以武汉市为例，提出了以投影改正平方和最小和投影变形大于 2.5 cm/km的点数最少为原则进行最佳中央子午线和抵偿高程面选取的方法［3］。武丰雷、刘曦灿等人提出建立斜轴墨卡托投影济南市坐标系［4］。

在讨论如何确定城市平面坐标系统选择时，相关文献基本采用两个方面来验证:一是采用有限的控制点作为计算结果验证，二是采用最小二乘法。由于城市区域范围较大，地形起伏分布不均匀，因此，有限的控制点无法全面准确概括。并且少量的地形突变会导致最小二乘法的最优评价结果畸变。本文提出的基于DEM的CGCS2000城市平面坐标系统确定方法可以兼顾解决这两方面问题。

2 城市平面坐标系统定义模型

一个城市坐标系统确定应根据以上基本原则，还得考虑城市具体地理范围和形状，最后来确定城市平面坐标系统的定义模型。根据《城市测量规范》，坐标系统应遵循以下原则:

(1)投影变形尽量小，投影长度变形值不大于2.5 cm/km(1/4 万)［5］。

(2)坐标系与国家(省)级控制点应联测，并建立严密的转换关系［6］。

2.1 高斯投影长度变形分析

由于定义国家大地坐标系的椭球面是一个凸起的、不可展平的曲面，当采用高斯正形投影将这个曲面上的元素投影到平面上时，投影后就会发生长度变形问题。高斯投影时，先把地面观测值归化至参考椭球面上，再把参考椭球面上的观测值归化至高斯平面上。两次归化产生的长度变形值称为投影变形。投影变形△D计算过程如下:

上式中各变量说明如下:

Hm:测量两端为高出参考椭球面的平均高程，单位:m;

H0:抵偿面高程，单位:m;

D:测量长度两端点平均高程面的水平距离，单位:m;

RA:归算边方向参考椭球法截弧的曲率半径，可取概值 6 371 000 m，单位:m;

ym:测距边两端点投影后近似横坐标平均值，单位:m;

Rm:测距边中点的平均曲率半径，可取概值6 371 000 m，单位:m。

2.2 采用格网占比法确定最佳关系

对于面积较大的城市或地区，在确定CGCS2000框架下的平面坐标系时，还应根据区域地理特点以及建设发展的实际情况，合理选取最佳中央子午线和抵偿高程面。不断调整高程参考面的大地高为H，使其投影变形满足要求，条件式可近似表示为:

若确定的中央子午线与高程抵偿面在城市坐标系统覆盖的范围内，能够较好地满足式(2)或式(3)的关系，则该中央子午线与高程抵偿面关系最佳。

3 基于DEM模型数据的计算模型

大部分学者计算最优值采用最小二乘法∑△D2=min［9，10］模型，由于选取的控制点有限，不能全面均匀代表该区域投影变形，因此，条件通过遍历搜索获得∑△D2=min的方案是否可以确定为最优，值得商榷。本文基于福州地区 1∶1万DEM数据，提取分布均匀的高程网格，利用以上的投影计算模型，计算最佳高程抵偿面。最佳高程抵偿面与中央子午线的计算采取遍历搜索寻优法:搜索范围在中央经度［Lmin，Lmax］，抵偿面高程［Hmin，Hmax］内，按适合精度的步长，以投影改正平方和最小为原则，以变形值小于 2.5 cm的格网占比(符合变形要求的格网个数/坐标定义覆盖区域的所有格网个数:m/n，if g(i)＜2.5)then m=m+1，n为所有格网数)作为参考观察指标，按逐步缩小搜索区间，直到理论最优值所在范围内达到预期设定要求，最后进行几种最优方案比选。本文将这种方法简称为格网占比法。

3.1 基于福州地区DEM提取500 m格网间距高程点

为了避免采用有限的控制点不够均匀，代表性不够，从而导致最优关系计算不够严密。因此，利用福州地区 1∶1万DEM数据，福州地区 1∶1万DEM数据格网尺寸为 12.5 m，高程精度 1 m。坐标系统为1980西安坐标系，1985黄海高程基准。根据福州山地地形的变化情况，选择边长 250 m格网来提取高程，基本可以满足计算需要。按一定的格网间距提取三维坐标点，进行最优关系计算。

3.2 最优中央子午线和抵偿高程面确定算法流程

(1)基础数据准备:按投影经度步长2分(约3 km)、抵偿面高程步长 50 m设置，生成所有计算方案表1、基于CGCS2000坐标系统的格网中心点数据表2。

(2)按遍历计算每个方案每个格网点的变形值△D，存入表3。

(3)完成一个方案的变形值计算后，将∑△D2和△D≤2.5 cm的格网占比结果计算存入表3的对应方案符合。

(4)最后遍历搜索法对比最优方案。

具体算法流程图如图1所示:

图1 计算流程图

3.3 最优中央子午线和抵偿高程面计算实验结果

以福州地区数据为例进行实验计算，选择精度118°40'～119°40'范围内，按 2 分步长递增，高程抵偿面从 0 m～600 m，按 20 m步长递增，对福州地区范围内的格网数据进行遍历计算。计算结果得到最优方案不同，如下:

(1)格网占比法:按△D≤2.5 cm的网格占比最大计算结果及对比数据。

变形值符合要求的格网占比计算结果及对比数据表 表1

变形值符合要求的格网占比模型计算结果表1列出最大的5个方案，最优的方案是中央子午线为119°34'，高程抵偿面为 50 m，符合每千米变形△D≤2.5 cm的网格占福州区域内格网总数的55.62%。但是每千米变形值平方和为39.78并非最小。空间分布专题图如图2(b)所示。

(2)最小二乘法:∑△D2最小值计算结果及对比数据。

每公里变形值平方和计算结果及对比数据表表2

∑△D2最小值模型计算结果表2列出最小的5个方案，最优的方案是中央子午线为119°30'，高程抵偿面为 150 m，符合每千米变形值平方和最小为38.00。但是符合每千米变形△D≤2.5 cm的网格占福州区域内格网总数为37.90%。空间分布专题图如图2(b)所示。

图2 格网占比法与最小二乘法最优方案对比图

从图2很明显看出左图的符合每千米变形△D≤2.5 cm的网格覆盖量比右图合理，左图福州市中心及沿海经济开发热点平原区域基本符合要求，而右图福州市中心城区大东边区域与沿海平潭岛及罗源区域大部分都不符合变形要求。显然，按格网占比法计算模型结果更为科学合理、更符合实际应用需要。按福州地方平面坐标系统选择的投影中央子午线为119°18'21″方 案［11］，计 算 结 果 为 符 合 每 千 米 变形△D≤2.5 cm的网格占比为46.04%，每千米变形值平方和值为64.15，并不是最优的方案。福州地区西北面基本为丘陵高山，许多高程基本在 400 m以上，经过遍历搜索计算，验证了该地形变化规律不适合采用高程抵偿面长度归化改正与投影变形进行抵消的计算模型。对于该高山区域的测图，只能采用抬高投影面的方法。

4 结论

本文采用DEM数据与地理信息系统计算方法，对CGCS2000城市平面坐标系统定义的两个计算模型进行研究，得出采用格网占比法方法得出的最优高程抵偿面与中央子午线的选择方案更具严密性，该方法对城市区域CGCS2000平面坐标系统定义具有较大的参考意义。

［1］陈俊勇，张鹏.国家测绘基准“十二五”重大项目的思考［J］.武汉大学学报·信息科学版，2009，34(10):1136～1138.

［2］陈俊勇.中国采用地心三维坐标系统对现有地图的影响［J］.测绘学报，2004，32(4):283 ～288.

［3］李江卫，白洁，郭际明等.基于CGCS2000的城市平面坐标系最佳选取［J］.城市勘测，2011，3(4):120～123.

［4］武丰雷，刘曦灿，王方等.济南市独立坐标系的改造探讨［J］.城市勘测，2010，2(2):105 ～106.

［5］CJJ T8－2011.城市测量规范［S］.

［6］吴云孙，晁定波，杨堂堂.浅谈广州市平面控制测量坐标系统的选取［J］.测绘通报，2005，23(8):42～43.

［7］陆鹏程，林冬伟.斜轴墨卡托投影模型及其应用分析［J］.铁道勘察，2010，12(4):26～29.

［8］吕忠刚，许世宁.关于抵偿高程面与移动中央子午线最佳选取问题的研究［J］.东北测绘，2003，25(2):3～8.

［9］谢征海，张泽烈.城市独立坐标系平面控制网之扩建及改造［J］.北京测绘，2000，23(4):22～23.

［10］孔建，姚宜斌，吴寒.整体最小二乘的迭代解法［J］.武汉大学学报·信息科学版，2010，35(6):711～714.

［11］高昭良.基于相似变换的城市坐标系统转换研究［J］.福州大学学报·自然科学版，2010，38(1):54～57.