CN109683184A

CN109683184A - 一种倾斜地面坐标测量方法

Info

Publication number: CN109683184A
Application number: CN201811595494.7A
Authority: CN
Inventors: 姚文杰; 李振亚; 刘凡
Original assignee: Harxon Corp
Current assignee: Harxon Corp
Priority date: 2018-12-25
Filing date: 2018-12-25
Publication date: 2019-04-26
Anticipated expiration: 2038-12-25
Also published as: CN109683184B

Abstract

本发明提供一种倾斜地面坐标测量方法，属于地面测绘技术领域。本发明的方法为通过RTK法测量待测点，直接计算待测点的坐标，以待测点的坐标的计算结果作为初值进行间接迭代计算，计算出实测大地纬度B，根据实测大地纬度B计算出实测高程H。本发明的方法的精度明显优于现有技术方法的精度，同时不受大地高H的影响。

Description

一种倾斜地面坐标测量方法

技术领域

本发明涉及地面测绘技术领域，具体涉及一种倾斜地面坐标测量方法。

背景技术

传统的倾斜测试地面点坐标的测量方法是通过磁传感器x和y轴延伸两个虚拟的点和GPS(Global Positioning System，全球定位系统)点三点定位的原理来计算获取，这种测试地面点的方法需要求解三元二次方程，算法复杂，并且在计算过程中也可能损失地面点坐标的精度。

现在常见的倾斜测试地面点坐标的测量方法是RTK法，RTK(Real-timekinematic，实时动态)载波相位差分技术，是实时处理两个测量站载波相位观测量的差分方法，将基准站采集的载波相位发给用户接收机，进行求差解算坐标。具体为通过RTK获取GPS定位信息和倾角姿态信息，然后利用矢量变换计算地面点坐标。但是此方法受限于大地高H的影响，当H大于一定阈值时，计算误差变大，影响地面点坐标精度。

发明内容

本发明的目的在于解决现有技术直接求解的方法易受到大地高H的影响的问题而设计。

为了解决上述技术问题，本发明所提供的技术方案包括：

一种倾斜地面坐标测量方法，所述方法为：通过RTK法测量待测点，直接计算待测点的坐标，以待测点的坐标的计算结果作为初值进行间接迭代计算，计算出实测大地纬度B，根据实测大地纬度B计算出实测高程H。

进一步的，所述方法中实测大地纬度的计算式为：

B＝atan2(Z+b*e’^2*sin(theta)^3,sqrtxy-a*e^2*cos(theta)^3)；

实测高程H的计算式为：

H＝sqrtxy*cos(B)+Z*sin(B)-N*(1-e^2*sin(B)^2)

其中，atan2为atan2(y，x)函数，表达的是坐标原点为起点，指向(x,y)的射线在坐标平面上与x轴正方向之间的角的角度，X、Y、Z为地心地固坐标系中的坐标值，e为椭球第一偏心率，e’为椭球第二偏心率，a、b分别为椭球的长半袖和短半轴，sqrtxy等于theta等于

进一步的，所述方法中通过RTK法测量待测点，直接计算待测点的坐标的过程具体包括：

步骤1：实时获取RTK天线中的GPS定位点的大地坐标B_gps,L_gps,H_gps；

步骤2：实时获取RTK倾角模块的姿态yaw,pitch,roll；

步骤3：将待测点坐标转换成以GPS点为坐标原点的东北天坐标系下的坐标值，得到E1,N1,U1；

步骤4：将GPS定位点的大地坐标系坐标转换到地心地固坐标系中的坐标值，得到X1,Y1,Z1；

步骤5：计算东北天坐标系在地心地固坐标系下的坐标投影dX,dY,dZ；

步骤6：计算待测点在地心地固坐标系下的坐标值X,Y,Z。

进一步的，所述方法中步骤3的结果具体为：

E1＝(sin(yaw)*sin(roll)+cos(roll)*sin(pitch)*cos(roll))*loc

N1＝(-sin(yaw)*cos(roll)+cos(roll)*sin(pitch)*sin(roll))*loc

U1＝cos(yaw)*cos(pitch)；

loc为测试杆长在U轴的坐标值。

进一步的，所述方法中步骤4的结果具体为：

X1＝(N+h)*cos(B_gps)*cos(L_gps)

Y1＝(N+h)*cos(B_gps)*sin(L_gps)

