双向解析光束法.docx
- 文档编号:1232561
- 上传时间:2022-10-19
- 格式:DOCX
- 页数:13
- 大小:19.23KB
双向解析光束法.docx
《双向解析光束法.docx》由会员分享,可在线阅读,更多相关《双向解析光束法.docx(13页珍藏版)》请在冰豆网上搜索。
双向解析光束法
双向解析光束法
光束法程序有问题,在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;
//计算片左的加密点的物方坐标
Class1xyz=newClass1(8,3,strXY);
Class1xyzt=xyz.Transpose();
//phi,omega,kappa为零,故旋转矩阵为单位阵
//像空间辅助坐标系中的坐标
Class1UVW=xyzt.Multiply(H);
//求地面摄影测量坐标系中的坐标
for(inti=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=newdouble[24];
Console.WriteLine("请输入右片加密点的相片坐标");
for(inti=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(inti=0;i<24;i++)
{
strxy[i]=strxy[i]*0.000025;
}
//计算片左的加密点的物方坐标
Class1XYZ=newClass1(8,3,strxy);
Class1XYZt=XYZ.Transpose();
//phi,omega,kappa为零,故旋转矩阵为单位阵
//像空间辅助坐标系中的坐标
Class1uvw=XYZt.Multiply(H);
//求地面摄影测量坐标系中的坐标
for(inti=0;i<8;i++)
{
uvw.SetElement(0,i,(uvw.GetElement(0,i)+Xsr+UVW.GetElement(0,i))/2);
uvw.SetElement(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=newdouble[9];
mtrRl[0]=Math.Cos(phil)*Math.Cos(kappal)-Math.Sin(phil)*Math.Sin(omegal)*Math.Sin(kappal);
mtrRl[1]=-Math.Cos(phil)*Math.Sin(kappal)-Math.Sin(phil)*Math.Sin(omegal)*Math.Cos(kappal);
mtrRl[2]=-Math.Sin(phil)*Math.Cos(omegal);
mtrRl[3]=Math.Cos(omegal)*Math.Sin(kappal);
mtrRl[4]=Math.Cos(omegal)*Math.Cos(kappal);
mtrRl[5]=-Math.Sin(omegal);
mtrRl[6]=Math.Sin(phil)*Math.Cos(kappal)+Math.Cos(phil)*Math.Sin(omegal)*Math.Sin(kappal);
mtrRl[7]=-Math.Sin(phil)*Math.Sin(kappal)+Math.Cos(phil)*Math.Sin(omegal)*Math.Cos(kappal);
mtrRl[8]=Math.Cos(phil)*Math.Cos(omegal);
//求右片旋转矩阵Rr
double[]mtrRr;
mtrRr=newdouble[9];
mtrRr[0]=Math.Cos(phir)*Math.Cos(kappar)-Math.Sin(phir)*Math.Sin(omegar)*Math.Sin(kappar);
mtrRr[1]=-Math.Cos(phir)*Math.Sin(kappar)-Math.Sin(phir)*Math.Sin(omegar)*Math.Cos(kappar);
mtrRr[2]=-Math.Sin(phir)*Math.Cos(omegar);
mtrRr[3]=Math.Cos(omegar)*Math.Sin(kappar);
mtrRr[4]=Math.Cos(omegar)*Math.Cos(kappar);
mtrRr[5]=-Math.Sin(omegar);
mtrRr[6]=Math.Sin(phir)*Math.Cos(kappar)+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);
//求系数阵A
double[]mtrA;
mtrA=newdouble[576];
for(inti=0;i<8;i++)
{
mtrA[0+48*i]=(mtrRl[0]*f+mtrRl[2]*xyz.GetElement(i,0))/(mtrRl[2]*(uvw.GetElement(0,i)-Xsl)+mtrRl[5]*(uvw.GetElement(1,i)-Ysl)+mtrRl[8]*(uvw.GetElement(2,i)-Zsl));
mtrA[1+48*i]=(mtrRl[3]*f+mtrRl[5]*xyz.GetElement(i,0))/(mtrRl[2]*(uvw.GetElement(0,i)-Xsl)+mtrRl[5]*(uvw.GetElement(1,i)-Ysl)+mtrRl[8]*(uvw.GetElement(2,i)-Zsl));
mtrA[2+48*i]=(mtrRl[6]*f+mtrRl[8]*xyz.GetElement(i,0))/(mtrRl[2]*(uvw.GetElement(0,i)-Xsl)+mtrRl[5]*(uvw.GetElement(1,i)-Ysl)+mtrRl[8]*(uvw.GetElement(2,i)-Zsl));
mtrA[3+48*i]=xyz.GetElement(i,1)*Math.Sin(omegal)-(xyz.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 双向 解析 光束