双向解析光束法上课讲义.doc

双向解析光束法双向解析光束法光束法程序有问题，在 Geteleme nt这个函数里便出现索引超限，这个问题一直 解决不了光束法的流程：1•根据同名像点对对相交理论求系数阵 A，系数阵B和常数阵La11=(a1f+a3x)/Z; a12=(b1f+b3x)/Z; a13=(c1f+c3x)/Z; a14=ysin(omega)-[x/f(xcos(kappa)-ysi n(kappa))+fcos(kappa)]cos(omega);a15=-fsi n(kappa)-x/f(xsi n(kappa)+ycos(kappa));a16=y;a21=(a2f+a3y)/Z; a22=(b2f+b3y)/Z; a23=(c2f+c3y); a24=-xsin(omega)-[x/f(xcos(kappa)-ysi n(kappa))-fsi n(kappa)]cos(omega)a25=-fcos(kappa)-y/f(xsi n(kappa)+ycos(kappa));a26=-x;2. 求方程的改化法方程求出外方位元素和物方坐标改正数3. 判断改正数的值，如果小于限差则输出结果光束法是最严密的一种方法的原因：在一张相片中，待定点与控制点的像点与摄影中心及相应地面点均构成一条光 束，该方法是以每张相片所组成的一束光线作为平差的基本单元，已共线条件 方程作为平差的基础方程，通过各个光束在空间中的旋转和平移，使模型之间 公共点的光线实现最佳交汇，并使整个区域纳入到已知的地面控制点坐标系 中，所以要建立全区域统一的误差方程，整体解求全区域内每张相片的六个外 方位元素及所有待定点坐标，光束法区域网平差是基于摄影时像点，物点和摄 站点三点共线提出来的。

由单张相片构成区域，其平差的数学模型是共线条件 方程，平差单元是单个光束，像点坐标是观测值，未知数是每张相片的外方位 元素及所有待定点坐标误差方程直接由像点坐标的观测值列出，能对像点坐 标进行系统误差改正光束法的程序代码为：〃计算像片外方位元素，架设 phi=0,omega=0,kappa=0求Xs,Ys,Zs〃求左片的Xs,Ys,ZsXsl = (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(strXI[0] - strXI[2], 2) + Math.Pow(strYI[0] - strYI[2], 2), 0.5);H = f * L / l;Zsl = (strZ[0] + strZ[2]) / 2 + H;〃计算片左的加密点的物方坐标Class1 xyz = new Class1(8, 3, strXY);Classi xyzt=xyz.Tra nspose ();〃phi,omega,kappa为零，故旋转矩阵为单位阵//像空间辅助坐标系中的坐标Classi 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,ZsXsr = (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 /I;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] = Con vert.ToDouble(Co nsole.ReadLi ne());//strxy[1 + 3 * i] = Con vert.ToDouble(Co nsole.ReadLi ne()); strxy[2 + 3 * i] = -f;}for (int i = 0; i < 24; i++){strxy[i] = strxy[i] * 0.000025;}〃计算片左的加密点的物方坐标Classi XYZ = new Class1(8, 3, strxy);Classi XYZt=XYZ .Tra nspose ();〃phi,omega,kappa为零，故旋转矩阵为单位阵//像空间辅助坐标系中的坐标Classi uvw = XYZt.Multiply(H);〃求地面摄影测量坐标系中的坐标for (int i = 0; i < 8; i++){uvw.SetEleme nt(0, i, (uvw.GetEleme nt(0, i) + Xs叶 UVW.GetEleme nt (0,i))/2);uvw.SetEleme nt(1, i, (uvw.GetEleme nt(1, i) + Ys叶 UVW.GetEleme nt(1.1) )/2);uvw.SetElement(2, i, (uvw.GetElement(2, i) + Zsr+UVW .GetElement(2.1) )/2);}//从此开始循环〃光束法进行平差//循环体do{〃求左片旋转矩阵Rdouble]] mtrRI;mtrRI = new double[9];mtrRI[0] = Math.Cos(phil) * Math.Cos(kappal) - Math.Sin(phil) * Math.S in (omegal) * Math.Si n( kappal);mtrRl[1] = -Math.Cos(phil) * Math.Si n( kappal) - Math.Si n( phil) * Math.Si n(omegal) * Math.Cos(kappal);mtrRl[2] = -Math.Sin(phil) * Math.Cos(omegal);mtrRl[3] = Math.Cos(omegal) * Math.Si n(kappal);mtrRl[4] = Math.Cos(omegal) * Math.Cos(kappal);mtrRl [5] = -Math.Si n(omegal);mtrRl[6] = Math.Sin(phil) * Math.Cos(kappal) + Math.Cos(phil) ” Math.S in (omegal) * Math.Si n( kappal);mtrRl[7] = -Math.Sin(phil) * Math.Sin(kappal) + Math.Cos(phil) ” Math.Si n(omegal) * Math.Cos(kappal);mtrRl[8] = Math.Cos(phil) * Math.Cos(omegal);〃求右片旋转矩阵Rrdouble[] mtrRr;mtrRr = new double[9];mtrRr[0] = Math.Cos(phir) * Math.Cos(k appar) - Math.Si n(phir) Math.Sin(omegar) * Math.Sin(kappar);mtrRr[1] = -Math.Cos(phir) * Math.Si n(kappar) - Math.Si n(phir)Math.Sin(omegar) * Math.Cos(kappar);mtrRr[2] = -Math.Si n(phir) * Math.Cos(omegar);mtrRr[3] = Math.Cos(omegar) * Math.Si n(k appar);mtrRr[4] = Math.Cos(omegar) * Math.Cos(k appar);mtrRr[5] = -Math.Si n(omegar);mtrRr ⑹=Math.Si n(phir) * Math.Cos(k appar) + Math.Cos(phir) Math.Sin(omegar) * Math.Sin(kappar);mtrRr[7] = -Math.Sin(phir) * Math.Sin(kappar) + Math.Cos(phir) Math.Sin(omegar) * Math.Cos(kappar);mtrRr[8] = Math.Cos(phir) * Math.Cos(omegar);//求系数阵Adouble]] mtrA;mtrA = new double[576];for (int i = 0; i < 8; i++){mtrA[0 + 48 * i] = (mtrRI[0] * f + mtrRI[2] * xyz.GetElement(i, 0)) / (mtrRI[2] * (uvw.GetElement(0, i) - Xsl) + mtrRI[5] * (uvw.GetElement(1, i) - Ysl) + mtrRI[8] * (uvw.GetElement(2, i) - Zsl));mtrA[1 + 48 * i] = (mtrRI[3] * f + mtrRI[5] * xyz.GetElement(i, 0)) / (mtrRI[2] * (uvw.GetElement(0, i) - Xsl) + mtrRI[5] * (uvw.GetElement(1, i) - Ysl) + mtrRI[8] * (uvw.GetElement(2, i) - Zsl));mtrA[2 + 48 * i] = (mtrRI[6] * f + mtrRI[8] * xyz.GetElement(i, 0)) / (mtrRI[2] * (uvw.GetElement(0, i) - Xsl) + mtrRI[5] * (uvw.GetElement(1, i) - Ysl) + mtrRI[8] * (uvw.GetElement(2, i) - Zsl));mtrA[3 + 48 * i] = xyz.GetElement(i, 1) * Math.Sin(omegal)-(xyz.GetElement(i, 0) * (xyz.GetEIement(i, 0) * Math.Cos(kappal)- xyz.。