独立坐标系的建立

# 只需 shift+回车 运行本单元格，就可以让jupyter notebook宽屏显示
from IPython.core.display import display, HTML
display(HTML('<style>.container { width:100% !important; }</style>'))

一、地方坐标系简介

地方坐标系（Local Coordinate System）是因建设、城市规划和科学研究需要而在局部地区建立的相对独立的坐标系统，是在局部地区建立平面控制网时，根据需要投影到任意选定面上和（或）采用地方子午线为中央子午线的一种直角坐标系，属于参心坐标系。地方独立坐标系通常采用的是高斯克吕格正形投影平面直角坐标系，然后把独立测量的工程控制网建立在当地海拔高程面，并与当地子午线为中央子午线作为投影变换的依据。

二、建立背景

• 进行工程测量建立平面控制网时，如局部地区没有已知控制点可利用，则选择网中某一点假定其坐标，选定某一边假定其坐标方位角，以此为起算数据推算网中各点的坐标。
• 在国家控制网未扩展到的地区，为了测绘地形图而布设控制网时，可在网中选一点观测其天文经纬度和至另一点的天文方位角，按统一投影带的划分，将该点的天文经纬度换算为平面直角坐标，将天文方位角换算为坐标方位角，以此为起算数据，推算网中各点的坐标 。
• 独立坐标系主要是根据城市或工程建设需求而建立的，其主要特点是限制长度变形，要求实地量测边长与坐标反算边长应满足2.5厘米/公里限差。然而，采用国家坐标系统在高海拔地区或离中央子午线较远地方不能满足这一要求，这就要考虑建立地方独立坐标系。

三、建立方法和步骤

3.1 建立方法概述

• 把中央子午线移到测区中央，归化高程面提高到该测区的平均高程面上，建立任意带高斯正形投影平面直角坐标系，这样可以使测区的长度变形在测区中央几乎为零。当测区高差起伏在100米范围内时可以保证离中央子午线40千米以内的地区其长度变形在每公里2.5厘米以内（可控制的东西宽度100千米）。这种地方独立坐标系最适合工程建设地区的需要，因此，在工程建设区域面积不是太大、东西跨度在80千米可以完全满足需要。
• 采用抵偿高程面的方法建立独立坐标系，即中央子午线保持不变，选择某一高程面作为归化高程面，使高程归化改正和高斯投影变形改正相互抵消，使测区中央的两项投影变形接近于零。
• 以上两种方法建立独立坐标系都变动了高程归化面，这将产生一个新椭球，这不仅要计算出新椭球参数，还要把本地区国家坐标系控制点转换到新产生的椭球面上作为独立坐标系的起算点，计算较复杂。为了避免这些复杂的计算，可以采用不变动高程归化面(长度仍然归化到国家坐标系参考椭球面)，只移动中央子午线的办法来建立独立坐标系。
• 利用GPS-RTK技术建立独立坐标系。近年来，随着测绘技术的发展，GPS相对定位技术已经成为改造和建立城市坐标系和控制网的主要技术手段。在常规测量中，这种独立坐标系只是一种高斯平面直角坐标系，而在采用GPS-RTK采集数据时，独立坐标系就是一种不同于国家坐标系的参心坐标系。

  建立地方独立坐标系的常规方法是以一个国家大地控制点和一条边的方位角作为起算数据，观测边长投影到某特定面（测区平均高程面、抵偿面）上。建立的方式的共同点都有自己的原点，自己的定向，即控制网便是作为独立坐标系建立的参考。

3.2 建立考虑的因素

  在独立坐标系的建立过程中，应主要注意中央子午线、投影面和坐标系参考椭球三个方面。

