CN111462322A - 一种基于dem的城市平面坐标系统建立方法 - Google Patents

Abstract

本发明公开了一种基于DEM的城市平面坐标系统建立方法，先进行统一3°投影带形变，若不能满足形变要求，则进行抵偿高程面统一3°带形变，还不能满足形变要求，则以此类推进行任意投影带形变、抵偿高程面的任意投影带形变、多个抵偿高程面的任意投影带形变、多个任意投影带形变，直至投影变形满足要求。本方法与传统的计算方法相比，计算效率更高、结果更可靠，同时可以定量计算区域满足该投影面形变要求所能控制的面积，结果可视化。

Description

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

技术领域

本发明属于测绘科学与工程技术领域，涉及一种基于DEM的城市平面坐标系统建立方法。

背景技术

坐标系统对测绘工作的开展和质量控制具有十分重要的作用，在城市规划、建设、测绘管理等工作中发挥着不可替代的作用。加强城市坐标系统的建设，有助于城市的规划发展及建设。坐标系统的建设，是测绘工作的基础。根据CJJ/T8《城市测量规范》的要求，城市平面控制网要采用统一坐标系统，必须具备以下条件：

(1)城市中心地区位于高斯正形投影统一3°带中央子午线附近。

(2)城市平均高程面必须接近国家参考椭球体面或平均海水面。

(3)城市所在的地区的国家网精度高于城市首级网的精度。

同时满足以上条件的城市并不多，有一些城市位于3°带中央子午线附近，但由于其城市平均高程太高，或者城市东西宽度较大远离中央子午线的区域投影变形也超出了限差的要求。城市地方平面坐标系统的长度变形应不大于2.5cm/km要求，即投影长度变形要小于1/40000。因此，当长度变形值大于2.5cm/km时，可根据具体情况与要求建议按照以下顺序选择坐标系统：

(1)抵偿高程面上的高斯正形投影3°带平面直角坐标系统；

(2)高斯正形投影任意带平面直角坐标系统；

(3)高斯正形投影任意带分带投影坐标系统；

(4)假定平面直角坐标系统；

投影变形控制的目的是使两点间坐标反算的长度和实测长度之比(投影长度比)接近于1，这样在使用这两点数据的时候可以不用进行任何归算，则可以直接利用实测长度进行城市建设中的大比例尺地形图测制、市政工程施工放样等有关建设活动。

因此，根据城市所处地理位置及城市建设发展需要，绝大部分城市都建立有自己的城市独立坐标系。现行的城市坐标系是在西安1980坐标系的基础上建设起来的，如今新旧坐标系过渡时间已结束，国家已明确要求，政府部门的地理信息成果及数据均需以2000国家大地坐标系提供使用。而现在须全面启用CGCS2000国家大地坐标系，则城市坐标系需要进一步调整完善，以满足国家要求及城市建设发展需要。

这就要求原先的城市坐标系统须在CGCS2000椭球下重新设计和建立。而城市坐标系统的建设，离不开投影带的分析。此外，大型工程建设项目，如水库、机场、公路等工程项目的建设也涉及需要建立其独立坐标系统，如何基于2000大地坐标系选择中央经线和投影面也是需要面临的问题。

现有的坐标系统中央经线和投影面的选择一般基于单点、控制网及经验，投影形变量大小主要顾及城市经济发展活跃区域。随着城市扩张及城乡一体化统筹管理，传统方法设计的地方坐标已不能适应城市发展的需要。例如云南省高原城市坐标系的选择，投影变形的大小在山区或城市边缘地区往往超过了规范不大于2.5cm/km的要求，形变差异的大小，以及超过规范要求的区域并不能量化指标地体现出来。因此，有必要重新设计一套城市或区域基于2000国家大地坐标系建立的新方法，能够实现参数化、自动化、可视化的计算。通过对高原地区城市坐标系统投影带及投影面可视化确定研究，可以帮助建立新坐标系下的城市坐标系，并直观、准确、量化的进行城市坐标系的精度分析及研究。同时，该方法同样也可以应用于大型工程建设、线性工程建设的独立坐标系选择和可视化计算。

