地理坐标到本地CAD坐标的坐标转换公式20140 209

地理全局坐标到本地局部CAD坐标的坐标转换公式

1、卫星位置用地球地心惯性坐标系ECI，由赤道面X轴（春分，由地心指向太阳），赤道面y轴（由地心指向春分向东90度方向），北极z轴（由地心指向地球北极），构成右手迪卡尔坐标系；

2、地面定位目标位置在导航电文中，使用WGS84地心地固坐标系ECEF，由赤道面X轴（由地心指向本初子午线，即0度经线方向），赤道面y轴（由地心指向东经90度方向），北极z轴（由地心指向地球北极），构成右手迪卡尔坐标系；

3、本地水平坐标系LTP（ENU），由本地切平面X，轴（由切点指向东，即东经方向），切平面y，轴（由切点指向北，即北纬方向），高度z，轴（由切点指向海拔高度方向），构成右手迪卡尔坐标系；

4、地固坐标系ECEF映射到经纬度坐标系的就是地理坐标系，由东经lon轴（由切点指向东经方向），北纬lat轴（由切点指向北纬方向）；

5、在由本地切平面不大时，切平面的EN坐标系和地理坐标系存在近似线性变换关系，使用线性坐标转换矩阵可满足工程应用精度；

6、在坐标转换精度要求不高时，且坐标转换范围不大时，可忽略地球扁率的影响，用球体近似表达地球表面；

7、EN（x，y）坐标系和地理坐标系（lon，lat）的近似坐标变换公式如下：

∆x=r*∆lon，其中，x的单位是m，lat的单位是弧度；r=R*COS(lat)，单位是m；

∆Y=R*∆lat，其中，y的单位是m，lat的单位是弧度，R是地球平均半径，6371000m；所以有

{∆x，∆y}T= RC{∆lon，∆lat }T

其中C为二阶对角线矩阵：

c11=cos(lat0)

c22= 1

c12= c21=0

且：{∆lon，∆lat }T= {lon- lon0，lat- lat0 }T，

得到

{∆x，∆y}T = RC{ lon- lon0，lat- lat0 }T

8、EN（∆x，∆y）坐标系到CAD坐标系（X，Y）的坐标变换公式如下：

{X，Y}T= K{∆x，∆y }T+S，这里{∆x，∆y}T = RC{ lon- lon0，lat- lat0 }T，S={ X0，Y0 }T

其中K为二阶旋转矩阵：

k11= cos(α-θ)

k12= -sin(α-θ)

k21= sin(α-θ)

k22= cos(α-θ)

通过基线测量得到两个旋转角：

tanθ=dy/dx= dlat /(dlon* cos(lat0))

tanα=DY/DX

式中

dlon = lon1- lon0

dlat = lat1- lat0

DX = X 1- X 0

DY = Y 1- Y 0

见坐标变换示意图，

9、综合7和8两步变换，得到地理坐标到CAD坐标系（X，Y）的坐标变换公式如下：

{X，Y}T= KRC{ lon- lon0，lat- lat0 }T +{ X0，Y0 }T，

地理坐标到CAD坐标变换示意图

10、考虑到CAD坐标系还需要归一化到切图坐标系下，对应的坐标变换公式如下：

{x，y}T= (60/L m)* ({ X，Y }T- { X C，Y C } T)，

即：{x，y}T= (60/L m)* (KRC{ lon- lon0，lat- lat0 }T +{ X0，Y0 }T- { X C，Y C } T)

其中：（X C，Y C）为楼面积外接矩形中心坐标，L m为楼面积外接矩形的最大边长，地理坐标单位是弧度。其中：

L m=Max{X max-X min，Y max-Y min};

X C=(X max-X min)/2+ X min;

Y C=(Y max-Y min)/2+ Y min;

最终坐标转换计算公式为：

x= (60/L m)* (6371000* (cos(α-θ) * cos(lat0) * (lon- lon0)- sin(α-θ)* (lat- lat0) )+ X0- X C)

y= (60/L m)* (6371000* (sin(α-θ) * cos(lat0) * (lon- lon0)+ cos(α-θ)* (lat- lat0) )+ Y0- Y C)

式中：L m，X C，Y C由楼面图数据计算得到，lon0，lat0，X0，Y0，lon1，lat1，X1，Y1由基线测量得到，α，θ由基线坐标数据计算得到。上式用于从GPS坐标向直角绘图坐标的转换。

