【Walk脚本】最小二乘求解方法——以布尔沙七参数为例

alert

测绘和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}

pengtp16

{:soso_e127:}这么复杂，想学都学习不了。

alert

pengtp16 发表于 2012-10-24 11:17
这么复杂，想学都学习不了。

这个很复杂？