双向解析光束法

双向解析光束法

光束法程序有问题，在Getelement这个函数里便出现索引超限，这个问题一直解决不了

光束法的流程：

1.根据同名像点对对相交理论求系数阵A，系数阵B和常数阵L

a11=(a1f+a3x)/Z; a12=(b1f+b3x)/Z; a13=(c1f+c3x)/Z; a14=ysin(omega)-[x/f(xcos(kappa)-ysin(kappa))+fcos(kappa)]cos(omega);

a15=-fsin(kappa)-x/f(xsin(kappa)+ycos(kappa));

a16=y;

a21=(a2f+a3y)/Z; a22=(b2f+b3y)/Z; a23=(c2f+c3y); a24=-xsin(omega)-[x/f(xcos(kappa)-ysin(kappa))-fsin(kappa)]cos(omega)

a25=-fcos(kappa)-y/f(xsin(kappa)+ycos(kappa));

a26=-x;

2.求方程的改化法方程求出外方位元素和物方坐标改正数

3.判断改正数的值，如果小于限差则输出结果

光束法是最严密的一种方法的原因：

在一张相片中，待定点与控制点的像点与摄影中心及相应地面点均构成一条光束，该方法是以每张相片所组成的一束光线作为平差的基本单元，已共线条件方程作为平差的基础方程，通过各个光束在空间中的旋转和平移，使模型之间公共点的光线实现最佳交汇，并使整个区域纳入到已知的地面控制点坐标系中，所以要建立全区域统一的误差方程，整体解求全区域内每张相片的六个外方位元素及所有待定点坐标，光束法区域网平差是基于摄影时像点，物点和摄站点三点共线提出来的。由单张相片构成区域，其平差的数学模型是共线条件方程，平差单元是单个光束，像点坐标是观测值，未知数是每张相片的外方位元素及所有待定点坐标。误差方程直接由像点坐标的观测值列出，能对像点坐标进行系统误差改正。

光束法的程序代码为：

//计算像片外方位元素，架设phi=0,omega=0,kappa=0求Xs,Ys,Zs

//求左片的Xs,Ys,Zs

Xsl = (strX[0] + strX[2]) / 2;

Ysl = (strY[0] + strY[2]) / 2;

L = Math.Pow(Math.Pow(strX[0] - strX[2], 2) + Math.Pow(strY[0] - strY[2], 2) + Math.Pow(strZ[0] - strZ[2], 2), 0.5);

l = Math.Pow(Math.Pow(strXl[0] - strXl[2], 2) + Math.Pow(strYl[0] - strYl[2], 2),

0.5);

H = f * L / l;

Zsl = (strZ[0] + strZ[2]) / 2 + H;

//计算片左的加密点的物方坐标

Class1 xyz = new Class1(8, 3, strXY);

Class1 xyzt=xyz.Transpose ();

//phi,omega,kappa为零，故旋转矩阵为单位阵

//像空间辅助坐标系中的坐标

Class1 UVW = xyzt.Multiply(H);

//求地面摄影测量坐标系中的坐标

for (int i = 0; i < 8; i++)

{

UVW.SetElement(0, i, UVW.GetElement(0, i) + Xsl);

UVW.SetElement(1, i, UVW.GetElement(1, i) + Ysl);

UVW.SetElement(2, i, UVW.GetElement(2, i) + Zsl);

}

//求右片的Xs,Ys,Zs

Xsr = (strX[1] + strX[3]) / 2;

Ysr = (strY[1] + strY[3]) / 2;

L = Math.Pow(Math.Pow(strX[1] - strX[3], 2) + Math.Pow(strY[1] - strY[3], 2) + Math.Pow(strZ[1] - strZ[3], 2), 0.5);

l = Math.Pow(Math.Pow(strXr[1] - strXr[3], 2) + Math.Pow(strYr[1] - strYr[3], 2),

0.5);

H = f * L / l;

Zsr = (strZ[1] + strZ[3]) / 2 + H;

//右片加密点

double[] strxy;

strxy = new double[24];

Console.WriteLine("请输入右片加密点的相片坐标");

for (int i = 0; i < 8; i++)

{

//strxy[0 + 3 * i] = Convert.ToDouble(Console.ReadLine());

//strxy[1 + 3 * i] = Convert.ToDouble(Console.ReadLine());

strxy[2 + 3 * i] = -f;

}

for (int i = 0; i < 24; i++)

{

strxy[i] = strxy[i] * 0.000025;

}

//计算片左的加密点的物方坐标

Class1 XYZ = new Class1(8, 3, strxy);

Class1 XYZt=XYZ .Transpose ();

//phi,omega,kappa为零，故旋转矩阵为单位阵

//像空间辅助坐标系中的坐标

Class1 uvw = XYZt.Multiply(H);

//求地面摄影测量坐标系中的坐标

for (int i = 0; i < 8; i++)

{

uvw.SetElement(0, i, (uvw.GetElement(0, i) + Xsr+UVW.GetElement