七参数求解

最小二乘求解方法——以布尔沙七参数为例
—Walkinfo—地信网论坛
测绘和GIS计算中经常要用到《最小二乘法》求解。

如坐标转换中的四参数、Bursa 七参数，等等。

若用matlib求解则需要学习其语法等，若使用他人的程序则需要求证自变量和因变量的关系，若自己编写c/c++程序则颇费周折。

求人不如求己，在Walk脚本中提供了类似于 matlib的矩阵运算功能。

以求解Bursa七参数为例，其方程如下：
r1 + 0 + 0 + 0 + r5*z0 - r6*y0 + r7*x0 = x1 - x0
0 + r2 + 0 - r4*z0 + 0 + r6*x0 + r7*y0 = y1 - y0
0 + 0 + r3 + r4*y0 - r5*x0 + 0 + r7*z0 = z1 - z0
式中，源srcC(x0,y0,z0), 目标tarC(x1,y1,z1)，r1,r2,...,r7为七参数。

bool solve7X(array &srcC, array &tarC, array &X)
{
int n=srcC.getSize();
//组成系数方程矩阵 A*X=L:
double A[3*n][7];
double L[3*n][1];
for (int ii=0, k=0; ii<n; ii++)
{
double x=srcC[ii][0], y=srcC[ii][1], z=srcC[ii][2]; A[k][0]=1.; A[k][1]=0.; A[k][2]=0.; A[k][3]=0.; A[k][4]=z; A[k][5]=-y; A[k][6]=x;
L[k++][0]=tarC[ii][0] - x;
A[k][0]=0.; A[k][1]=1.; A[k][2]=0.; A[k][3]=-z; A[k][4]=0.; A[k][5]=x; A[k][6]=y;
L[k++][0]=tarC[ii][1] - y;
A[k][0]=0.; A[k][1]=0.; A[k][2]=1.; A[k][3]=y; A[k][4]=-x; A[k][5]=0.; A[k][6]=z;
L[k++][0]=tarC[ii][2] - z;
}
//按最小二乘法组成法方程： (A'*A)*X = (A'*L) -- A'为A的转置矩阵
// 设： N=(A'*A), W=(A'*L) 法方程为 N*X=W
array At=A.mtTranspose(); // A'
array N=At.mtMultiply(A); // (A'*A)
array W, W0=At.mtMultiply(L); // (A'*L)
for (ii=0; ii<W0.getSize(); ii++) W.add(W0[ii][0]);
//高斯主元素消元法求X
int o1 = X.mtGauss(N,W);
return (o1 == 1);
}
//加入自己的重合点坐标，运行
void main()
{
//目标点坐标
array tarC={// N, E, H
{2457352.726, 524988.625, 0.},
{2459025.942, 534157.249, 0.},
{2460338.594, 529155.636, 0.}
};
//源点坐标
array srcC={// N, E, H
{2457426.677, 525103.512, 0.},
{2459100.114, 534273.255, 0.},
{2460412.919, 529271.031, 0.}
};
array X;
if (solve7X(srcC, tarC, X))
trace("\n最小二乘解 X = " + toString(X)+"\n");
}
运行结果：
最小二乘解 X = [7]{227.642618576805, -54.7105352482602, 0.0, 0.0, 0.0, 1.66072305203775e-6, -0.000122372906272742}。