每栋楼的位置特征由以下9个独立参量确定：lon0，lat0，X0，Y0，α，θ，L m，X C，Y C。上式中的地理坐标单位为弧度，若地理坐标以度为单位，则：

x= (60/L m)* (6371000* (cos(α-θ) * cos(lat0) * (lon- lon0)*π/180- sin(α-θ)* (lat- lat0) *π/180 )+ X0- X C)；

y= (60/L m)* (6371000* (sin(α-θ) * cos(lat0) * (lon- lon0) *π/180+ cos(α-θ)* (lat- lat0)

*π/180 )+ Y0- Y C)；

上式用于GPS输出的地理坐标转换为CAD切图坐标，在室内定位终端程序中使用。

11、同理，从CAD坐标系到地理坐标系的坐标变换公式如下：

由：{x，y}T= (60/L m)* ({ X，Y }T- { X C，Y C } T)，

有：(L m /60)*{x，y}T+ { X C，Y C } T = { X，Y }T= KRC{ lon- lon0，lat- lat0 }T +{ X0，Y0 }T，

(L m /60)*{x，y}T+ { X C，Y C } T - { X0，Y0 }T = KRC{ lon- lon0，lat- lat0 }T，

(L m /60)*{x，y}T+ { X C - X0，Y C - Y0} T = KRC{ lon- lon0，lat- lat0 }T，

(L m /60)*{x，y}T+ { X C - X0，Y C - Y0} T = RKC{ lon- lon0，lat- lat0 }T，

1/R* ( (L m /60)*{x，y}T+ { X C - X0，Y C - Y0} T) = KC{ lon- lon0，lat- lat0 }T，

K-1/R*( (L m /60)*{x，y}T+ { X C - X0，Y C - Y0} T) = C{ lon- lon0，lat- lat0 }T，

C-1 K-1/R* ((L m /60)*{x，y}T+ { X C - X0，Y C - Y0} T) = { lon- lon0，lat- lat0 }T，

变形得：{ lon，lat }T= C-1 K-1/R * ((L m /60)*{x，y}T+ { X C-X0，Y C -Y0}T)+ { lon0，lat0 }T，

其中K-1为K的逆矩阵：

k11= cos(α-θ)

k12= sin(α-θ)

k21= -sin(α-θ)

k22= cos(α-θ)

其中C -1为C的逆矩阵：

c11=1/cos(lat0)

c22= 1

c12= c21=0

将上式展开，最终得到从CAD切图坐标系到地理坐标系的坐标转换计算公式为：lon = ( (L m /60*x + X C - X0) cos(α-θ)+ (L m /60*y + Y C - Y0) sin(α-θ) )/6371000/cos(lat0) + lon0

lat = (-(L m /60*x + X C - X0) sin(α-θ) + (L m /60*y + Y C - Y0) cos(α-θ) )/6371000 + lat0

上式计算结果的单位为弧度，若以度为单位，则：

lon = ( (L m /60*x + X C - X0) cos(α-θ)+ (L m /60*y + Y C - Y0) sin(α-θ) )/6371000/cos(lat0)*180/π + lon0 lat = (-(L m /60*x + X C - X0) sin(α-θ) + (L m /60*y + Y C - Y0) cos(α-θ) )/6371000*180/π + lat0

12、准备CAD切图坐标系下的地理底图数据

1）、计算裁剪矩形角点地理坐标：

利用上式，分别代入裁剪四边形角点的CAD切图坐标，计算裁剪四边形角点地理坐标，得到地理坐标裁剪四边形（Excel表1）；

裁剪四边形角点的CAD切图坐标为：

P1={-90,-90};

P2={-90, 90};

P3={ 90, 90};

P4={ 90, -90};

2）、用实测基线在ArcMap中确定高德地图数据偏差：

首先在ArcMap中，根据基线起点、终点的地标特征，在高德地图上确定基线对应位置，再读取基线在高德地图上的起点和终点地理坐标值{lon g0，lat g0}，{lon g1，lat g1}；

确定平移偏差值：Dlon g= lon g0 - lon 0，Dlat g= lat g0 - lat 0，

确定旋转偏差角：为γ-β，

其中