• 中央子午线： 独立坐标系的中央子午线既可与国家坐标系标准带的中央子午线重合，也可与其不重合。但当测区离标准带中央子午线较远时，可选取过测区中心点或过某点的经线作为独立坐标的新中央子午线。
• 投影面： 在实际工程中，若变换中央子午线还是不能有效解决投影变形的问题，就要考虑建立合适的投影面参数。一般情况下可选择研究区的平均高程面或者抵偿高程面作为独立坐标的投影面。
• 参考椭球： 地方独立坐标系参考椭球相应的参数设置在原则上应使得椭球面与投影面的拟合程度达到最好，投影长度的变形值降低到最小，同时也要满足方便与我国国家坐标系统进行相互的坐标换算。

3.3 基于坐标椭球参数建立独立坐标系

  一般情况下，独立坐标系采用国家坐标系椭球参数，(如基于2000国家大地坐标系建立的独立坐标系统，称为2000独立坐标系)，根据城市或区域中心的地理位置设定高斯投影中央子午线，或以测区平均高程面作为坐标投影面，通过抬高或降低坐标投影面的方法解决长度变形问题；有些独立坐标进行加常数或者平移旋转变换等。以重庆市为例，我们知道重庆位于东经105°17′至110°11′、北纬28°10′至32°13′之间，重庆市城区的经纬度北纬29.35°东经106.33°。最常用的就是根据具体地理位置设定105°、108°或111°为高斯投影中央子午线。为了使横坐标y不出现负值，则无论3°或6°带，每带的纵坐标轴要西移500km，即在每带的东坐标上加500 km。为了指明该点属于何带，还规定在横坐标y值之前，要写上带号。未加500km和带号的横坐标值称为自然值，加上500km和带号的横坐标值称为通用值。

3.4 坐标转换

同椭球系统内坐标转换分为两种情况，一种是同椭球之间的转换，另外则是不同椭球之间的坐标转换。同椭球一种分为空间大地坐标系与空间直角坐标系间的转换，即$(B,L,H)$与$(X,Y,Z)$之间的转换。具体公式如下：

$(B,L,H)\Rightarrow (X,Y,Z)$：
 
$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}(N+H)cosBcosL\\ (N+H)cosBsinL\\ \left[N(1-e^2)+H\right]sinB\end{array}\right]\text{\hspace{1em}(3.4.1)}$$
 
$(X,Y,Z)\Rightarrow (B,L,H)$：
 
$$\left[\begin{array}{c}L\\ B\\ H\end{array}\right]=\left[\begin{array}{c}\arctan \frac{Y}{X}\\ \arctan \frac{1}{\sqrt{X^2+Y^2}}\left(Z+\frac{a{e}^2\tan B}{\sqrt{1+{\tan }^{2}B-e^2{\tan }^{2}B}}\right)\\ \frac{\sqrt{X^2+Y^2}}{\cos B}-N\end{array}\right]\text{\hspace{1em}(3.4.2)}$$
 
式$(3.4.1)$和$(3.4.2)$的公式相关符号说明：

符号	说明及计算方法
$X,Y,Z$	直角坐标系下的坐标
$B,L,H$	大地坐标系下的坐标，分别表示大地纬度、大地经度和大地高
$a$	参考椭球长半轴
$b$	参考椭球短半轴
$f$	扁率的计算：$f=\frac{a-b}{a}$
$e$	第一偏心率的计算：$e=\sqrt{\frac{a^2-b^2}{a^2}}$
$e^{\prime }$	第二偏心率的计算：$e^{\prime }=\sqrt{\frac{a^2-b^2}{a^2}}$
$W$	辅助函数$W=\sqrt{1-e^{2}si{n}^{2}B}$
$V$	辅助函数$V=\sqrt{1+e^{\prime 2}co{s}^{2}B}$
$N$	卯酉圈曲率半径$N=\frac{a}{W}$
$L_0$	中央子午线经度
$l^{\prime \prime }$	$l^{\prime \prime }=L-L_0$
$t$	$t=\tan B$
${\eta }^2$	${\eta }^2=e^{\prime 2}{\cos }^{2}B$
${\rho }^{\prime \prime }$	$\rho ^{\prime \prime}=\frac{180}{\pi}*3600 \approx 206265 $