Z1＝(N*(1-e2)+h)*sin(L_gps)；

h为大地高H_gps，N为卯酉圈曲率半径。

进一步的，所述方法中步骤5的结果具体为：

D＝M*[E,N,U]'

dX＝D1

dY＝D2

dZ＝D3

M为转换矩阵，phi为大地纬度，lambda为大地经度，

M＝[-sin(lambda),-sin(phi)*cos(lambda),cos(phi)*cos(lambda)；

cos(lambda),-sin(phi)*sin(lambda),cos(phi)*sin(lambda)；

0,cos(phi),sin(phi)]。

进一步的，所述方法中步骤6的结果具体为：

X＝X1+dX；Y＝Y1+dY；Z＝Z1+dZ。

本发明的上述技术方案的有益效果如下：

本发明提供一种经过改善的方法，增加间接迭代算法，利用RTK法直接求解的结果作为间接迭代法的初值，进行迭代运算。本发明的方法的精度明显优于现有技术方法的精度，同时不受大地高H的影响。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本发明的实施例，并与说明书一起用于解释本发明的原理。

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

图1为实施例1中RTK法测量并直接求解待测点坐标的流程图。

图2为实施例1中东北天坐标系在地心地固坐标系下的坐标投影示意图。

图3为实施例1中间接迭代法的推导过程示意图。

图4为实施例2中的高度误差示意图。

图5为实施例2中的纬度误差示意图。

具体实施方式

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

实施例1

如图1所示，首先通过RTK法测量并直接求解待测点坐标，步骤包括步骤一至步骤六，具体如下：

步骤一：实时获取RTK天线中的GPS定位点的大地坐标(B_gps,L_gps,H_gps)。

步骤二：实时获取RTK倾角模块的姿态(yaw,pitch,roll)。

步骤三：将待测点坐标转换成以GPS点为坐标原点的东北天坐标系(ENU)下的坐标值，得到(E1,N1,U1)。loc为测试杆长在U轴的坐标值。

E1＝(sin(yaw)*sin(roll)+cos(roll)*sin(pitch)*cos(roll))*loc

N1＝(-sin(yaw)*cos(roll)+cos(roll)*sin(pitch)*sin(roll))*loc

U1＝cos(yaw)*cos(pitch)

步骤四：将GPS定位点的大地坐标系(BLH)坐标转换到地心地固坐标系(ECEF)中的坐标值。得到(X1,Y1,Z1)，h为大地高H_gps，N为卯酉圈曲率半径。

X1＝(N+h)*cos(B_gps)*cos(L_gps)

Y1＝(N+h)*cos(B_gps)*sin(L_gps)

Z1＝(N*(1-e2)+h)*sin(L_gps)

步骤五：如图2所示，计算东北天坐标系(ENU)在地心地固坐标系(ECEF)下的坐标投影(dX,dY,dZ)。

D＝M*[E,N,U]'

dX＝D1

dY＝D2

dZ＝D3

M为转换矩阵，phi为大地纬度，lambda为大地经度。

M＝[-sin(lambda),-sin(phi)*cos(lambda),cos(phi)*cos(lambda)；

cos(lambda),-sin(phi)*sin(lambda),cos(phi)*sin(lambda)；

0,cos(phi),sin(phi)]

步骤六：计算待测点在地心地固坐标系(ECEF)下的坐标值(X,Y,Z)。

X＝X1+dX

Y＝Y1+dY

Z＝Z1+dZ。

然后以上述结果为初值，通过间接迭代法计算实测坐标值，具体步骤及原理如下：

计算待测点在大地坐标系(BLH)下的坐标值。本发明在此处改善了传统算法，并增加了间接迭代方法进行计算。

当：

DB＞DELTA

B1＝atan(Z/sqrtxy*(1+a*e^2*sin(B0)/(Z*sqrt(1-e^2*sin(B0)^2))))

DB＝abs(B1-B0)

B0＝B1

B＝B0

得到：

B＝atan2(Z+b*e’^2*sin(theta)^3，sqrtxy-a*e^2*cos(theta)^3)

L＝atan2(Y，X)

H＝sqrtxy*cos(B)+Z*sin(B)-N*(1-e^2*sin(B)^2)。

如图3所示，上述算法流程的数学计算过程如下：

计算大地经度L，下式中N为卯酉圈曲率半径，e为椭球第一偏心率，e’为椭球第二偏心率。