当城市地区高程大于160m或其平面位置距离3°带中央子午线的东西方向距离(横坐标)大于45km，其长度变形均超过规定的1/40000，这时就应该采取适当的措施。传统的独立坐标系建立，投影面的选择一般基于联测控制点的平均高程或者控制点控制范围的平均高程进行确定。这种方法主要根据高斯投影坐标及投影面高程计算抵偿地带高程和相应横坐标区间的关系，以及投影变形限差要求(2.5cm/km)，大致确定该地区的投影中央经线和投影面，能够满足区域大部分的计算要求。

李祖峰的“轴线投影变形最小任意带高斯正形投影参数确定”中，针对任意带高斯正形投影，将投影参考位置选择在测区中央并不能使投影变形在该准则下达到最优化的问题，提出了一种基于选定轴线投影变形最小准则确定任意带投影参考位置的方法、并确定了中央子午线至测区投影参考位置的方法。

这些方法，一般只考虑了投影范围内的主要控制点高程(海拔高)，以及大致的测区平均高程，无法对整个区域的投影变形情况，如精度、分布、比例情况进行详细的掌握和调整。传统的方法，只能满足范围内大部分的投影变形要求，无法进行大区域多投影带、多投影面计算，无法实现自动计算。此外，投影变形计算时，传统方法一般直接利用海拔高计算，高程异常较大的区域，对计算精度也有一定的影响。最重要的是对于范围较大的计算区域，如多投影带、多投影面计算时，传统的方法无法准确、直观的选择投影面，无法定量的计算各投影面所能控制的区域面积；无法对确定的投影面所控制的范围进行可视化、直观的表达。

由于以上原因不能采用统一3°带高斯正形投影平面直角坐标系或者抵偿坐标系时，则可以采用任意带高斯正形投影平面直角坐标系，并用城市平均高程面进行高程归化，以减小长度变形。如果城市的东西方向跨度很大，当高斯正形投影任意带投影的一个投影带不能保证将投影长度变形控制在2.5cm/km以内时，可以采用分带投影的方法。

以上则为传统的城市(区域/独立)坐标建立时投影带选择及投影面选择的基本方法，这种方法主要根据高斯投影坐标及投影面高程计算抵偿地带高程和相应横坐标区间的关系，以及投影变形限差要求(2.5cm/km),大致确定该地区的投影中央经线和投影面，能够满足区域大部分的计算要求。当城市区域范围较大，以及高海拔区域时，涉及多投影带或者多投影面选择时，如何评价和计算不同投影带或投影面的形变情况指标、控制范围时，传统的方法显得无能为力。

传统的独立坐标系建立，投影面的选择一般基于联测控制点的平均高程或者控制点控制范围的平均高程进行确定。这种方法主要根据高斯投影坐标及投影面高程计算抵偿地带高程和相应横坐标区间的关系，以及投影变形限差要求(2.5cm/km),大致确定该地区的投影中央经线和投影面，能够满足区域大部分的计算要求。

综上所述，在实现本发明的过程中，发明人发现现有技术中存在如下缺点：

1、虽然传统的方法能够满足区域大部分投影变形(2.5cm/km)的要求，但无法定量计算区域满足该投影面形变要求所能控制的面积，特别是在高海拔地区或地形起伏变化较大的区域。

2、当城市行政区域范围较大，单个投影面无法满足城市投影变形要求(2.5cm/km)时，传统的方法无法准确、直观的选择投影面；同时，也是无法定量的计算各投影面所能控制的区域面积；无法对确定的投影面所控制的范围进行可视化、直观的表达。

3、传统的计算方法，投影变形计算时，实测边长归算至参考椭球面时一般采用海拔高，而不是采用大地高进行计算，但在中国西部区域高程异常较大(小于-30m及以下)，对投影变形的影响较大，对计算的结果有较大影响，因此有必要利用大地高来进行投影变形的计算。

4、传统的计算方法，计算量大，计算效率低，量化指标不明显，特别是在城市大范围、高程起伏较大的地区，无法进行大区域多投影带、多投影面计算，无法实现自动计算。

发明内容

为实现上述目的，本发明提供一种基于DEM的城市平面坐标系统建立方法，与传统的计算方法相比，计算效率更高、结果更可靠，同时可以定量计算区域满足该投影面形变要求所能控制的面积，结果可视化。

本发明所采用的技术方案是，一种基于DEM的城市平面坐标系统建立方法，按照如下方法进行：

步骤S1：坐标读取和转换

读取DEM数据，通过DEM数据计算每一格网的纬度B和经度L，得到大地坐标为(B,L)；