import math 
180 * 3600 / math.pi

206264.80624709636

另一种为空间大地坐标系与高斯平面直角坐标系之间的转换，采用高斯正反算的方法进行计算。高斯正算是指将大地经度和大地纬度换算为高斯平面坐标的计算；高斯反算是指将高斯平面坐标换算为大地经度和大地纬度的计算。其换算公式的主要参数有中央子午线、投影高、投影纬度、东北平移量、投影比例。高斯正反算公式如下：

子午线弧长计算公式：

$$\begin{array}{rl}X& =a_{0}B-\sin B\cos B\left[(a_2-a_4+a_6)+(2a_4-\frac{16}{3}{a}_6)si{n}^{2}B+\frac{16}{3}{a}_{6}si{n}^{4}B\right]\\ & =a_{0}B-\frac{a_2}{2}sin2B+\frac{a_4}{4}sin4B-\frac{a_6}{6}sin6B+\frac{a_8}{8}sin8B\end{array}$$

式中：$a_0,a_2,a_4,a_6,a_8$为基本常量，计算公式如下：

$$\begin{array}{rl}{a}_0& =m_0+\frac{1}{2}{m}_2+\frac{3}{8}{m}_4+\frac{5}{16}{m}_6+\frac{35}{128}{m}_8\\ a_2& =\frac{1}{2}{m}_2+\frac{1}{2}{m}_4+\frac{15}{32}{m}_6+\frac{7}{16}{m}_8\\ a_4& =\frac{1}{8}{m}_4+\frac{3}{16}{m}_6+\frac{7}{32}{m}_8\\ a_6& =\frac{1}{32}{m}_6+\frac{1}{16}{m}_8\\ a_8& =\frac{1}{128}{m}_8\end{array}$$
式中：$m_0,m_2,m_4,m_6,m_8$为基本常量，计算公式为：

$$m_o=a(1-e^2)；m_2=\frac{3}{2}{e}^2{m}_0；m_4=5e^2{m}_2；m_6=\frac{7}{6}{e}^2{m}_4；m_8=\frac{9}{8}{e}^2{m}_6$$
 
高斯正算$(B,L)\Rightarrow (X,Y)$：
 
$$\begin{array}{rl}x& =X+\frac{N}{2{\rho }^{\prime \prime 2}}\sin B\cos B\cdot l^{\prime \prime 2}+\frac{N}{24{\rho }^{\prime \prime 4}}\sin B{\cos }^{3}B(5-t^2+9{\eta }^2+4{\eta }^4)\cdot l^{\prime \prime 4}+\frac{N}{720{\rho }^{\prime \prime 6}}\sin B{\cos }^{5}B(61-58t^2+t^4)\cdot l^{\prime \prime 6}\\ y& =\frac{N}{{\rho }^{\prime \prime }}\cos B\cdot l^{\prime \prime }+\frac{N}{6{\rho }^{\prime \prime 3}}{\cos }^{3}B(1-t^2+{\eta }^2)\cdot l^{\prime \prime 3}+\frac{N}{120{\rho }^{\prime \prime 5}}{\cos }^{5}B(5-18t^2+t^4+14{\eta }^2-58{\eta }^2{t}^2)\cdot l^{\prime \prime 5}\end{array}$$
 
高斯反算 $(X,Y)\Rightarrow (B,L)$：
 
$$\begin{array}{rl}B& =B_f-\frac{t_f}{2M_f{N}_f}{y}^2+\frac{t_f}{24M_f{N}_f^3}(5+3t_f^2+{\eta }_f^2-9{\eta }_f^2{t}_f^2)y^4-\frac{t_f}{720M_f{N}_f^5}y(61+90t_f^2+45t_f^4)y^6\\ l& =\frac{y}{N_{f}cos{B}_f}-\frac{y^3}{6N_f^{5}cos{B}_f}(1+2t_f^2+{\eta }_f^2)+\frac{y^5}{120N_f^{5}cos{B}_f}(5+28t_f^2+24t_f^4+6{\eta }_f^2+8{\eta }_f^2{t}_f^2)\\ L& =l+L_0\end{array}$$
 
