一种基于地理信息应用城市独立坐标系参数获取方法

技术领域

本发明涉及利用高斯投影技术，通过估算城市独立坐标系参数，从而实现地理信息在各坐标系间相互转换的方法，具体是一种基于地理信息应用城市独立坐标系参数获取方法。

背景技术

随着计算机及信息技术的进步，地理信息产业取得了长足的发展。作为地理信息产业核心的资产──地理信息数据的重要性也日益凸显。我国现有的地理空间信息横向上分布在不同委办局，纵向上又分散在市、区（县）、乡（镇）等多级行政级别，呈现出典型的“多重分布”的复杂特性，在坐标体系上，我国先后建成了北京54坐标系、西安80坐标系、2000国家坐标系，国际上较通用的有WGS84坐标系，而在城市规划、建设、管理等实际工作中，又多采用城市独立坐标系，导致在地理信息应用过程中，需要通过坐标转换来整合不同空间数据。

目前，坐标转换理论很成熟，相关的方法研究也比较多，最常用的有四参数、七参数。其中四参数是平面仿射变换，比较简单，缺点是没有考虑椭球体的变形，不能进行大范围高精度的坐标转换。而七参数转换基于参考椭球，理论严谨，转换精度高，控制范围大，是地理信息应用中常用的一种方法。但七参数法需要知道各城市独立坐标系的参数，但这些数据涉密，在实际应用中很难获取。

建立城市独立坐标系的目的是为了满足投影变形小于2.5cm/km的要求。我国大地测量法和有关测量规范中明确规定，国家大地测量控制网依高斯投影方法按3°带或6°带进行分带和计算，并把观测成果归算到参考椭球面上。这样的规定，不但符合高斯投影的分带原则和计算方法，与国际惯例相一致，而且也便于大地测量成果的统一、使用和互算。但由于我国国土辽阔，地形复杂，我国的北京54坐标系、西安80坐标系的参考椭球长半径比国际椭球半径大，与地球表面存在一定的套合差，所以国家在对原有北京54坐标系、西安80坐标系进行优化的基础上，建立2000国家坐标系。同时根据高斯投影变化规律，离中央子午线越远，变形越大，按照国家坐标分带不能满足变形小于2.5cm/km的要求，因此我国大部分城市和地区布设城市测量控制网都基于北京54坐标系、西安80坐标系或2000国家坐标系参考椭球，采用任意带的独立高斯投影平面直角坐标系，归算测量结果的参考面也是根据需要选定的方法来建立城市独立坐标系。

在计算城市独立坐标系参数中，通常利用高斯投影理论，进行大地坐标和高斯平面直角坐标间的相互转换，其中城市独立坐标系采用高斯平面直角坐标，北京54坐标系、西安80坐标系或2000国家坐标系采用大地坐标，方法如下：

1、大地坐标（B，L，H）：大地坐标是采用大地纬度B、经度L和大地高程H来描述空间位置的。纬度是空间的点与参考椭球面的法线与赤道面的夹角；经度是空间的点与参考椭球的自转轴所在的面与参考椭球的起始子午面的夹角λ；大地高程是空间的点沿着参考椭球的法线方向到参考椭球面的距离，参见图1。

2、高斯平面直角坐标（x，y，H）：其中，x，y表示平面坐标，H表示大地高程。椭球体中心O在椭圆柱中心轴上，椭球体南北极与椭圆柱相切，并使某一子午线与椭圆柱相切。此子午线称中央子午线。然后将椭球体面上的点、线按正形投影条件投影到椭圆柱上，再沿椭圆柱N、S点母线割开，并展成平面，即成为高斯投影平面，参见图2。在此平面上：

①中央子午线是直线，其长度不变形，离开中央子午线的其他子午线是弧形，凹向中央子午线。离开中央子午线越远，变形越大。

②投影后赤道是一条直线，赤道与中央子午线保持正交。

③离开赤道的纬线是弧线，凸向赤道。

3、空间直角坐标（X，Y，Z）的坐标原点位于参考椭球的中心，Z轴指向参考椭球的北极，X轴指向起始子午面与赤道的交点，Y轴位于赤道面上切按右手系于X轴呈90度夹角，某点中的坐标可用该点在此坐标系的各个坐标轴上的投影来表示。参见图1。

大地坐标向高斯平面坐标转换（高斯投影正算公式）：

$>x=X+\mathrm{Nt}{\cos }^{2}B\frac{1^2}{{\rho }^2}[0.5+\frac{1}{24}(5-t^2+9{\eta }^2+4{\eta }^4)\mathrm{co}{s}^{2}B\frac{1^2}{{\rho }^2}+\frac{1}{720}(61-58t^2+t^4){\cos }^{4}B\frac{1^4}{{\rho }^4}]---\left(1\right)$>