步骤S2：高程异常和大地高计算

通过大地坐标(B,L)，求得该点对应的高程异常ε，然后根据高程异常ε计算大地高H；

步骤S3：统一3°投影带形变计算

步骤S3.1根据区域范围内的最小经度minL和最大经度maxL计算中央子午线L₀,然后根据中央子午线L₀计算其所在统一3°带的中央子午线L'₀；

步骤S3.2按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x,y)，中央子午线取L'₀；

步骤S3.3计算每个格网的纬度B对应的卯酉圈曲率半径N；

步骤S3.4利用每个格网的横坐标y、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ，见公式(5)：

步骤S3.5利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影统一3°带平面坐标系”；

步骤S4：抵偿高程面统一3°带形变计算

若按照步骤S3计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则使用抵偿高程面的高斯正形投影统一3°带平面坐标系，计算过程如下：

步骤S4.1利用区域范围大地高的最小值minH和最大值maxH，计算抵偿高程面的区域范围[H_min，H_max],H_min为抵偿高程面的最小值，H_max为抵偿高程面的最大值，见公式(6)和(7)：

式中，step为步长；

步骤S4.2计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(8)：

式中，投影面高程H_d为抵偿高程面的区间范围[H_min，H_max]以step为步长的取值；

步骤S4.3利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的高斯正形投影统一3°带平面坐标系”；

步骤S5：任意投影带形变计算

若按照步骤S4计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则使用任意带高斯正形投影平面坐标系，计算过程如下：

步骤S5.1按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀；

步骤S5.2利用每个格网的横坐标y1、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ₁，见公式(9)：

步骤S5.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影任意带平面坐标系”；

步骤S6：抵偿高程面的任意投影带形变计算

若按照步骤S5计算的大部分区域无法满足投影变形的要求，即大部分区域形变大于2.5cm/km，则使用抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S6.1按照高斯正形投影公式，将每个格网大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀；

步骤S6.2重复步骤S4.1和步骤S4.2，将公式(8)中y变为y₁，计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(10)：

步骤S6.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的任意带高斯正形投影平面坐标系”；

步骤S7：多个抵偿高程面的任意投影带形变计算

若步骤S6计算的大部分区域无法满足投影形变的要求，即大部分区域形变大于2.5cm/km，则使用多个抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S7.1重复步骤S6.1、步骤S6.2和步骤S6.3，然后利用步骤S10将格网形变标识区域转化为带属性0和1的矢量多边形；

步骤S7.2利用步骤S11或步骤S12进行叠置分析计算，得到符合要求的投影面高程，并计算其所能覆盖的区域面积和比例，即采用“多个抵偿高程面的任意带高斯正形投影平面坐标系”；

步骤S8：多个任意投影带形变计算

若步骤S7计算的大部分区域无法满足投影形变的要求，即大部分区域形变大于2.5cm/km，则使用多个任意投影带形变计算，计算过程如下：

步骤S8.1根据投影范围的经度跨度，选择2个三等分处的经线为中央经线L'_i(i＝1,2)；

步骤S8.2按照步骤S3至步骤S7的顺序，重复计算；

步骤S8.3若步骤S8.2计算的投影变形无法满足要求，则选择3个四等分处的经线为中央经线L'_i(i＝1,2,3)，以此类推，重复步骤S8.2和步骤S8.3，直至计算的投影变形满足要求为止；

步骤S9：格网形变标识及占比统计计算

步骤S9.1对各格网的综合投影变形δ与限差δ_限进行比较和标识，若δ≤δ_限，则标识该格网形变值M＝0，否则M＝1；

步骤S9.2计算标识为0的格网数量同所有格网数的比例，若这一比例达到预定的比例，则投影变形符合要求；

步骤S10：格网形变标识区域转化为矢量多边形

步骤S10.1分别将标识为0和1的格网，输出为栅格文件；

步骤S10.2分别将栅格文件转化为矢量文件，即分别得到标识为0和1覆盖范围的多边形，属性为0的多边形即满足形变要求的范围，属性为1的多边形即不满足形变要求的范围；

步骤S11：多个分区范围的多边形叠置分析计算

步骤S11.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min,H_max]共计算n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)；

步骤S11.2读取需要考虑投影形变的多个分区范围面，记为subA_j(j＝1，2，...，k)；

