一种基于DEM的城市平面坐标系统建立方法

王友昆，谢正明，张君华，余章蓉

(1.武汉大学 测绘学院，湖北 武汉 430079；2．昆明市测绘研究院，云南 昆明 650051； 3．云南国土资源职业学院，云南 昆明 652501；4．昆明理工大学 津桥学院，云南 昆明 650106)

城市平面坐标系统建立，在高斯统一3°带平面直角坐标系或高斯任意带平面直角坐标系无法满足长度变形规定时(≤25 mm/km)，可采用城市平均高程面作为投影面(抵偿高程面)。传统的投影面选择，一般基于控制点的平均高程或者控制点控制范围的平均高程确定[1-4]。该方法计算量小，速度快，但无法准确计算和分析符合长度变形规定的区域范围。随着计算机性能的提升，高分辨率数字高程模型(Digital Elevation Model,DEM)已经用于确定城市平均高程面[5]，长度变形计算[6]，以及投影面的可视化选择[7]。这些方法大多应用于平原等低海拔区域，主要涉及中央经线(投影带)的选取，没有涉及多个投影面的选择，也没有充分利用GIS空间分析和统计分析的方法对不同投影方案进行量化指标分析和统计。本文针对云南高海拔区域，城市地方坐标系的建立需考虑任意带和多个投影面的特点[8]，利用DEM成果进行长度变形计算、空间分析和可视化输出，并以昆明市城市平面坐标系建立为例进行详细说明。

1 长度变形原理

地面上实测的边长，内业计算时应经过高程归化至椭球面和高斯投影至高斯平面两个步骤，两者对于实测边长的影响分别为缩短和伸长[9]。《城市测量规范》规定，城市平面坐标系建立时，坐标反算的边长与实测边长应尽可能相符，即边长归算至参考椭球面上的高程归化和高斯正形投影的距离改化的总和(即长度变形)应小于25 mm/km，如式(1)所示。

(1)

式中:H为归化高程；R为地球平均曲率半径；ym为投影边长端点横坐标平均值。

按照式(1)计算分析，当城市地区高程大于160 m或其平面位置距离3°带中央子午线的东西方向距离(横坐标)大于45 km，其长度变形均超过规定的1/40 000[10]。

若式(1)无法满足长度变形要求，可采用任意带平面直角坐标系统。仍然无法满足时，可采用城市平均高程面作为投影面(抵偿高程面)，人为地改变归化高程使它与高斯投影的长度改化相抵偿，如式(2)所示。

(2)

式中:Hd为抵偿面高程值。

为提高长度变形计算的精确性，可采用投影点所在卯酉圈半径N替代式(2)中的地球平均曲率半径R，N的算式为：

(3)

式中:a为椭球长半轴；e为椭球第一偏心率；B为投影点纬度。

2 DEM计算方法

2.1 计算流程

长度变形原理揭示了城市平面坐标建立时需要满足的条件和计算方法。利用DEM成果进行城市平面坐标系建立，主要流程包括数据准备、预处理、变形计算、分析统计和成果输出等步骤(见图1)。

图1 计算流程图

2.2 数据准备

需要收集各类范围线、DEM成果数据和似大地水准面精化模型。范围线包括城市各级行政区划面、建成区和重点规划区范围，用于对DEM原始数据的裁剪和分析长度变形对这些区域的影响。DEM可利用1∶5万、1∶1万等高精度的DEM数据，若无法获取时可利用开源数据替代(分辨率优于100 m，精度优于20 m)，如航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission,SRTM)的DEM成果。似大地水准面精化模型主要用于将DEM的正常高成果改化为大地高成果。

2.3 预处理

1)统一坐标系统。将以上数据统一为2000国家大地坐标系或WGS84(两者椭球定义的参数差异极小，对长度变形影响不大)。

2)DEM裁剪。DEM按照格网存储，为减少计算量应利用市级行政区划对其进行裁剪，仅保留行政区划范围内的DEM成果，区域外格网属性值设置为NoData。