$>y=N\cos B\frac{1}{\rho }[1+\frac{1}{6}(1-t^2+{\eta }^2)\mathrm{co}{s}^{2}B\frac{1^2}{{\rho }^2}+\frac{1}{120}(5-18t^2+t^4+14{\eta }^2-58{\eta }^2{t}^2)N{\cos }^{4}B\frac{1^4}{{\rho }^4}]---\left(2\right)$>

X=a(1-e²)(A′arcB-B′sin2B+C′sin4B-D′sin6B+E′sin8B-F′sin10B+G′sin12B)

$>A^{\prime }=1+\frac{3}{4}{e}^2+\frac{45}{64}{e}^4+\frac{175}{256}{e}^6+\frac{11025}{16384}{e}^8+\frac{43659}{65536}{e}^{10}+\frac{693693}{1048576}{e}^{12}$>

$>B^{\prime }=\frac{3}{8}{e}^2+\frac{15}{32}{e}^4+\frac{525}{1024}{e}^6+\frac{2205}{4096}{e}^8+\frac{72765}{131072}{e}^{10}+\frac{297297}{524288}{e}^{12}$>

$>C^{\prime }=\frac{15}{256}{e}^4+\frac{105}{1024}{e}^6+\frac{2205}{16384}{e}^8+\frac{10395}{65536}{e}^{10}+\frac{1486485}{8388608}{e}^{12}$>

$>D^{\prime }=\frac{35}{3072}{e}^6+\frac{105}{4096}{e}^8+\frac{10395}{262144}{e}^{10}+\frac{55055}{1048576}{e}^{12}$>

$>E^{\prime }=\frac{315}{131072}{e}^8+\frac{3465}{524288}{e}^{10}+\frac{99099}{8388608}{e}^{12}$>

$>F^{\prime }=\frac{693}{1310720}{e}^{10}+\frac{9009}{5242880}{e}^{12}$>

$>G^{\prime }=\frac{1001}{8388608}{e}^{12}$>

其中：

ρ=180×3600/π

$>f=\frac{a-b}{a}$>

$>b=a\sqrt{1-e^2}$>

e─第一偏心率

$>e=\frac{\sqrt{a^2-b^2}}{a}$>e²=2f-f²

e′─第二偏心率

$>e^{\prime }=\frac{\sqrt{a^2-b^2}}{b}$>

η²=e'²cos²B

t=tgB

$>V=\sqrt{1+e^{\prime 2}{\cos }^{2}B}$>V²=1+η²

$>W=\sqrt{1-e^2{\sin }^{2}B}$>

B为纬度，单位弧度

$>c=\frac{a^2}{b}$>

M─子午圈曲率半径$>M=\frac{a(1-e^2)}{W^3}=\frac{c}{V^3}$>

N─卯酉圈曲率半径$>N=\frac{a}{W}=\frac{c}{V}$>

在进行城市独立坐标系与国家通用坐标系进行转换时，必须要知道城市独立坐标系的参数，但城市独立坐标系的参数涉密，在实际应用中很难获取。因此要实现城市独立坐标系与国家通用坐标系之间的相互转换，亟需一种能获得城市独立坐标系的参数的方法。

发明内容

本发明的目的在于利用高斯投影技术，计算城市独立坐标系参数，包括参考椭球扁率，参考椭球长半径，中央子午线，北方向偏移量和东方向偏移量，在避免城市独立坐标系参数的公布的前提下，而实现城市独立坐标系与北京54坐标系、西安80坐标系或2000国家坐标系成果之间的相互转换。

通常城市独立坐标系的参数包括参考椭球扁率e，参考椭球长半径a，中央子午线L₀，北方向偏移量X₀和东方向偏移量Y₀，其中参考椭球扁率e已知，与国家通用坐标系参考椭球保持一致，其他几个参数未知。所以确定一个城市独立坐标系参数，就是要计算出参考椭球长半径a，中央子午线L₀，北方向偏移量X₀和东方向偏移量Y₀四个参数。已知某城市独立坐标系和国家通用坐标系两种坐标系下的同名点成果各三个，通过对高斯投影正算公式进行变换，得到一个关于中央子午线L₀的一元高次方程，可以通过迭代算法解算出L₀，进而可以计算出椭球长半径，北方向偏移量和东方向偏移量，从而获取城市独立坐标系参数。因而本发明的技术方案是：

1）利用GPS联测某城市独立坐标系四等以上的控制点（为了提高坐标计算的精度，所以同名点成果的测量等级尽量高），获得至少三个该城市独立坐标系同名点成果P₁、P₂、P₃，在与测得所述城市独立坐标系同名点成果相同的三点上利用GPS与国家通用坐标系联测四等以上控制点，得到国家通用坐标系同名点成果P₁′、P₂′、P₃′；