步骤S11.3遍历各分区范围面subA_j，并分别对n个投影面的变形多边形范围prj_i进行空间叠置分析，叠置分析采用相交的策略，将prj_i与subA_j的相交面积记为Inter sec tion_i,并计算相交面同对应分区的范围面的面积占比，取占比最大的投影面高程作为该行政区范围的投影面；

步骤S11.4计算Inter sec t ion_i同各分区的面积占比，统计各投影面能达到的形变要求比例；

步骤S12：一个分区范围的多边形叠置分析计算

步骤S12.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min,H_max]共计算了n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)；

步骤S12.2计算各投影面属性为0的多边形面积占整个投影区域的比例，以占比最大的矢量多边形所在的投影面，作为主投影面，记为prj_max；

步骤S12.3遍历各投影面prj_i(i＝1,2,...,n,i≠max)，分别同主投影面的多边形prj_max进行空间叠置分析，叠置分析采用删除策略，即prj_i面删除面prj_max，剩余的面则为第i个投影面排除主投影面后符合形变要求的面积，记为Erase_i；

步骤S12.4对删除后的面Erase_i按照面积从大到小进行排序，面积最大的多边形确定为次投影面prj_max-1；

步骤S12.5重复步骤S12.3和步骤S12.4，其中遍历各投影面prj_i(i＝1,2,...,n,i≠max,max-1)，prj_max-1替代prj_max作为主投影面，并以此类推，分别确定出后续的投影面prj_max-2，prj_max-3...prj₁；

步骤S12.6根据叠置分析计算的投影面prj_max-j(j＝0,1,...,i-1)的面积，分别计算所占整个分区面积的比例大小，确定整个分区投影面的个数。

进一步的，步骤S1中所述纬度B和经度L的计算过程为：通过DEM数据的左上角经纬度及格网分辨率大小，计算每一格网的纬度B和经度L，即可得到大地坐标(B,L)；如果DEM数据的参考椭球不是城市坐标系统参考椭球，采用布尔沙模型将计算的纬度B和经度L转换为对应参考椭球的大地坐标(B,L)。

进一步的，步骤S1中所述DEM数据为一个二维数组，行号代表纬度的索引，列号代表经度的索引，数组内的取值为该格网的平均海拔高程值h。

进一步的，步骤S2中所述高程异常ε求解过程为：如果有可利用的区域似大地水准面模型，通过大地坐标(B,L)，查询该点对应的高程异常ε值；如果没有可利用的区域似大地水准面模型，则利用全球超高阶地球重力场模型，通过大地坐标(B,L)，计算对应的高程异常ε值。

进一步的，步骤S2中所述大地高H计算见公式(1)：

H＝h+ε (1)

式中，H为大地高，h为海拔高程值，ε为高程异常。

进一步的，步骤S3.1中所述中央子午线L₀的计算见公式(2)，所述所在统一3°带的中央子午线L'₀,见公式(3)：

L₀＝(min L+max L)/2 (2)

L'₀＝(int)(L₀/3)*3 (3)

进一步的，步骤S3.3中所述卯酉圈曲率半径N计算见公式(4)：

式中，a为椭球长半轴，e为椭球第一偏心率。

本发明的有益效果是：

1、本发明充分利用DEM规则格网数据易于读取、计算和分析的特点，基于DEM进行城市平面坐标系统的建立适用于大范围的城市投影形变分析，较传统的单点计算效率更高、结果更可靠。

2、本发明提出了一整套城市坐标系统建立不同形变方案的计算方法，涵盖了不同投影形变的场景，能够适用于全国各个城市的平面坐标系建立时投影形变的计算，特别是解决了需要确定多个投影抵偿面或多个中央经线的高海拔地区或地形起伏变化较大的区域。

3、本发明充分利用空间分析的原理，对输出的形变多边形范围、各行政区域范围进行空间叠置分析，较传统的方法更能有效、快速地计算和统计不同方案场景的形变量，定量地计算各投影面所能控制的区域面积和比例。

4、本发明将DEM计算分析的投影形变量以栅格图形或矢量多边形可视的方式进行显示，让计算者能够直观地对不同投影方案的形变结果同实际地形情况(DEM)、行政区划范围等要素进行叠加和比选，更加有利于投影面的合理选择。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例城市平面坐标系统建立的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例处理对象为区域范围内的DEM数据(分辨率优于100m，精度优于20m)，坐标系统为2000国家大地坐标系或WGS84坐标系统，高程为正常高(海拔高)，数据格式为常用的DEM数据格式。一种基于DEM的城市平面坐标系统建立方法，具体实施流程如下(如图1)：

