[发明专利]一种工测数据与三维数字地球高程拟合的方法

权利要求书

1.一种工测数据与三维数字地球高程拟合的方法，其特征在于，包括：

步骤101，工测数据相对坐标与工测控制点投影坐标相互修正；

步骤102，根据工测投影坐标和三维数字地球坐标计算七参数；

步骤103，使用七参数将工测投影坐标中每个工测点的坐标转换到目标数字地球的坐标系统下获得新的工测投影坐标，以新的工测投影坐标进行地形拟合，将工测数据高程调整到三维数字地球高程附近；

步骤104，以新的工测投影坐标修正三维数字地球高程坐标；

所述步骤101，进一步包括：

获得工测点相对坐标；

判断工测点相对坐标与工测控制点投影坐标是否吻合；如果不吻合，则根据工测点相对坐标推算新的投影坐标；

利用平差法修正新的投影坐标；

所述获得工测相对坐标，包括：

已知有n个桩分别为Z₁、Z₂、……、Z_n；档距分别为l₁、l₂、……、l_n-1；转角分别为α₁、α₂、……、α_n；高程分别为h₁、h₂、……、h_n；

测得桩的投影坐标分别为(x₁,y₁,z₁)、(x₂,y₂,z₂)、……、(x_n,y_n,z_n)；

所述根据工测点相对坐标推算新的投影坐标，利用平差法修正新的投影坐标，包括：

步骤1：生成中间坐标：

中间坐标记为(x′₁,y′₁,z′₁)、(x′₂,y′₂,z′₂)、……、(x′_n,y′_n,z′_n)，则生成方法如下：

记ΔH＝z₁-h₁，中间坐标递推公式如下：

$\left\{\begin{array}{c}(x_1^{\′},y_1^{\′})=(x_1,y_1)\\ z_1^{\′}=h_1+\mathrm{\Δ}H=z_1\\ (u_1,v_1)=\frac{(x_2-x_1,y_2-y_1)}{l}\end{array}\right.$

$\begin{array}{c}\left\{\begin{array}{c}(x_2^{\′},y_2^{\′})=(x_1^{\′},y_1^{\′})+(u_1,v_1)\·l_1\\ z_2^{\′}=h_2+\mathrm{\Δ}H\\ (u_2,v_2)=(u_1,v_1)\·\cos \left({\mathrm{\α}}_2\right)+(-v_1,u_1)\·sin\left({\mathrm{\α}}_2\right)\end{array}\right.\\ \mathrm{......}\\ \left\{\begin{array}{c}(x_n^{\′},y_n^{\′})=(x_{n-1}^{\′},y_{n-1}^{\′})+(u_{n-1},v_{n-1})\·l_{n-1}\\ z_n^{\′}=h_n+\mathrm{\Δ}H\\ (u_n,v_n)=(u_{n-1},v_{n-1})\·\cos \left({\mathrm{\α}}_n\right)+(-v_{n-1},u_{n-1})\·\sin \left({\mathrm{\α}}_n\right)\end{array}\right.\end{array}$

其中，ΔH为第一个桩投影坐标高程与工测坐标高程之差，u_i,v_i为中间坐标下第i档的单位方向向量即第i个桩到第i+1个桩的单位方向向量，u₁,v₁由投影坐标确定，第i档的单位方向向量由第i-1档的单位方向向量及第i档起始桩转角确定；

步骤2、进行平差处理生成新投影坐标：

新投影坐标记为(x″₁,y″₁,z″₁)、(x″₂,y″₂,z″₂)、……、(x″_n,y″_n,z″_n)，则新投影坐标与中间坐标的转换公式如下：

$\left\{\begin{array}{c}\left(\begin{array}{c}{x}_i^{\′\′}\\ y_i^{\′\′}\end{array}\right)=\left(\begin{array}{cc}\cos \left(\mathrm{\θ}\right)& -\sin \left(\mathrm{\θ}\right)\\ \sin \left(\mathrm{\θ}\right)& \cos \left(\mathrm{\θ}\right)\end{array}\right)\left(\begin{array}{c}{x}_i^{\′}\\ y_i^{\′}\end{array}\right)+\left(\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\end{array}\right)\\ z_i^{\′\′}=z_i^{\′}+\mathrm{\Δ}h\end{array}\right.$

其中i∈{1,2,...,n}，θ、Δx、Δy、Δh为待求参数；

为了求解θ、Δx、Δy、Δh，需使的以下两个函数值最小：

$\left\{\begin{array}{c}f(\mathrm{\θ},\mathrm{\Δ}x,\mathrm{\Δ}y)=\underset{i=1}{\overset{n}{\Σ}}\[{(x_i^{\′\′}-x_i)}^2+{(y_i^{\′\′}-y_i)}^2\]\\ g\left(\mathrm{\Δ}h\right)=\underset{i=1}{\overset{n}{\Σ}}{(z_i^{\′\′}-z_i)}^2\end{array}\right.$

使用最小二乘法求解θ、Δx、Δy、Δh，其中θ是坐标系旋转角度，Δx、Δy、Δh为坐标的平移值；

所述步骤102根据工测投影坐标和三维数字地球坐标计算七参数，进一步包括：

工测投影坐标中工测点的WGS84经纬度高程为(B₁,L₁,H₁)、(B₂,L₂,H₂)、……、(B_n,L_n,H_n)；

把测得的WGS84经纬度高程转换为WGS84椭球下的空间直角坐标(a₁,b₁,c₁)、(a₂,b₂,c₂)、……、(a_n,b_n,c_n)；

把新投影坐标转换为其所对应椭球下的空间直角坐标(a′₁,b′₁,c′₁)、(a′₂,b′₂,c′₂)、……、(a′_n,b′_n,c′_n)；

七参数转换公式为：

其中i∈{1,2,...,n}，Dx、Dy、Dz、Rx、Ry、Rz、m为待求参数；

为了求解Dx、Dy、Dz、Rx、Ry、Rz、m，需使的以下函数值最小：

$h(Dx,Dy,Dz,Rx,Ry,Rz,m)=\underset{i=1}{\overset{n}{\Σ}}\[{(a_i^{\′\′}-a_i)}^2+{(b_i^{\′\′}-b_i)}^2+{(c_i^{\′\′}-c_i)}^2\]$

使用最小二乘法求解七参数，包括：Dx、Dy、Dz、Rx、Ry、Rz、m。