2）记

$>F\left(l\right)=\frac{1}{2}t{\cos }^2{\mathrm{Bl}}^2+\frac{1}{24}t(5-t^2+9{\eta }^2+4{\eta }^4){\cos }^{4}B{l}^4+\frac{1}{720}t(61-58t^2+t^4){\cos }^6{\mathrm{Bl}}^6$>

高斯投影正算公式简化为：

x=X+NF(l)         （3）

$>F_1\left(l\right)=\cos \mathrm{Bl}+\frac{1}{6}(1-t^2+{\eta }^2){\cos }^3{\mathrm{Bl}}^3+\frac{1}{120}(5-18t^2+t^4+14{\eta }^2-58{\eta }^2{t}^2){\cos }^5{\mathrm{Bl}}^5$>

高斯投影正算公式简化为：

y=NF₁(l)         （4）

其中，经差l=L-L₀，以弧度为单位，L₀为中央子午线，L为国家通用坐标系同名点成果的经度，B为国家通用坐标系同名点成果的纬度，x，y为城市独立坐标系同名点成果的平面坐标。

3）由于城市独立坐标系参考椭球通常基于国家通用坐标系参考椭球，采用膨胀法建设，所以城市独立坐标系参考椭球的扁率e总是与对应国家通用坐标系参考椭球的扁率相同，但是二者参考椭球的长半径a不同。这这种情况下城市独立坐标系的经度变化△L=0，纬度变化△B≈0，将P₁和P₁′、P₂和P₂′、P₃和P₃′任意两对同名点成果分别带入式（3）、式（4）整理可得：

记任意两对同名点成果为：国家通用坐标系的记为（B₁，L₁，H₁）和（B₂，L₂，H₂），城市独立坐标系的记为（x₁，y₁，H₁）和（x₂，y₂，H₂）。

x₁=X₀+X₁+N₁F(L₁-L₀)       （5）

y₁=Y₀+N₁F₁(L₁-L₀)           （6）

x₂=X₀+X₂+N₂F(L₂-L₀)       （7）

y₂=Y₀+N₂F₁(L₂-L₀)        （9）

x₁-x₂=X₁-X₂+N₁F(L₁-L₀)-N₂F(L₂-L₀)    （10）

$>\frac{x_1-x_2}{y_1-y_2}=\frac{X_1-X_2+N_{1}F(L_1-L_0)-N_{2}F(L_2-L_0)}{N_1{F}_1(L_1-L_0)-N_2{F}_1(L_2-L_0)}---\left(11\right)$>

式（11）是一个关于L₀的一元高次方程，利用高斯迭代算法可以解出L₀。

将L₀代入式（10）可以求解出参数椭球长半径a

将a，L₀分别代入式（5）、式（6）求解出北方向偏移量X₀和东方向偏移量Y₀。

这样便求得了城市独立坐标系的各个参数。

4）国家通用坐标系参数是已知的，可查阅的，根据步骤3）中获得的城市独立坐标系参数（中央子午线L₀、椭球长半径a、北方偏移量X₀和东方偏移量Y₀），利用P₁、P₂、P₃和P₁′、P₂′、P₃′的两套坐标系同名点成果，基于Bursa七参数转换模型，计算出坐标转换参数，转换参数包括3个平移参数，3个旋转参数和一个比例因子，从而实现坐标转换。

Bursa七参数转换模型：

$>\left(\begin{array}{c}{X}_T\\ Y_T\\ Z_T\end{array}\right)=\left(\begin{array}{c}\mathrm{\Delta X}\\ \mathrm{\Delta Y}\\ \mathrm{\Delta Z}\end{array}\right)+\left(\begin{array}{ccc}0& -Z_S& Y_S\\ Z_S& 0& -X_S\\ -Y_S& X_S& 0\end{array}\right)\left(\begin{array}{c}{ϵ}_X\\ {ϵ}_Y\\ {ϵ}_Z\end{array}\right)+m\begin{array}{}\end{array}\left(\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right)+\left(\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right)$>

式中，[△X △Y △Z]^T为3个平移参数，[ε_X ε_Z ε_Z]^T为3个旋转参数和m为比例因子，[X_S Y_S Z_S]^T，[X_T Y_T Z_T]^T分别为城市独立坐标系同名点成果和国家通用坐标系同名点成果。

具体地，上述国家通用坐标系包括北京54坐标系、西安80坐标系和2000国家坐标系。国家通用坐标系用其他已知坐标系替换，也能使用本发明的方法。

