坐标系转换--平面多项式拟合公式转换模型变换关系理解

平面多项式拟合公式转换模型变换关系理解

{ X 2 = X 1 + Δ X Y 2 = Y 1 + Δ Y (1) \tag{1} \begin{dcases} X_2 = X_1 + \varDelta{X} \\ Y_2 = Y_1 + \varDelta{Y} \\ \end{dcases} {X2=X1+ΔXY2=Y1+ΔY(1)

式中: 式中: 式中:
X 1 , Y 1 :原平面直角坐标系; X_1,Y_1:原平面直角坐标系; X1,Y1:原平面直角坐标系;
X 2 , Y 2 :新平面直角坐标系; X_2,Y_2:新平面直角坐标系; X2,Y2:新平面直角坐标系;
Δ X , Δ Y :坐标转换改正量,计算公式为: \varDelta{X},\varDelta{Y}:坐标转换改正量,计算公式为: ΔX,ΔY:坐标转换改正量,计算公式为:
Δ X = a 0 + a 1 X + a 2 Y + a 3 X 2 + a 4 X Y + a 5 Y 2 + a 6 X 3 + a 7 X 2 Y + a 8 X Y 2 + a 9 Y 3 + ⋯ Δ Y = b 0 + b 1 X + b 2 Y + b 3 X 2 + b 4 X Y + b 5 Y 2 + b 6 X 3 + b 7 X 2 Y + b 8 X Y 2 + b 9 Y 3 + ⋯ \varDelta{X}=a_0+a_1X+a_2Y+a_3X^2+a_4XY+a_5Y^2+a_6X^3+a_7X^2Y+a_8XY^2+a_9Y^3+\cdots \\ \varDelta{Y}=b_0+b_1X+b_2Y+b_3X^2+b_4XY+b_5Y^2+b_6X^3+b_7X^2Y+b_8XY^2+b_9Y^3+\cdots \\ ΔX=a0+a1X+a2Y+a3X2+a4XY+a5Y2+a6X3+a7X2Y+a8XY2+a9Y3+ΔY=b0+b1X+b2Y+b3X2+b4XY+b5Y2+b6X3+b7X2Y+b8XY2+b9Y3+
其中, a i , b i 为系数 ( i = 0 , 1 , 2 , ⋯   ) ,通过最小二乘求解。 其中,a_i,b_i为系数(i=0,1,2,\cdots),通过最小二乘求解。 其中,ai,bi为系数(i=0,1,2,),通过最小二乘求解。

( 1 ) 式转为矩阵方程式为: (1)式转为矩阵方程式为: (1)式转为矩阵方程式为:
[ X 2 Y 2 ] = [ X 1 Y 1 ] + [ 1 X Y X 2 X Y Y 2 X 3 X 2 Y X Y 2 Y 3 ⋯ 0 0 0 0 0 0 0 0 0 0 ⋯ 0 0 0 0 0 0 0 0 0 0 ⋯ 1 X Y X 2 X Y Y 2 X 3 X 2 Y X Y 2 Y 3 ⋯ ] [ a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 ⋮ b 0 b 1 b 2 b 3 b 4 b 5 b 6 b 7 b 8 b 9 ⋮ ] \begin{bmatrix} X_2 \\ Y_2 \end{bmatrix} = \begin{bmatrix} X_1 \\ Y_1 \end{bmatrix} + \begin{bmatrix} 1 & X & Y & X^2 & XY & Y^2 & X^3 & X^2Y & XY^2 & Y^3 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 1 & X & Y & X^2 & XY & Y^2 & X^3 & X^2Y & XY^2 & Y^3 & \cdots \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \\ a_8 \\ a_9 \\ \vdots \\ b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \\ b_5 \\ b_6 \\ b_7 \\ b_8 \\ b_9 \\ \vdots \end{bmatrix} [X2Y2]=[X1Y1]+[10X0Y0X20XY0Y20X30X2Y0XY20Y30010X0Y0X20XY0Y20X30X2Y0XY20Y3] a0a1a2a3a4a5a6a7a8a9b0b1b2b3b4b5b6b7b8b9

基于最小二乘与多对同名点对计算参数

设存在 n 对同名点对: ( X a , Y a ) 1 → ( X b , Y b ) 1 , ⋯ , ( X a , Y a ) n → ( X b , Y b ) n . 设存在n对同名点对:(X_a,Y_a)_1 \rarr (X_b,Y_b)_1,\cdots,(X_a,Y_a)_n \rarr (X_b,Y_b)_n. 设存在n对同名点对:(Xa,Ya)1(Xb,Yb)1(Xa,Ya)n(Xb,Yb)n.
令 令
θ = [ a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 ⋮ b 0 b 1 b 2 b 3 b 4 b 5 b 6 b 7 b 8 b 9 ⋮ ] \theta = \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \\ a_8 \\ a_9 \\ \vdots \\ b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \\ b_5 \\ b_6 \\ b_7 \\ b_8 \\ b_9 \\ \vdots \end{bmatrix} θ= a0a1a2a3a4a5a6a7a8a9b0b1b2b3b4b5b6b7b8b9
v i = ( X b − X a , Y b − Y a ) i T , v_i=(X_b - X_a,Y_b - Y_a)^T_i, vi=(XbXa,YbYa)iT
P i = [ 1 X Y X 2 X Y Y 2 X 3 X 2 Y X Y 2 Y 3 ⋯ 0 0 0 0 0 0 0 0 0 0 ⋯ 0 0 0 0 0 0 0 0 0 0 ⋯ 1 X Y X 2 X Y Y 2 X 3 X 2 Y X Y 2 Y 3 ⋯ ] i , P_i= \begin{bmatrix} 1 & X & Y & X^2 & XY & Y^2 & X^3 & X^2Y & XY^2 & Y^3 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 1 & X & Y & X^2 & XY & Y^2 & X^3 & X^2Y & XY^2 & Y^3 & \cdots \\ \end{bmatrix}_i, Pi=[10X0Y0X20XY0Y20X30X2Y0XY20Y30010X0Y0X20XY0Y20X30X2Y0XY20Y3]i
i = 1 , ⋯   , n i=1,\cdots,n i=1,,n

根据式 ( 3 ) ,代入样本值得到方程组如下: 根据式(3),代入样本值得到方程组如下: 根据式(3),代入样本值得到方程组如下:
{ P 1 θ = v 1 P 2 θ = v 2 ⋮ P n θ = v n \begin{dcases} P_1\theta = v_1 \\ P_2\theta = v_2 \\ \vdots \\ P_n\theta = v_n \end{dcases} P1θ=v1P2θ=v2Pnθ=vn
则变换为矩阵方程为: 则变换为矩阵方程为: 则变换为矩阵方程为:
v = P θ v = P\theta v=
P = [ P 1 P 2 ⋮ P n ] , v = [ v 1 v 2 ⋮ v n ] P= \begin{bmatrix} P_1 \\ P_2 \\ \vdots \\ P_n \end{bmatrix}, v= \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ v_n \end{bmatrix} P= P1P2Pn v= v1v2vn

考虑 v = P θ 无解,需要从 P 的列空间中找出最接近 v 的向量 u ( u 可以理解为 v 在 P 的列空间中的投影,理解如下图所示:) 考虑v = P\theta无解,需要从P的列空间中找出最接近v的向量u(u可以理解为v在P的列空间中的投影,理解如下图所示:) 考虑v=无解,需要从P的列空间中找出最接近v的向量uu可以理解为vP的列空间中的投影,理解如下图所示:)
在这里插入图片描述

如上图所示, p 是 b 在 [ a 1 a 2 ] 列空间中的投影。 如上图所示,p是b在\begin{bmatrix} a_1 & a_2 \end{bmatrix} 列空间中的投影。 如上图所示,pb[a1a2]列空间中的投影。
令 e = v − u ,最小二乘就是找到 ∥ e ∥ 2 最小的点,最小二乘就是指向量长度的最小平方。 令e=v-u,最小二乘就是找到\parallel e \parallel^2最小的点,最小二乘就是指向量长度的最小平方。 e=vu,最小二乘就是找到e2最小的点,最小二乘就是指向量长度的最小平方。

由上可知, u 位于 P 的列空间中,即 u 是 P 的各列的线性组合: 由上可知,u位于P的列空间中,即u是P的各列的线性组合: 由上可知,u位于P的列空间中,即uP的各列的线性组合:
令 P 的列空间为 P = [ C 1 C 2 ⋯ C m ] 令P的列空间为 P= \begin{bmatrix} C_1 & C_2 & \cdots & C_m \end{bmatrix} P的列空间为P=[C1C2Cm]
故存在 u = C 1 θ 1 ~ + C 2 θ 2 ~ + ⋯ + C m θ m ~ 故存在 u=C_1\tilde{\theta_1} + C_2\tilde{\theta_2} + \cdots + C_m\tilde{\theta_m} 故存在u=C1θ1~+C2θ2~++Cmθm~
即 P θ ~ = u 有解。 即P\tilde{\theta}=u有解。 Pθ~=u有解。

e = v − u = v − P θ ~ e=v-u=v-P\tilde{\theta} e=vu=vPθ~
e 正交于 P 的列空间,存在: e正交于P的列空间,存在: e正交于P的列空间,存在:
e ⊥ C 1 , e ⊥ C 2 , ⋯   , e ⊥ C m e \perp C_1,e \perp C_2,\cdots,e \perp C_m eC1,eC2,,eCm

由向量点积关系式可得: 由向量点积关系式可得: 由向量点积关系式可得:

⇒ { C 1 T ( v − P θ ~ ) = 0 C 2 T ( v − P θ ~ ) = 0 ⋮ C m T ( v − P θ ~ ) = 0 \Rarr \begin{dcases} C_1^T(v-P\tilde{\theta})=0 \\ C_2^T(v-P\tilde{\theta})=0 \\ \vdots \\ C_m^T(v-P\tilde{\theta})=0 \end{dcases} C1T(vPθ~)=0C2T(vPθ~)=0CmT(vPθ~)=0

⇒ [ C 1 T C 2 T C 3 T ⋮ C m T ] ( v − P θ ~ ) = [ 0 0 0 ⋮ 0 ] \Rarr \begin{bmatrix} C_1^T \\ C_2^T \\ C_3^T \\ \vdots \\ C_m^T \end{bmatrix} (v-P\tilde{\theta})= \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} C1TC2TC3TCmT (vPθ~)= 0000

∵ P = [ C 1 C 2 ⋯ C m ] \because P= \begin{bmatrix} C_1 & C_2 & \cdots & C_m \end{bmatrix} P=[C1C2Cm]
∴ P T = [ C 1 T C 2 T ⋮ C m T ] \therefore P^T = \begin{bmatrix} C_1^T \\ C_2^T \\ \vdots \\ C_m^T \end{bmatrix} PT= C1TC2TCmT

⇒ P T ( v − P θ ~ ) = 0 \Rarr P^T(v-P\tilde{\theta})=0 PT(vPθ~)=0
⇒ P T P θ ~ = P T v \Rarr P^TP\tilde{\theta}=P^Tv PTPθ~=PTv
⇒ θ ~ = ( P T P ) − 1 P T v \Rarr \tilde{\theta}=(P^TP)^{-1}P^Tv θ~=(PTP)1PTv

即 θ ~ = ( P T P ) − 1 P T v 为基于最小二乘计算出来的最接近实际参数的转换值 即\tilde{\theta}=(P^TP)^{-1}P^Tv为基于最小二乘计算出来的最接近实际参数的转换值 θ~=(PTP)1PTv为基于最小二乘计算出来的最接近实际参数的转换值

  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
本软件是“测量计算工具包软件”的全面升级版。升级后的软件强化了坐标转换的功能,精简了其他不大使用的功能,软件名称更改为“坐标转换”,2013是全面升级后的第一个版本。 为适应国家测绘局地理信息办公室《2000国家大地坐标系推广使用技术指南》(以下简称《指南》)和《大地测量控制点坐标转换技术规程》(以下简称《规程》)的要求,坐标转换2013除保留原有的布尔沙模型和二维四参数模型外,增加了三维七参数、二维七参数、三维四参数和多项式拟合模型。另外,在转换参数的表达形式上也进行了调正,将“尺度比”改为“尺度变化”,与《指南》和《规程》保持一致。 升级后的坐标转换软件对程序界面和代码也进行了优化,参数的数值表示方式由固定宽度改为科学表示方式,使得其计算精度更高。 升级前的“椭球间的坐标转换”对应于升级后的“布尔沙模型”,升级前的“多公共点相似变换”对应于升级后的“二维四参数模型”。这两种模型升级前的转换参数完全可以用于升级后的软件,仅需将将“尺度比”换算为“尺度变化”即可,换算公式为:尺度变化D=尺度比K-1。 如果用户拥有转换区域的公共点(《指南》和《规程》叫“重合点”)的话,建议用升级后的软件重新计算转换参数。 必须说明的是,不同的转换模型转换参数是不能互换的。 本软件的所有转换模型的计算公式都来源于《指南》和《规程》,仅对“多项式拟合公式的表达形式进行了格式上的统一。 坐标转换2014版增加了GPS高程拟合和墨卡托投影正反算转换
在二维坐标系统中进行多项式拟合,可以使用numpy库的polyfit函数。假设已知一组原始坐标点,需要对其进行多项式拟合,并将其转换到另一个坐标系中。以下是具体实现步骤: 1. 读取原始坐标数据,并将其存储为numpy数组。 2. 根据需要进行多项式拟合,可以使用numpy的polyfit函数进行拟合。该函数的参数包括输入数据、多项式阶数和权重。 3. 计算拟合后的坐标点在新坐标系中的位置。假设需要将原始坐标系转换为新坐标系,可以通过以下公式进行转换: ``` x_new = a*x + b*y + c y_new = d*x + e*y + f ``` 其中,a、b、c、d、e、f是转换矩阵的元素,可以通过求解线性方程组得到。 4. 绘制拟合后的曲线或散点图,并将其转换到新坐标系中。 下面是代码示例: ```python import numpy as np import matplotlib.pyplot as plt # 读取原始坐标数据 data = np.loadtxt('coords.txt') x = data[:, 0] y = data[:, 1] # 进行多项式拟合 z = np.polyfit(x, y, 3) p = np.poly1d(z) # 定义转换矩阵 a = 1.0 b = 0.5 c = 10.0 d = -0.5 e = 1.0 f = 5.0 # 计算拟合后的坐标点在新坐标系中的位置 x_new = a*x + b*y + c y_new = d*x + e*y + f # 绘制拟合曲线或散点图,并将其转换到新坐标系中 xp = np.linspace(x.min(), x.max(), 100) plt.plot(x_new, p(x), '-', xp*a+b*p(xp)+c, d*xp+e*p(xp)+f, '-') plt.show() ``` 其中,coords.txt是存储原始坐标数据的文本文件,每行包括一个点的x坐标和y坐标。np.polyfit函数的第三个参数3表示进行三次多项式拟合,可以根据实际情况进行调整。转换矩阵的元素可以根据实际情况进行修改。最后的绘图代码可以根据需要进行修改,例如添加标题、坐标轴标签等。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值