我试图从结果中扣除 2D 转换参数。
给出的是未知 XY 坐标系中的大量样本以及它们在 WGS84(经度、纬度)中的对应物。由于面积很小,我们可以假设目标系统也是平坦的。
可悲的是,我不知道使用了哪种比例、旋转、翻译顺序,我什至不确定是否有 1 个或 2 个翻译。
我试图创建一个冗长的方程系统,但结果太复杂了,我无法处理。基本几何也让我失望了,因为变换的顺序是未知的,我必须检查每一个可能的组合顺序。
是否有解决这个问题的系统方法?
我试图从结果中扣除 2D 转换参数。
给出的是未知 XY 坐标系中的大量样本以及它们在 WGS84(经度、纬度)中的对应物。由于面积很小,我们可以假设目标系统也是平坦的。
可悲的是,我不知道使用了哪种比例、旋转、翻译顺序,我什至不确定是否有 1 个或 2 个翻译。
我试图创建一个冗长的方程系统,但结果太复杂了,我无法处理。基本几何也让我失望了,因为变换的顺序是未知的,我必须检查每一个可能的组合顺序。
是否有解决这个问题的系统方法?
计算比例因子很容易,只需选择任意两个点并在您的 XY 空间和 WGS84 空间中找到它们之间的距离,它们的比率就是您的比例因子。
旋转和平移有点棘手,但当您了解到应用任意数量的旋转或平移(仅限二维!)的结果可以减少到以某个未知角度围绕某个未知点的单个旋转时,难度几乎没有那么大.
突然你有 N 个点来确定 3 个未知数,即旋转轴(x 和 y 坐标)和旋转角度。
计算旋转如下所示:
Pr = R*(Pxy - Paxis_xy) + Paxis_xy
Pr
是您在 XY 空间中的旋转点,然后需要将其转换为 WGS84 空间(如果您的坐标系的轴不同)。
R
是熟悉的旋转矩阵,具体取决于您的旋转角度。
Pxy
是您在 XY 空间中的未旋转点。
Paxis_xy
是 XY 空间中的旋转轴。
要真正找到 3 个未知数,您需要通过您找到的比例因子来取消缩放 WGS84 点(或等效地缩放您的 XY 点)并移动您的点,以便两个坐标系具有相同的原点。
一、求旋转角度:取两对对应的点P1, P1'
,P2, P2'
写出
P1' = R(P1-A) + A
P2' = R(P2-A) + A
为了简洁起见,我在哪里换A = Paxis_xy
了。减去这两个方程得到:
P2'-P1' = R(P2-P1)
B = R * C
Bx = cos(a) * Cx - sin(a) * Cy
By = cos(a) * Cx + sin(a) * Cy
By + Bx = 2 * cos(a) * Cx
(By + Bx) / (2 * Cx) = cos(a)
...
(By - Bx) / (2 * Cy) = sin(a)
a = atan2(sin(a), cos(a)) <-- to get the right quadrant
你有你的角度,你也可以做一个快速检查,cos(a) * cos(a) + sin(a) * sin(a) == 1
以确保你得到了所有的计算正确或者你的系统真的是一个保持方向的等距(只包括平移和旋转)。
既然我们知道a
我们知道R
,所以要找到A
我们这样做:
P1` = R(P1-A) + A
P1' - R*P1 = (I-R)A
A = (inverse(I-R)) * (P1' - R*P1)
其中2x2 矩阵的求逆很容易。
编辑:上面有一个错误,或者更具体地说是一种需要单独处理的情况。
平移和旋转的一种组合不会减少为单个旋转,而是单个平移。您可以将其视为固定点(操作后有多少点不变)。
平移没有固定点(所有点都改变了),旋转有 1 个固定点(轴不改变)。事实证明,两次旋转留下 1 个固定点,一次平移和一次旋转留下 1 个固定点,这(有一点证据表明固定点的数量告诉您执行的操作)是这些任意组合导致的原因单次旋转。
这对你来说意味着如果你的角度出来了,0
那么使用上面的方法也会给你A = 0
,这可能是不正确的。在这种情况下,您必须这样做A = P1' - P1
。
如果我正确理解了这个问题,那么您有n个点 (X 1 ,Y 1 ),...,(X n ,Y n ),对应的点,例如 (x 1 ,y 1 ),...,( x n ,y n ) 在另一个坐标系中,而前者据说是从后者通过旋转、缩放和平移获得的。
请注意,此数据不能确定旋转/缩放的固定点,或“应该”应用操作的顺序。另一方面,如果你事先知道这些或任意选择它们,你会发现一个旋转、平移和缩放因子可以按预期转换数据。
例如,您可以选择一个任意点,例如 p 0 = [X 1 , Y 1 ] T(列向量)作为旋转和缩放的固定点,然后从其他两个点的坐标中减去其坐标得到 p 2 = [X 2 -X 1 , Y 2 -Y 1 ] T和 p 3 = [X 3 -X 1 , Y 3 -Y 1 ] T。还取列向量 q 2 = [x 2 -x 1 , y 2 -y 1 ] T, q 3 = [x 3 -x 1 , y 3 -y 1 ] T。现在 [p 2 p 3 ] = A*[q 2 q 3 ],其中 A 是表示旋转缩放的未知 2x2 矩阵。你可以解决它(除非你不走运并选择退化点)作为 A = [p 2 p 3 ] * [q 2 q 3 ] -1其中-1表示矩阵逆(2x2 矩阵 [q 2 q 3])。现在,如果坐标系之间的变换真的是旋转缩放平移,那么所有点都应该满足 P k = A * (Q k -q 0 ) + p 0,其中 P k = [X k , Y k ] T , Q k = [x k , y k ] T , q 0 =[x 1 , y 1 ] T , k=1,..,n。
如果你愿意,你可以很容易地从 A 的分量中确定缩放和旋转参数,或者结合 b = -A * q 0 + p 0得到 P k = A*Q k + b。
上述方法对噪声或选择退化点反应不佳。如有必要,可以通过应用例如主成分分析来解决此问题,如果 MATLAB 或其他一些线性代数工具可用,这也只是几行代码。