步骤S1:坐标读取和转换

步骤S1.1读取数字高程模型(DEM)数据。DEM数据读取为一个二维数组，行号代表纬度的索引，列号代表经度的索引，数组内的取值为该格网的海拔高程值h。

步骤S1.2计算格网纬度和经度。通过DEM数据的左上角经纬度及格网分辨率大小，计算每一格网的纬度B和经度L，则大地坐标为(B,L)。

步骤S1.3坐标转换。如果DEM数据的参考椭球不是城市坐标系统参考椭球，采用布尔沙模型将计算的纬度B和经度L转换为对应参考椭球的大地坐标(B,L)。

步骤S2：高程异常和大地高计算

步骤S2.1如果有可利用的区域似大地水准面模型，通过大地坐标(B,L)，查询该点对应的高程异常ε值；如果没有可利用的区域似大地水准面模型，则利用全球超高阶地球重力场模型(如EGM2008模型)，通过大地坐标(B,L)，计算对应的高程异常ε值。

步骤S2.2计算大地高,见公式(1)：

H＝h+ε (1)

式中，H为大地高，h为海拔高程值(正常高)，ε为高程异常。

步骤S3：统一3°投影带形变计算

步骤S3.1根据区域范围内的最小经度min L和最大经度max L计算中央子午线L₀,见公式(2)；再计算其所在统一3°带的中央子午线L'₀,见公式(3)。

L₀＝(min L+max L)/2 (2)

L'₀＝(int)(L₀/3)*3 (3)

步骤S3.2按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x,y)，中央子午线取L'₀。

步骤S3.3为保证投影综合长度变形更加精确，采用卯酉圈曲率半径N代替综合形变公式中的地球平均曲率半径，因此需要计算每个格网纬度B对应的卯酉圈曲率半径N,见公式(4)。

式中，a为椭球长半轴，e为椭球第一偏心率。

步骤S3.4利用每个格网的横坐标y、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ，见公式(5)。

步骤S3.5利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影统一3°带平面坐标系”。

步骤S4:抵偿高程面统一3°带形变计算

若按照步骤S3计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则尝试使用抵偿高程面的高斯正形投影统一3°带平面坐标系，计算过程如下：

步骤S4.1利用区域范围大地高的最小值minH和最大值maxH，计算抵偿高程面的区域范围[H_min，H_max],H_min为抵偿高程面的最小值，H_max为抵偿高程面的最大值，见公式(6)和(7)。

式中，step为步长，根据计算的精细程度一般取50m或100m。

步骤S4.2计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(8)：

式中，投影面高程H_d为抵偿高程面的区间范围[H_min，H_max]以step为步长的取值。

步骤S4.3利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的高斯正形投影统一3°带平面坐标系”。

步骤S5：任意投影带形变计算

若按照步骤S4计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则尝试使用任意带高斯正形投影平面坐标系，计算过程如下：

步骤S5.1按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀。

步骤S5.2利用每个格网的横坐标y1、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ₁，见公式(9)：

步骤S5.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影任意带平面坐标系”。

步骤S6：抵偿高程面的任意投影带形变计算

若按照步骤S5计算的大部分区域无法满足投影变形的要求，即大部分区域形变大于2.5cm/km，则尝试使用抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S6.1按照高斯正形投影公式，将每个格网大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀。

步骤S6.2重复步骤S4.1和步骤S4.2，将公式(8)中y变为y₁，计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(10)：

步骤S6.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的任意带高斯正形投影平面坐标系”。

步骤S7：多个抵偿高程面的任意投影带形变计算

若步骤S6计算的大部分区域无法满足投影形变的要求，则尝试使用多个抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S7.1重复步骤S6.1、步骤S6.2和步骤S6.3，然后利用步骤S10将格网形变标识区域转化为带属性0和1的矢量多边形。

步骤S7.2利用步骤S11或步骤S12进行叠置分析计算，得到符合要求的投影面高程，并计算其所能覆盖的区域面积和比例，即采用“多个抵偿高程面的任意带高斯正形投影平面坐标系”。

步骤S8：多个任意投影带形变计算

