双向解析光束法

1、双向解析光束法光束法程序有问题，在Getelement这个函数里便出现索引超限，这个问题一直解决不了光束法的流程：1. 根据同名像点对对相交理论求系数阵A，系数阵B和常数阵La11=(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+c3

2、y); 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. 判断改正数的值，如果小于限差则输出结果光束法是最严密的一种方法的原因：在一张相片中，待定点与控制点的像点与摄影中心及相应地面点均构成一条光束，该方法是以每张相片所组成的一束光线作为平差的基本单元，已共线条件方程作为平差的基础方程，通过各个光束在空间中的旋转和平移，使模型之间公共点的光线实现最佳交

3、汇，并使整个区域纳入到已知的地面控制点坐标系中，所以要建立全区域统一的误差方程，整体解求全区域内每张相片的六个外方位元素及所有待定点坐标，光束法区域网平差是基于摄影时像点，物点和摄站点三点共线提出来的。由单张相片构成区域，其平差的数学模型是共线条件方程，平差单元是单个光束，像点坐标是观测值，未知数是每张相片的外方位元素及所有待定点坐标。误差方程直接由像点坐标的观测值列出，能对像点坐标进行系统误差改正。光束法的程序代码为： /计算像片外方位元素，架设phi=0,omega=0,kappa=0求Xs,Ys,Zs /求左片的Xs,Ys,Zs Xsl = (strX0 + strX2) / 2; Ys

4、l = (strY0 + strY2) / 2; L = Math.Pow(Math.Pow(strX0 - strX2, 2) + Math.Pow(strY0 - strY2, 2) + Math.Pow(strZ0 - strZ2, 2), 0.5); l = Math.Pow(Math.Pow(strXl0 - strXl2, 2) + Math.Pow(strYl0 - strYl2, 2), 0.5); H = f * L / l; Zsl = (strZ0 + strZ2) / 2 + H; /计算片左的加密点的物方坐标 Class1 xyz = new Class1(8, 3,

5、 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)

6、 + Zsl); /求右片的Xs,Ys,Zs Xsr = (strX1 + strX3) / 2; Ysr = (strY1 + strY3) / 2; L = Math.Pow(Math.Pow(strX1 - strX3, 2) + Math.Pow(strY1 - strY3, 2) + Math.Pow(strZ1 - strZ3, 2), 0.5); l = Math.Pow(Math.Pow(strXr1 - strXr3, 2) + Math.Pow(strYr1 - strYr3, 2), 0.5); H = f * L / l; Zsr = (strZ1 + strZ3) /