本发明的有益效果是：利用高斯投影技术，只需城市独立坐标系和国家通用坐标系三个同名点成果，就可以近似计算出城市独立坐标系参数，实现两个坐标系间成果的高精度转换，完全满足地理信息坐标转换的需要，避免了独立坐标系参数的公布。本发明可广泛用于不同坐标系的地理空间数据整合和坐标换算，以及在数字城市和智慧城市建设中，进行多源数据空间套合和分析等。

附图说明

图1为大地坐标系示意图；

图2为高斯投影示意图。

具体实施方式

由于城市独立坐标系的参数是涉密的，很难获取，所以不能将本发明方法获得的城市独立坐标系参数与实际的城市独立坐标系参数进行比较，来验证方法的有效性。因此将用本方法获得的城市独立坐标系参数，利用七参数转换模型，将国家通用坐标系的同名点成果转换为独立坐标系转换成果，然后将独立坐标系转换成果与实际测得的独立坐标系同名点成果进行比较，来验证本发明方法的有效性。具体参见以下实施例。

实施例1，某城市独立坐标系与北京54坐标系的转换

首先，利用GPS分别与北京54坐标与某城市独立坐标联测，测得同名点成果8个，如下表1所示。

表1某城市独立坐标系与北京54坐标系同名点成果表

其次，利用1、2点成果根据式11，式10，式5及式6，计算出的坐标参数。

其中：

中央子午线：108:45:05.029260

椭球长半径：6379489.37490557

北方偏移量：-6.15235514800548米

东方偏移量：500144.292887685米

第三，利用该坐标系参数，将1、2、3点成果带入七参数转换模型，计算出七参数为：

X偏移：-5.90702772431541

Y偏移：77.6861151843914

Z偏移：50.5046847679187

X旋转：000:00:00.361474

Y旋转：-000:00:00.356650

Z旋转：-000:00:00.976716

比例因子：0.000180875238540795

第四，利用该七参数进行坐标转换后，比较转换结果与测量值，最大差值为0.036米，完全满足地理信息坐标转换的需要，见表2。

表2某城市独立坐标系与北京54坐标系坐标转换成果较差统计表

实施例2，某城市独立坐标系与西安80坐标系的转换

首先，利用GPS分别与西安80坐标系与某城市独立坐标联测，测得同名点成果5个，如下表3所示。

表3某城市独立坐标系与西安80坐标系的转换同名点成果表

其次，利用1、2点成果根据式11，式10，式5及式6，计算出的坐标参数。

其中：

中央子午线：105:53:58.943278

椭球长半径：6378423.47268279

北方偏移量：-3199971.81273551米

东方偏移量：49971.3935643691米

第三，利用该坐标系参数，将1、2、3点成果带入七参数转换模型，计算出七参数为：

X偏移：-2.75906461090835

Y偏移：-7.97653608962355

Z偏移：-7.82526087134284

X旋转：000:00:00.086500

Y旋转：-000:00:00.710559

Z旋转：-000:00:00.215598

比例因子：0.0000459872988582905

第四，利用该七参数进行坐标转换后，转换结果与测量值进行比较，最大差值为0.0272米，完全满足地理信息坐标转换的需要，见表4。

表4某城市独立坐标系与西安80坐标系的转换坐标转换成果较差统计表

实施例3，某城市独立坐标系与2000国家坐标系的转换

首先，利用GPS分别与2000国家坐标与某城市独立坐标联测，测得同名点成果9个，如下表5所示。

表5独立坐标系与2000国家坐标系同名点成果表

其次，利用1、2点成果根据式11，式10，式5及式6，计算出的坐标参数。

其中：

中央子午线：108:22:17.492463

椭球长半径：6378421.52585621

北方偏移量：-1.9293381461545米

东方偏移量：499965.255096663米

第三，利用该坐标系参数，将1、2、3点成果带入七参数转换模型，计算出七参数为：

X偏移：-0.0144736003057915

Y偏移：1.45046010742954

Z偏移：3.63460416544149

X旋转：000:00:00.087127

Y旋转：-000:00:00.012301

Z旋转：-000:00:00.039853

比例因子：0.0000441302554666525

第四，利用该七参数进行坐标转换后，转换结果与测量结果进行比较，最大差值为0.042米，完全满足地理信息坐标转换的需要，见表6。

表6某城市独立坐标系与2000国家坐标系坐标转换成果较差统计表

本发明利用高斯投影理论，在已知至少三个某城市独立坐标系与国家通用坐标系的同名点，未知该独立坐标系参数的情况下，用于计算该城市的近似坐标系参数，实现地理信息成果的高精度转换，避免了城市独立坐标系参数的公布，促进了地理信息系统的应用与推广。