若步骤S7计算的大部分区域无法满足投影形变的要求，则尝试使用多个任意投影带形变计算，计算过程如下：

步骤S8.1根据投影范围的经度跨度，选择2个三等分处的经线为中央经线L'_i(i＝1,2)。

步骤S8.2按照步骤S3至步骤S7的顺序，重复计算。

步骤S8.3若步骤S8.2计算的投影变形无法满足要求，则选择3个四等分处的经线为中央经线L'_i(i＝1,2,3)，以此类推，重复步骤S8.2和步骤S8.3，直至计算的投影变形满足要求为止。城市的东西跨度，一般2个合适的中央经线均能满足要求，当为线性工程时，可能存在3个及以上的投影带。

步骤S9：格网形变标识及占比统计计算

步骤S9.1对各格网的综合投影变形δ与限差δ_限(规范要求一般取值

即2.5cm/km——每千米形变量不大于2.5cm)进行比较和标识，若δ≤δ_限，则标识该格网形变值M＝0，否则M＝1。

步骤S9.2计算标识为0的格网数量同所有格网数的比例，若这一比例达到预定的比例(如95％以上，即大部分区域投影综合形变均满足规范要求)，则认为该方案的投影变形符合要求。

步骤S10：格网形变标识区域转化为矢量多边形

步骤S10.1分别将标识为0和1的格网，输出为栅格文件。

步骤S10.2分别将栅格文件转化为矢量文件，即分别得到标识为0和1覆盖范围的多边形，属性为0的多边形即满足形变要求的范围，属性为1的多边形即不满足形变要求的范围。

步骤S11：多个分区范围的多边形叠置分析计算

步骤S11.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min,H_max]共计算n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)。

步骤S11.2读取需要考虑投影形变的多个分区范围面(例如昆明市单个投影面无法满足整个昆明市投影变形需求的时候，则需要考虑分区域投影变形计算分析，这个时候则需要各分区的范围面，用于空间分析计算)，数量记为k(如城市所辖县、区行政区划面或特殊建设区域的范围)，记为subA_j(j＝1，2，...，k)。

步骤S11.3遍历各分区范围面subA_j，并分别对n个投影面的变形多边形prj_i进行空间叠置分析，叠置分析采用相交的策略，将prj_i与subA_j的相交面积记为Inter sec t ion_i,并计算相交面同对应分区的范围面的面积占比，取占比最大的投影面高程作为该行政区范围的投影面。

步骤S11.4计算Inter sec t ion_i同各分区的面积占比，统计各投影面能达到的形变要求比例(例如，石林县的县城大地高范围为1500-2000m，当投影面取1600m时，达到形变要求的Inter sec t ion_i面积占比达到90％，则我们可以认为在投影面1600m时能够满足石林县的大部分区域的投影形变要求)。

步骤S12：一个分区范围的多边形叠置分析计算

步骤S12.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min,H_max]共计算了n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)。

步骤S12.2计算各投影面属性为0的多边形面积占整个投影区域的比例，以占比最大的矢量多边形所在的投影面，作为主投影面，记为prj_max。

步骤S12.3遍历各投影面prj_i(i＝1,2,...,n,i≠max)，分别同主投影面的多边形prj_max进行空间叠置分析，叠置分析采用删除策略，即prj_i面删除面prj_max，剩余的面则为第i个投影面排除主投影面后符合形变要求的面积，记为Erase_i。

步骤S12.4对删除后的面Erase_i按照面积从大到小进行排序，面积最大的多边形确定为次投影面prj_max-1。

步骤S12.5重复步骤S12.3和步骤S12.4，其中遍历各投影面prj_i(i＝1,2,...,n,i≠max,max-1)，prj_max-1替代prj_max作为主投影面，并以此类推，分别确定出后续的投影面prj_max-2，prj_max-3...prj₁。

步骤S12.6根据叠置分析计算的投影面prj_max-j(j＝0,1,...,i-1)的面积，分别计算所占整个分区面积的比例大小，确定整个分区投影面的个数(例如，如prj_max占60％，prj_max-1占20％，prj_max-2占15％，那么我们认为选择这三个投影面所能覆盖整个区域的面积达到了95％，则可以满足整个区域的投影要求)。