$B_f$为底点维度，即是当$x=X$时的子午线弧长所对应的纬度，按照子午线弧长$X$公式迭代进行计算。
开始时设$B_f^1=\frac{X}{a_0}$，然后按照下式迭代计算：
 
$$\begin{array}{rl}{B}_f^{i+1}& =\frac{X(X-F(B_f^i))}{a_0}\\ F(B_f^i)& =-\frac{a_2}{2}sin2B_f^i+\frac{a_4}{4}sin4B_f^i-\frac{a_6}{6}sin6B_f^i+\frac{a_8}{8}sin8B_f^i\end{array}$$

重复迭代至$B^{i+1}_f - B^i_f < \epsilon $为止。
 
其余参数为：

$$\begin{array}{rl}{N}_f& =a(1-e^{2}si{n}^2{B}_f{)}^{-\frac{1}{2}}\\ M_f& =a(1-e^2)(1-e^{2}si{n}^2{B}_f{)}^{-\frac{3}{2}}\\ t_f& =tan{B}_f\\ {\eta }_f^2& =e^{\prime 2}co{s}^2{B}_f\end{array}$$
 
不同椭球间的坐标转换的方法是：七参数（空间直角坐标系统间转换）、四参数+高程拟合参数（平面直角坐标系统间转换）。常见的Bursa七参数坐标转换公式如下：

$${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_C=\left[\begin{array}{ccccccc}1& 0& 0& 0& -Z_D& Y_D& X_D\\ 0& 1& 0& Z_D& 0& -X_D& Y_D\\ 0& 0& 1& -Y_D& X_D& 0& Z_D\end{array}\right]\left[\begin{array}{c}{T}_X\\ T_Y\\ T_Z\\ {\epsilon }_X\\ {\epsilon }_Y\\ {\epsilon }_Z\\ m\end{array}\right]+{\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_D\text{\hspace{1em}(1)}$$

四、地方坐标系建立意义与弊端

4.1 地方坐标系建立意义

  在实际测量作业中，我们通常依据不同的用途和工程项目，采用不同的坐标系来满足工程项目的需要。高斯一克吕格投影分带有效的限制了长度变形，但是在投影带的边缘地区，其长度变形仍然达到了很大的数值 。为了达到城市和工程建设的要求，我们就必须对长度变形加以限制，为此考虑建立独立坐标系， 目的是减小高程归化与投影长度变形产生的影响，将它们控制在一个微小的范围，使计算的长度在实际应用时(如工程放样时)不需要做任何的改正。

4.2 地方坐标系建立弊端

• 起算点坐标从国家坐标的参考椭球高斯成果直接搬至地方独立坐标系的投影面，这在理论上不严密，同时因起算点不同，整个网成果不同；
• 与国家大地控制点不能严格转换，不利于资源共享；
• 不能充分利用国家大地控制点提高网的精度，对于带状控制网（公路、输电线路等）尤为突出。由此，应该建立一种既与国家坐标系有严密换算公式，又能保证投影变形在规定范围的地方独立坐标系统。

  在城市范围内布设控制网时，应考虑不仅要满足大比例尺测图的需要，还要满足一般工程放样的需要，通常情况下要求控制网由平面直角坐标反算的长度与实测的长度尽可能地相符，而国家坐标系的坐标成果则往往无法满足这些要求，这是因为国家坐标系每个投影带都是按照一定的间隔划分，由西向东有规律地分布，其中央子午线不可能恰好落在每个城市的中央。为了减小长度投影变形所产生的影响，使由控制点的平面直角坐标反算出来的长度在实际利用时不需要做任何改正，方便测绘实际作业，根据《城市测量规范》的要求，需要建立有别与国家统一坐标系统的城市独立坐标系统。