X＝(N+H)cos B cos L

Y＝(N+H)cos B sinL

Z[N(1-e²)+H]sin B

计算大地纬度B：

a和b分别代表椭球的长短半轴，Ф表示地心纬度，u表示归化纬度，U_p代替U_P’。

r＝X²+Y²

r＝a cos u

Z＝b sin u

r_M＝e acosu_p：

Z_M＝-e’²b sin u_p：

r_M：＝e²a cos³u_Q

Z_M：＝-e’²b sin u_Q

ΔB表示纬度误差，上式中，当H≈2a＞1000km，sin³Bcos³B＝1/8，ΔB≈0″.0018，误差很大。这本质是数学问题所导致。针对这个问题，采用直接解作为初值，结合间接迭代法进行迭代计算，公式如下：

计算大地高H

实施例2

本发明方法在实际测量中使用，通过与测点坐标的已知的实际值比较误差情况。如图4、图5所示，从图中可以看出，本发明提供的测量方法，有效解决了空间直角坐标和大地坐标转换过程中精度受大地高H影响的问题，大地纬度的解算精度提高到10^-7"，大地高的解算精度提高到10^-5m。

以上所述仅是本发明的具体实施方式，使本领域技术人员能够理解或实现本发明。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所申请的原理和新颖特点相一致的最宽的范围。

Claims

1.一种倾斜地面坐标测量方法，其特征在于，所述方法为：

以RTK法测量待测点，直接计算待测点的坐标；

以所述待测点的坐标的计算结果作为初值进行间接迭代计算，计算出实测大地纬度B；

根据实测大地纬度B计算出实测高程H。

2.如权利要求1所述的倾斜地面坐标测量方法，其特征在于，所述方法中实测大地纬度的计算式为：

B＝atan2(Z+b*e’^2*sin(theta)^3,sqrtxy-a*e^2*cos(theta)^3)；

实测高程H的计算式为：

H＝sqrtxy*cos(B)+Z*sin(B)-N*(1-e^2*sin(B)^2)

其中，atan2为atan2(y，x)函数，表达的是坐标原点为起点，指向(x,y)的射线在坐标平面上与x轴正方向之间的角的角度，X、Y、Z为地心地固坐标系中的坐标值，e为椭球第一偏心率，e’为椭球第二偏心率，a、b分别为椭球的长半轴和短半轴，sqrtxy等于theta等于

3.如权利要求2所述的倾斜地面坐标测量方法，其特征在于，所述方法中通过RTK法测量待测点，直接计算待测点的坐标的过程具体包括：

步骤2：实时获取RTK倾角模块的姿态yaw,pitch,roll；

步骤6：计算待测点在地心地固坐标系下的坐标值X,Y,Z。

4.如权利要求3所述的倾斜地面坐标测量方法，其特征在于，所述方法中步骤3的结果具体为：

E1＝(sin(yaw)*sin(roll)+cos(roll)*sin(pitch)*cos(roll))*loc

N1＝(-sin(yaw)*cos(roll)+cos(roll)*sin(pitch)*sin(roll))*loc

U1＝cos(yaw)*cos(pitch)；

loc为测试杆长在U轴的坐标值。

5.如权利要求4所述的倾斜地面坐标测量方法，其特征在于，所述方法中步骤4的结果具体为：

X1＝(N+h)*cos(B_gps)*cos(L_gps)

Y1＝(N+h)*cos(B_gps)*sin(L_gps)

Z1＝(N*(1-e2)+h)*sin(L_gps)；

h为大地高H_gps，N为卯酉圈曲率半径。

6.如权利要求5所述的倾斜地面坐标测量方法，其特征在于，所述方法中步骤5的结果具体为：

D＝M*[E,N,U]'

dX＝D1

dY＝D2

dZ＝D3

M为转换矩阵，phi为大地纬度，lambda为大地经度，

M＝[-sin(lambda),-sin(phi)*cos(lambda),cos(phi)*cos(lambda)；

cos(lambda),-sin(phi)*sin(lambda),cos(phi)*sin(lambda)；

0,cos(phi),sin(phi)]。

7.如权利要求6所述的倾斜地面坐标测量方法，其特征在于，所述方法中步骤6的结果具体为：X＝X1+dX；Y＝Y1+dY；Z＝Z1+dZ。