7、 2 + H; /右片加密点 double strxy; strxy = new double24; Console.WriteLine("请输入右片加密点的相片坐标"); for (int i = 0; i < 8; i+) /strxy0 + 3 * i = Convert.ToDouble(Console.ReadLine(); /strxy1 + 3 * i = Convert.ToDouble(Console.ReadLine(); strxy2 + 3 * i = -f; for (int i = 0; i < 24; i+) strxyi = st

8、rxyi * 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 (0,i)/2); uvw.S

9、etElement(1, i, (uvw.GetElement(1, i) + Ysr+UVW.GetElement (1,i)/2); uvw.SetElement(2, i, (uvw.GetElement(2, i) + Zsr+UVW .GetElement (2,i)/2); /从此开始循环 /光束法进行平差 /循环体 do /求左片旋转矩阵R double mtrRl; mtrRl = new double9; mtrRl0 = Math.Cos(phil) * Math.Cos(kappal) - Math.Sin(phil) * Math.Sin(omegal) * Math.

10、Sin(kappal); mtrRl1 = -Math.Cos(phil) * Math.Sin(kappal) - Math.Sin(phil) * Math.Sin(omegal) * Math.Cos(kappal); mtrRl2 = -Math.Sin(phil) * Math.Cos(omegal); mtrRl3 = Math.Cos(omegal) * Math.Sin(kappal); mtrRl4 = Math.Cos(omegal) * Math.Cos(kappal); mtrRl5 = -Math.Sin(omegal); mtrRl6 = Math.Sin(phil

11、) * Math.Cos(kappal) + Math.Cos(phil) * Math.Sin(omegal) * Math.Sin(kappal); mtrRl7 = -Math.Sin(phil) * Math.Sin(kappal) + Math.Cos(phil) * Math.Sin(omegal) * Math.Cos(kappal); mtrRl8 = Math.Cos(phil) * Math.Cos(omegal); /求右片旋转矩阵Rr double mtrRr; mtrRr = new double9; mtrRr0 = Math.Cos(phir) * Math.Co

12、s(kappar) - Math.Sin(phir) * Math.Sin(omegar) * Math.Sin(kappar); mtrRr1 = -Math.Cos(phir) * Math.Sin(kappar) - Math.Sin(phir) * Math.Sin(omegar) * Math.Cos(kappar); mtrRr2 = -Math.Sin(phir) * Math.Cos(omegar); mtrRr3 = Math.Cos(omegar) * Math.Sin(kappar); mtrRr4 = Math.Cos(omegar) * Math.Cos(kappar

13、); mtrRr5 = -Math.Sin(omegar); mtrRr6 = Math.Sin(phir) * Math.Cos(kappar) + Math.Cos(phir) * Math.Sin(omegar) * Math.Sin(kappar); mtrRr7 = -Math.Sin(phir) * Math.Sin(kappar) + Math.Cos(phir) * Math.Sin(omegar) * Math.Cos(kappar); mtrRr8 = Math.Cos(phir) * Math.Cos(omegar); /求系数阵A double mtrA; mtrA =

14、 new double576; for (int i = 0; i < 8; i+) mtrA0 + 48 * i = (mtrRl0 * f + mtrRl2 * xyz.GetElement(i, 0) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrA1 + 48 * i = (mtrRl3 * f + mtrRl5 * xyz.GetElement(i, 0) / (mtrRl2

15、* (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrA2 + 48 * i = (mtrRl6 * f + mtrRl8 * xyz.GetElement(i, 0) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrA3

16、 + 48 * i = xyz.GetElement(i, 1) * Math.Sin(omegal) - (xyz.GetElement(i, 0) * (xyz.GetElement(i, 0) * Math.Cos(kappal) - xyz.GetElement(i, 1) * Math.Sin(kappal) / f) + f * Math.Cos(kappal) * Math.Cos(omegal); mtrA4 + 48 * i = -f * Math.Sin(kappal) - xyz.GetElement(i, 0) * (xyz.GetElement(i, 0) * Mat

17、h.Sin(kappal) + xyz.GetElement(i, 1) * Math.Cos(kappal) / f; mtrA5 + 48 * i = xyz.GetElement(i, 1); for (int j = 0; j < 6; j+) mtrAj + 6 + 48 * i = 0; mtrA12 + 48 * i = (mtrRl1 * f + mtrRl2 * xyz.GetElement(i, 1) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + m

18、trRl8 * (uvw.GetElement(2, i) - Zsl); mtrA13 + 48 * i = (mtrRl4 * f + mtrRl5 * xyz.GetElement(i, 1) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrA14 + 48 * i = (mtrRl7 * f + mtrRl8 * xyz.GetElement(i, 1) / (mtrRl2 * (uvw

19、.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrA15 + 48 * i = -xyz.GetElement(i, 0) * Math.Sin(omegal) - (xyz.GetElement(i, 0) * (xyz.GetElement(i, 0) * Math.Cos(kappal) - xyz.GetElement(i, 1) * Math.Sin(kappal) / f - f * Math.Sin(kappal)

20、 * Math.Cos(omegal); mtrA16 + 48 * i = -f * Math.Cos(kappal) - xyz.GetElement(i, 1) * (xyz.GetElement(i, 0) * Math.Sin(kappal) + xyz.GetElement(i, 1) * Math.Cos(kappal) / f; mtrA17 + 48 * i = -xyz.GetElement(i, 0); for (int j = 0; j < 12; j+) mtrAj + 18 + 48 * i = 0; mtrA30 + 48 * i = (mtrRr0 * f

21、 + mtrRr2 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrA31 + 48 * i = (mtrRr3 * f + mtrRr5 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) +

22、mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrA32 + 48 * i = (mtrRr6 * f + mtrRr8 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrA33 + 48 * i = XYZ.GetElement(i, 1) * Math.Sin(omegar) - (XYZ.GetElement(i

23、, 0) * (XYZ.GetElement(i, 0) * Math.Cos(kappar) - XYZ.GetElement(i, 1) * Math.Sin(kappar) / f) + f * Math.Cos(kappar) * Math.Cos(omegar); mtrA34 + 48 * i = -f * Math.Sin(kappar) - XYZ.GetElement(i, 0) * (XYZ.GetElement(i, 0) * Math.Sin(kappar) + XYZ.GetElement(i, 1) * Math.Cos(kappar) / f; mtrA35 +

24、48 * i = XYZ.GetElement(i, 1); for (int j = 0; j < 6; j+) mtrAj + 36 + 48 * i = 0; mtrA42 + 48 * i = (mtrRr1 * f + mtrRr2 * XYZ.GetElement(i, 1) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrA43 + 48 * i = (mtrRr4 * f

25、+ mtrRr5 * XYZ.GetElement(i, 1) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrA44 + 48 * i = (mtrRr7 * f + mtrRr8 * XYZ.GetElement(i, 1) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + m

26、trRr8 * (uvw.GetElement(2, i) - Zsr); mtrA45 + 48 * i = -XYZ.GetElement(i, 0) * Math.Sin(omegar) - (XYZ.GetElement(i, 0) * (XYZ.GetElement(i, 0) * Math.Cos(kappar) - XYZ.GetElement(i, 1) * Math.Sin(kappar) / f - f * Math.Sin(kappar) * Math.Cos(omegar); mtrA46 + 48 * i = -f * Math.Cos(kappar) - XYZ.G

27、etElement(i, 1) * (XYZ.GetElement(i, 0) * Math.Sin(kappar) + XYZ.GetElement(i, 1) * Math.Cos(kappar) / f; mtrA47 + 48 * i = -XYZ.GetElement(i, 0); for (int i = 8; i < 12; i+) mtrA0 + 48 * i = (mtrRl0 * f + mtrRl2 * strXli - 8) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 *

28、(strZi - 8 - Zsl); mtrA1 + 48 * i = (mtrRl3 * f + mtrRl5 * strXli - 8) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl); mtrA2 + 48 * i = (mtrRl6 * f + mtrRl8 * strXli - 8) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl);

29、 mtrA3 + 48 * i = strYli - 8 * Math.Sin(omegal) - (strXli - 8 * (strXli - 8 * Math.Cos(kappal) - strYli - 8 * Math.Sin(kappal) / f) + f * Math.Cos(kappal) * Math.Cos(omegal); mtrA4 + 48 * i = -f * Math.Sin(kappal) - strXli - 8 * (strXli - 8 * Math.Sin(kappal) + strYli - 8 * Math.Cos(kappal) / f; mtr

30、A5 + 48 * i = strYli - 8; for (int j = 0; j < 6; j+) mtrAj + 6 + 48 * i = 0; mtrA12 + 48 * i = (mtrRl1 * f + mtrRl2 * strYli - 8) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl); mtrA13 + 48 * i = (mtrRl4 * f + mtrRl5 * strYli - 8) / (mtrRl2 * (strXi - 8 -

31、Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl); mtrA14 + 48 * i = (mtrRl7 * f + mtrRl8 * strYli - 8) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl); mtrA15 + 48 * i = -strXli - 8 * Math.Sin(omegal) - (strXli - 8 * (strXli - 8 * Math.Cos(kappal

32、) - strYli - 8 * Math.Sin(kappal) / f - f * Math.Sin(kappal) * Math.Cos(omegal); mtrA16 + 48 * i = -f * Math.Cos(kappal) - strYli - 8 * (strXli - 8 * Math.Sin(kappal) + strYli - 8 * Math.Cos(kappal) / f; mtrA17 + 48 * i = -strXli - 8; for (int j = 0; j < 12; j+) mtrAj + 18 + 48 * i = 0; mtrA30 +

33、48 * i = (mtrRr0 * f + mtrRr2 * strXri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); mtrA31 + 48 * i = (mtrRr3 * f + mtrRr5 * strXri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); mtrA32 + 48 * i = (mtrRr6

34、 * f + mtrRr8 * strXri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); mtrA33 + 48 * i = strYri - 8 * Math.Sin(omegar) - (strXri - 8 * (strXri - 8 * Math.Cos(kappar) - strYri - 8 * Math.Sin(kappar) / f) + f * Math.Cos(kappar) * Math.Cos(omegar); mtrA34 +

35、 48 * i = -f * Math.Sin(kappar) - strXri - 8 * (strXri - 8 * Math.Sin(kappar) + strYri - 8 * Math.Cos(kappar) / f; mtrA35 + 48 * i = strYri - 8; for (int j = 0; j < 6; j+) mtrAj + 36 + 48 * i = 0; mtrA42 + 48 * i = (mtrRr1 * f + mtrRr2 * strYri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi

36、 - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); mtrA43 + 48 * i = (mtrRr4 * f + mtrRr5 * strYri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); mtrA44 + 48 * i = (mtrRr7 * f + mtrRr8 * strYri - 8) / (mtrRr2 * (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mt

37、rRr8 * (strZi - 8 - Zsr); mtrA45 + 48 * i = -strXri - 8 * Math.Sin(omegar) - (strXri - 8 * (strXri - 8 * Math.Cos(kappar) - strYri - 8 * Math.Sin(kappar) / f - f * Math.Sin(kappar) * Math.Cos(omegar); mtrA46 + 48 * i = -f * Math.Cos(kappar) - strYri - 8 * (strXri - 8 * Math.Sin(kappar) + strYri - 8

38、* Math.Cos(kappar) / f; mtrA47 + 48 * i = -strXri - 8; /将数组A转化为矩阵 Class1 A = new Class1(48, 6, mtrA); /求系数阵B double mtrB; mtrB = new double1728; for (int i = 0; i < 1728; i+) mtrBi = 0;/给B矩阵赋值 for (int i = 0; i < 8; i+) mtrB0 + 147 * i = -mtrA0 + 48 * i; mtrB1 + 147 * i = -mtrA1 + 48 * i; mtrB

39、2 + 147 * i = -mtrA2 + 48 * i; mtrB36 + 147 * i = -mtrA12 + 48 * i; mtrB37 + 147 * i = -mtrA13 + 48 * i; mtrB38 + 147 * i = -mtrA14 + 48 * i; mtrB72 + 147 * i = -mtrA30 + 48 * i; mtrB73 + 147 * i = -mtrA31 + 48 * i; mtrB74 + 147 * i = -mtrA32 + 48 * i; mtrB108 + 147 * i = -mtrA42 + 48 * i; mtrB109 +

40、 147 * i = -mtrA43 + 48 * i; mtrB110 + 147 * i = -mtrA44 + 48 * i; /转化矩阵B Class1 B = new Class1(48, 36, mtrB); /求常数项 double mtrL; mtrL = new double48; for (int i = 0; i < 8; i+) mtrL0 + 4 * i = xyz.GetElement(i, 0) + f * (mtrRl0 * (uvw.GetElement(0, i) - Xsl) + mtrRl3 * (uvw.GetElement(1, i) - Ys

41、l) + mtrRl6 * (uvw.GetElement(2, i) - Zsl) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrL1 + 4 * i = xyz.GetElement(i, 1) + f * (mtrRl1 * (uvw.GetElement(0, i) - Xsl) + mtrRl4 * (uvw.GetElement(1, i) - Ysl) + mtrRl7 * (u

42、vw.GetElement(2, i) - Zsl) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) + mtrRl5 * (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl); mtrL2 + 4 * i = XYZ.GetElement(i, 0) + f * (mtrRr0 * (uvw.GetElement(0, i) - Xsr) + mtrRr3 * (uvw.GetElement(1, i) - Ysr) + mtrRr6 * (uvw.GetElement(2,

43、 i) - Zsr) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); mtrL3 + 4 * i = XYZ.GetElement(i, 1) + f * (mtrRr1 * (uvw.GetElement(0, i) - Xsr) + mtrRr4 * (uvw.GetElement(1, i) - Ysr) + mtrRr7 * (uvw.GetElement(2, i) - Zsr) / (mt

44、rRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr); for (int i = 8; i < 12; i+) mtrL0 + 4 * i = strXli - 8 + f * (mtrRl0 * (strXi - 8 - Xsl) + mtrRl3 * (strYi - 8 - Ysl) + mtrRl6 * (strZi - 8 - Zsl) / (mtrRl2 * (strXi - 8 - Xsl) + mtrRl5 * (strYi - 8 - Ysl) + mt