上述过程可以通过平台软件实现可视化的计算和操作，特别在多投影带、多投影面，多个计算区域范围时，能够快速准确的计算，将计算的结果(符合形变要求的区域，格网属性＝0；不符合形变要求的区域，格网属性＝1)呈现在计算者面前。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims (7)

1.一种基于DEM的城市平面坐标系统建立方法，其特征在于，按照如下方法进行：

步骤S1：坐标读取和转换

读取DEM数据，通过DEM数据计算每一格网的纬度B和经度L，得到大地坐标为(B,L)；

步骤S2：高程异常和大地高计算

通过大地坐标(B,L)，求得该点对应的高程异常ε，然后根据高程异常ε计算大地高H；

步骤S3：统一3°投影带形变计算

步骤S3.1根据区域范围内的最小经度minL和最大经度maxL计算中央子午线L₀,然后根据中央子午线L₀计算其所在统一3°带的中央子午线L'₀；

步骤S3.2按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x,y)，中央子午线取L'₀；

步骤S3.3计算每个格网的纬度B对应的卯酉圈曲率半径N；

步骤S3.4利用每个格网的横坐标y、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ，见公式(5)：

步骤S3.5利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影统一3°带平面坐标系”；

步骤S4：抵偿高程面统一3°带形变计算

若按照步骤S3计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则使用抵偿高程面的高斯正形投影统一3°带平面坐标系，计算过程如下：

步骤S4.1利用区域范围大地高的最小值minH和最大值maxH，计算抵偿高程面的区域范围[H_min，H_max],H_min为抵偿高程面的最小值，H_max为抵偿高程面的最大值，见公式(6)和(7)：

式中，step为步长；

步骤S4.2计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(8)：

式中，投影面高程H_d为抵偿高程面的区间范围[H_min，H_max]以step为步长的取值；

步骤S4.3利用过程步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的高斯正形投影统一3°带平面坐标系”；

步骤S5：任意投影带形变计算

若按照步骤S4计算的大部分区域投影变形无法满足要求，即大部分区域形变大于2.5cm/km，则使用任意带高斯正形投影平面坐标系，计算过程如下：

步骤S5.1按照高斯正形投影，将每个格网的大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀；

步骤S5.2利用每个格网的横坐标y1、大地高H和卯酉圈曲率半径N,计算每个格网的综合投影变形δ₁，见公式(9)：

步骤S5.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“高斯正形投影任意带平面坐标系”；

步骤S6：抵偿高程面的任意投影带形变计算

若按照步骤S5计算的大部分区域无法满足投影变形的要求，即大部分区域形变大于2.5cm/km，则使用抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S6.1按照高斯正形投影公式，将每个格网大地坐标(B,L)高斯正算为高斯平面坐标(x1,y1)，中央子午线取L₀；

步骤S6.2重复步骤S4.1和步骤S4.2，将公式(8)中y变为y₁，计算不同投影面高程H_d下的每个格网的投影形变量δ_d，见公式(10)：

步骤S6.3利用步骤S9进行标识分析计算，符合要求，则结束计算，即采用“抵偿高程面的任意带高斯正形投影平面坐标系”；

步骤S7：多个抵偿高程面的任意投影带形变计算

若步骤S6计算的大部分区域无法满足投影形变的要求，即大部分区域形变大于2.5cm/km，则使用多个抵偿高程面的任意带高斯正形投影平面坐标系，计算过程如下：

步骤S7.1重复步骤S6.1、步骤S6.2和步骤S6.3，然后利用步骤S10将格网形变标识区域转化为带属性0和1的矢量多边形；

步骤S7.2利用步骤S11或步骤S12进行叠置分析计算，得到符合要求的投影面高程，并计算其所能覆盖的区域面积和比例，即采用“多个抵偿高程面的任意带高斯正形投影平面坐标系”；

步骤S8：多个任意投影带形变计算

若步骤S7计算的大部分区域无法满足投影形变的要求，即大部分区域形变大于2.5cm/km，则使用多个任意投影带形变计算，计算过程如下：

步骤S8.1根据投影范围的经度跨度，选择2个三等分处的经线为中央经线L'_i(i＝1,2)；

步骤S8.2按照步骤S3至步骤S7的顺序，重复计算；

步骤S8.3若步骤S8.2计算的投影变形无法满足要求，则选择3个四等分处的经线为中央经线L'_i(i＝1,2,3)，以此类推，重复步骤S8.2和步骤S8.3，直至计算的投影变形满足要求为止；