3)大地高计算。DEM成果为正常高，而长度变形计算式(1)和式(2)中H为大地高。我国西部区域高程异常达到-40 m，对长度变形计算影响较大。因此，需要利用区域似大地水准面模型，将DEM正常高改算为大地高，如式(4)所示。

H=h+ε.

(4)

式中:H为大地高；h为正常高；ε为高程异常。

如果没有可利用的区域似大地水准面模型，也可利用公开的全球超高阶地球重力场模型EGM2008，其中国区域精度优于1 m[11-13]，能够满足投影计算的要求。

4)高程统计。对区域范围的DEM高程进行统计，便于计算投影面的大致高程。

2.4 变形计算

变形计算遵循图 1中方案1至方案6的顺序，依次尝试计算，直至找到符合长度形变要求的城市平面坐标系统。如首先采用方案1进行变形计算，即采用统一3°带的平面直角坐标系进行DEM变形计算，如果方案1不符合变形要求，则尝试方案2，以此类推。是否符合变形要求，根据分析统计结果来判定。

1)高斯投影。按照长度变形计算的原理，需要对DEM成果进行高斯投影，中央经线L0按照图1的变形计算顺序设置，计算DEM格网的横坐标y(应减去加常数500 km)。

2)长度变形计算。DEM为规则格网，格网自带坐标及高程，为快速计算每一格网的长度变形提供了便利。横坐标y和卯酉圈半径N可用一维数组存储，大地高可用二维数组存储。利用成熟的矩阵库(如Numpy)，可大大提高长度变形的计算效率。

3)格网标识。对长度变形值小于1/40 000的格网，标识为0，否则标识为1，并将标识后的DEM格网输出为栅格，用于后续的分析和统计。

2.5 分析统计

分析统计主要是对DEM变形标识的格网，利用GIS原理对其分析和统计，作为城市平面坐标系统建立投影带和投影面确定的依据。

1)格网统计。通过统计符合长度变形要求的格网比例来衡量区域长度变形是否符合要求。假设占比达到预定的比例(如90%以上，即大部分区域长度变形满足1/40 000的要求)，则认为该方案符合要求。

2)栅矢转换。为直观地展示符合长度变形要求的范围，以及多个投影面或关注区域(城市重点经济建设区域)的长度变形情况，需要将变形标识为0的栅格转换为矢量多边形。

3)叠置分析。如果城市平面坐标系涉及多个投影面，应利用县(区)行政区域范围或重点经济建设区域范围，同符合长度变形要求的矢量面进行叠置分析，以此作为选择投影面的依据。对于某投影带，n个投影面计算的符合长度变形的多边形记为prji(i=1,2,…,n)，其叠置分析主要步骤如下：

a)将prji中占比面积最大的作为主投影面，记为prjmax；

b)遍历prji(i=1,2,…,n,i≠max)，分别同prjmax叠置分析，采用prji删除prjmax，即第i个投影面排除主投影面后符合变形要求的范围，记为Erasei；

c)将Erasei面积最大的作为次投影面prjmax-1；

d)重复步骤b)和c)，分别确定出后续的投影面prjmax-2，prjmax-3…prj1。

以上叠置分析过程中，还应同时考虑投影面的数量、各县(市、区)建成区以及重点经济建设区范围，以及单个符合长度变形区域范围面积是否最大化等因素。

2.6 成果输出

将DEM计算过程的所有成果，包括裁剪和投影后的DEM成果，各类矢量多边形范围线，不同计算方案的长度变形计算结果和统计分析结果，以数据库和图、表的方式输出，用于成果报告的编制。

3 实例应用

3.1 需求分析

由于历史原因，昆明市1987昆明坐标系、2004昆明坐标系以及2013昆明坐标系等多个坐标系统并存，给城市建设和成果共享带来诸多不便。按照国家及昆明市政府相关要求，昆明市启动了建立基于CGCS2000椭球的昆明2000坐标系。昆明市辖7个区、3个县、代管1个县级市和3个自治县，总面积约2.1万km2，东西跨度150 km，南北跨度240 km，海拔范围746～4 247 m，地势地貌呈现北高南低、分布不均匀，且各县(市、区)建成区高程各不相同(见图2)，黑色实线为昆明市主城区及各县(市、区)建成区范围。

图2 昆明市DEM

3.2 计算过程

1)数据准备。DEM采用STRM1数据，分辨率为30 m，高程精度为1 m；范围线收集了昆明市及各县(市、区)行政区域范围，以及各县(市、区)建成区范围线；区域似大地水准面精化模型利用EGM2008成果；以上成果坐标系统统一为CGCS2000。

2)数据预处理。DEM经过裁剪、大地高计算、高程统计[14]等预处理。如表1所示，昆明市高程值跨度较大，1 500～2 500 m占86.4%，高程大于2 500 m的比例达到9.7%。

表1 昆明市DEM高程值区间统计

3)变形计算和分析统计。按照图1变形计算的顺序，经试算和分析，采用统一3°带等计算方案(方案1—方案4)均无法满足变形要求，可采用多个投影面的任意带(单中央经线)计算方案。为便于后期成果的接边和管理，方案优先选择顺序：①满足昆明市主城区长度变形要求；②满足各县(市、区)建成区长度变形要求；③符合要求的长度变形区域范围能够连片且面积最大化；④投影面数量尽量少；⑤多个投影面符合变形要求的区域面积之和达到一定的比例(如占昆明市域面积90%)。按照以上原则，利用程序对不同投影带(102°42′～103°09′为区间，经差1′为间隔)、不同投影面(700～4 300 m为区间，10 m为间隔)的长度变形结果分析，选择最优的投影带和投影面(见图3)。由于涉及大量空间分析计算，为减少计算量，初算可采用经差6′或3′，投影面间隔100 m或50 m，确定大致区间范围后，再缩小计算间隔，减少冗余计算，提高计算效率。

图3 投影面面积对比图

图3显示不同投影带满足原则：①符合变形区域范围的最大面积(虚线星号)和单个连片最大面积(实线圆点)对比图。若仅考虑覆盖面积最大，随着中央经线增大，覆盖范围逐渐增大，L0=102°56′时达到最大值8 615.88 km2(单个连片最大面积仅为5 598.20 km2),之后逐渐减小。若综合考虑原则③(单个连片最大面积)，L0=102°47′时达到最大值8 025.46 km2。此时，投影面为1 890 m，作为第1投影面(记为Prj1)更为合理。之后，将其他投影面多边形同Prj1进行叠置分析(删减策略)，按照原则②③④进行逐一分析和比选，确定第2投影面。以此类推，直至多个投影面覆盖面积之和满足原则⑤。

4)成果输出。利用ArcGIS制图模板[15]进行成果输出，只需修改各图层数据源即可达到统一的输出标准，便于对不同计算方案结果比选。图 4显示了满足原则①②③④不同投影带(102°42′～102°57′，3′为间隔)的计算结果。红色实线为昆明市域范围，黑色实线为各县(市、区)建成区范围，黑色虚线为投影中央经线，绿色、蓝色、粉色填充区域分别为第1、第2、第3投影面覆盖范围。图4(d)—图4(f)可见，第1投影面覆盖面积逐渐增大，但单个连片面积逐渐减小，结果同图 3是一致的。

如果考虑原则⑤，需设置第4和第5投影面，如表 2所示，此时符合长度变形要求的总面积达到19 256.36 km2，占市域面积的91.7%。未覆盖区域海拔大多在3 000 m以上，且位于昆明市生态红线保护范围内，经济建设较少，可在实际工程建设时再行考虑。

图4 投影面范围分布图

表2 投影面面积指标统计表

4 结束语

本文利用DEM数据易于读取、计算和分析的特点，应用于高原城市平面坐标系统建立时投影带和投影面的选择，并充分利用GIS空间分析的优势，对长度变形区域和行政区域、重点建设区域进行空间叠置分析，定量地分析和统计不同投影方案的计算结果，并实现了成果的可视化输出。同时，实现了软件程序化批量计算，大大提高了计算的效率。这种方法使计算者能够直观、便利、及时地调整计算方案，较传统方法更加有效和可靠，能更优地选择出投影带和投影面。同理，该方法还可应用于大型区域工程和线状工程的平面坐标系建立。