步骤S9：格网形变标识及占比统计计算

步骤S9.1对各格网的综合投影变形δ与限差δ_限进行比较和标识，若δ≤δ_限，则标识该格网形变值M＝0，否则M＝1；

步骤S9.2计算标识为0的格网数量同所有格网数的比例，若这一比例达到预定的比例，则投影变形符合要求；

步骤S10：格网形变标识区域转化为矢量多边形

步骤S10.1分别将标识为0和1的格网，输出为栅格文件；

步骤S10.2分别将栅格文件转化为矢量文件，即分别得到标识为0和1覆盖范围的多边形，属性为0的多边形即满足形变要求的范围，属性为1的多边形即不满足形变要求的范围；

步骤S11：多个分区范围的多边形叠置分析计算

步骤S11.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min，H_max]共计算n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)；

步骤S11.2读取需要考虑投影形变的多个分区范围面，记为subA_j(j＝1，2，...，k)；

步骤S11.3遍历各分区范围面subA_j，并分别对n个投影面的变形多边形范围prj_i进行空间叠置分析，叠置分析采用相交的策略，将prj_i与subA_j的相交面积记为Intersection_i，并计算相交面同对应分区的范围面的面积占比，取占比最大的投影面高程作为该行政区范围的投影面；

步骤S11.4计算Intersection_i同各分区的面积占比，统计各投影面能达到的形变要求比例；

步骤S12：一个分区范围的多边形叠置分析计算

步骤S12.1对不同投影面计算的矢量多边形进行分析，假设投影面的区域范围[H_min，H_max]共计算了n个投影面，各投影面属性为0的多边形范围记为prj_i(i＝1，2，...，n)；

步骤S12.2计算各投影面属性为0的多边形面积占整个投影区域的比例，以占比最大的矢量多边形所在的投影面，作为主投影面，记为prj_max；

步骤S12.3遍历各投影面prj_i(i＝1，2，...，n，i≠max)，分别同主投影面的多边形prj_max进行空间叠置分析，叠置分析采用删除策略，即prj_i面删除面prj_max，剩余的面则为第i个投影面排除主投影面后符合形变要求的面积，记为Erase_i；

步骤S12.4对删除后的面Erase_i按照面积从大到小进行排序，面积最大的多边形确定为次投影面prj_max-1；

步骤S12.5重复步骤S12.3和步骤S12.4，其中遍历各投影面prj_i(i＝1，2，...，n，i≠max，max-1)，prj_max-1替代prj_max作为主投影面，并以此类推，分别确定出后续的投影面prj_max-2，prj_max-3...prj₁；

步骤S12.6根据叠置分析计算的投影面prj_max-j(j＝0，1，...，i-1)的面积，分别计算所占整个分区面积的比例大小，确定整个分区投影面的个数。

2.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S1中所述纬度B和经度L的计算过程为：通过DEM数据的左上角经纬度及格网分辨率大小，计算每一格网的纬度B和经度L，即可得到大地坐标(B,L)；如果DEM数据的参考椭球不是城市坐标系统参考椭球，采用布尔沙模型将计算的纬度B和经度L转换为对应参考椭球的大地坐标(B,L)。

3.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S1中所述DEM数据为一个二维数组，行号代表纬度的索引，列号代表经度的索引，数组内的取值为该格网的平均海拔高程值h。

4.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S2中所述高程异常ε求解过程为：如果有可利用的区域似大地水准面模型，通过大地坐标(B,L)，查询该点对应的高程异常ε值；如果没有可利用的区域似大地水准面模型，则利用全球超高阶地球重力场模型，通过大地坐标(B,L)，计算对应的高程异常ε值。

5.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S2中所述大地高H计算见公式(1)：

H＝h+ε (1)

式中，H为大地高，h为海拔高程值，ε为高程异常。

6.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S3.1中所述中央子午线L₀的计算见公式(2)，所述所在统一3°带的中央子午线L'₀,见公式(3)：

L₀＝(minL+maxL)/2 (2)

L'₀＝(int)(L₀/3)*3 (3)

7.根据权利要求1所述一种基于DEM的城市平面坐标系统建立方法，其特征在于，步骤S3.3中所述卯酉圈曲率半径N计算见公式(4)：

式中，a为椭球长半轴，e为椭球第一偏心率。

Images