空间后方交会报告.doc

任务已知f=153.24mm，m=10000，限差0.1各点坐标点号像点坐标地面坐标x（mm）y（mm）X（m）Y（m）Z（m）1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31求近似垂直摄影情况下后方交会解设计任务1、 确定未知数的初始值：0 =0 =0 = 0 , 内方位元素x0=0.0，y0=0.0，f=153.24mm。Zs0=H+ 14i=14Zi=mf+14i=14Zi=3049.3 m ；Xs0= 14i=14Xi = 38437m ； Ys0= 14i=14Yi = 27963.16m2、 计算旋转矩阵R 利用角元素的近似值计算方向余弦值，组成R阵根据摄影测量学P32中的公式（3-9），初步计算R阵R00=cos()*cos(K)-sin()*sin(W)*sin(K); R01=-cos()*sin(K)-sin()*sin(W)*cos(K); R02=-sin()*cos(W); R10=cos(W)*sin(K); R11=cos(W)*cos(K); R12=-sin(W); R20=sin()*cos(K)+cos()*sin(W)*sin(K); R21=-sin()*sin(K)+cos()*sin(W)*cos(K); R22=cos()*cos(W);得初始R阵1000100013、 逐点计算近似值（x），（y）：带入摄影测量学P61的公式（5-1）；得(x1)1(y1)1(x2)1(y2)1(x3)1(y3)1(x4)1(y4)1-0.33147-0.48258-0.053220.221960.15351-0.700980.133020.157564、 组成误差方程式：按（5-8）；（5-9b）、（5-4）式逐点计算误差方程式的系数和常数项根据Lx=x-(x)；Ly=y-(y)得Lx1Ly1Lx2Ly2Lx3Ly3Lx4Ly40.245320.41359-0.00018-0.13975-0.168290.62347-0.12256-0.09313解得A阵为-0.000100.000056219-0.201678-0.038785-0.068990-0.00010.0000450209-0.00000387855-0.18430.08615-0.000100.0000348473-0.171848-0.28648-0.082210-0.0001-0.0000536479-0.0000028648-0.1973440.0534534-0.000100.000009645-0.154666-0.00739096-0.076630-0.00010.0000500065-0.000000739096-0.191560.01478-0.00010-0.00000682589-0.153954-0.004397920.064430-0.0001-0.0000420452-0.000000439792-0.18033-0.01046根据公式X=（ATA）-1ATL解除X阵，检验1 、1 、1 并不在限差0.1之内;利用解得的1 、1 、1， 带入2中循环；知道解出合格结果！5、 结果，通过C+程序解得结果6、 源程序代码#include #include #include #include #include using namespace std;const int n=6;void inverse (double cnn);templatevoid transpose (T1*mat1,T2*mat2,int a,int b);templatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);templatevoid input (T*mat,int a,int b);templatevoid output(T*mat,char*s,int a,int b);int main()ofstream outFile;cout.precision(5);double x0=0.0, y0=0.0; double fk=0.15324; /内方位元素double m=10000; /估算比例尺 double B45=0.0,R33,XG61,AT68,ATA66,ATL61;input (B,4,5); /从文件中读取控制点的影像坐标和地面坐标，存入数组Bdouble Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0; double X,Y,Z,L81,A86; /确定未知数的出始值for(int i=0;i4;i+)Xs=Xs+Bi2; Ys=Ys+Bi3; Zs=Zs+Bi4;Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk;/coutXs Ys Zsendl;/int f=0;do/迭代计算f+;/组成旋转矩阵 R00=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R01=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R02=-sin(Q)*cos(W); R10=cos(W)*sin(K); R11=cos(W)*cos(K); R12=-sin(W); R20=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K); R21=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K); R22=cos(Q)*cos(W); for(int t=0;t3;t+) for(int j=0;j3;j+)cout Rtj; coutendl; coutendl; /计算系数阵和常数项 for(int i=0,k=0,j=0;i=3;i+,k+,j+) X=R00*(Bi2-Xs)+R10*(Bi3-Ys)+R20*(Bi4-Zs);Y=R01*(Bi2-Xs)+R11*(Bi3-Ys)+R21*(Bi4-Zs);Z=R02*(Bi2-Xs)+R12*(Bi3-Ys)+R22*(Bi4-Zs); / Lj0=Bi0-(x0-fk*X/Z); /l=x-(x) /插入coutLj0 ;/Lj+10=Bi1-(y0-fk*Y/Z); /l=y-(y)/插入coutLj+10endl;/j+;Ak0=(R00*fk+R02*(Bi0-x0)/Z;Ak1=(R10*fk+R12*(Bi0-x0)/Z;Ak2=(R20*fk+R22*(Bi0-x0)/Z;Ak3=(Bi1-y0)*sin(W)-(Bi0-x0)*(Bi0-x0)*cos(K)-(Bi1-y0)*sin(K)/fk+fk*cos(K)*cos(W);Ak4=-fk*sin(K)-(Bi0-x0)*(Bi0-x0)*sin(K)+(Bi1-y0)*cos(K)/fk;Ak5=Bi1-y0;Ak+10=(R01*fk+R02*(Bi1-y0)/Z;Ak+11=(R11*fk+R12*(Bi1-y0)/Z;Ak+12=(R21*fk+R22*(Bi1-y0)/Z;Ak+13=-(Bi0-x0)*sin(W)-(Bi1-y0)*(Bi0-x0)*cos(K)-(Bi1-y0)*sin(K)/fk-fk*sin(K)*cos(W);Ak+14=-fk*cos(K)-(Bi1-y0)*(Bi0-x0)*sin(K)+(Bi1-y0)*cos(K)/fk;Ak+15=-(Bi0-x0);k+;transpose(A,AT,6,8);multi(AT,A,ATA,6,8,6);inverse(ATA);multi(AT,L,ATL,6,8,1);multi(ATA,ATL,XG,6,6,1);Xs=Xs+XG00; Ys=Ys+XG10; Zs=Zs+XG20;Q=Q+XG30; W=W+XG40; K=K+XG50;/coutQ*W*K=6.0/206265.0|XG40=6.0/206265.0|XG50=6.0/206265.0);cout迭代次数为：fendl;/精度评定double AXG81,V81,VT18,VTV11,m0,D66;multi(A,XG,AXG,8,6,1); for( i=0;i8;i+) /计算改正数Vi0=AXGi0-Li0; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV00/2;for(i=0;i6;i+)for(int j=0;j6;j+)Dij=m0*ATAij;/屏幕输出误差方程系数阵、常数项、改正数output(A,误差方程系数阵A为：,8,6);output(L,常数项L为：,8,1);output(XG,改正数为：,6,1);outFile.open(aim.txt,ios:app); /打开并添加aim.txt文件 outFile.precision(10);/以文件的形式输出像片外方位元素、旋转矩阵、方差阵outFile一、像片的外方位元素为：endlendl;outFilesetw(10)Xs=Xssetw(10)Ys=Yssetw(10)Zs=Zsendl;outFilesetw(20)航向倾角为：Qsetw(10)旁向倾角为：Wsetw(10)像片旋角为：Kendl;outFileendl二、旋转矩阵R为:endlendl;for( i=0;i3;i+)for(int j=0;j3;j+)outFilesetw(25)Rijsetw(25);outFileendl;outFileendl;outFilesetw(0)三、精度评定结果为：endl;outFile.precision(5);for(i=0;i6;i+)for(int j=0;j6;j+)outFilesetw(14)Dijsetw(14);outFileendl;outFile.close();return 0;templatevoid transpose(T1*mat1,T2*mat2,int a,int b) int i,j; for(i=0;ib;i+) for(j=0;ja;j+)mat2ji=mat1ij;return;templatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resultij=0;for(k=0;kb;k+)resultij+=mat1ik*mat2kj;return;template void input (T*mat,int a,int b) ifstream inFile;inFile.open(控制点坐标.txt);while(!inFile.eof()for (int i=0;ia;i+) for(int j=0;jmatij;inFile.close();return;templatevoid output(T*mat,char*s,int a,int b) coutsetw(15)sendl;for(int i=0;ia;i+)for(int j=0;jb;j+) coutsetw(13)matij;coutendl;return;void inverse(double cnn) int i,j,h,k; double p; double qn12; for(i=0;in;i+)/构造高斯矩阵 for(j=0;jn;j+) qij=cij; for(i=0;in;i+) for(j=n;j12;j+) if(i+6=j) qij=1; else qij=0; for(h=k=0;kn-1;k+,h+)/消去对角线以下的数据 for(i=k+1;in;i+) if(qih=0) continue; p=qkh/qih; for(j=0;j0;k-,h-) / 消去对角线以上的数据 for(i=k-1;i=0;i-) if(qih=0) continue; p=qkh/qih; for(j=0;j12;j+) qij*=p; qij-=qkj;for(i=0;in;i+)/将对角线上数据化为1 p=1.0/qii; for(j=0;j12;j+) qij*=p; for(i=0;in;i+) /提取逆矩阵 for(j=0;jn;j+) cij=qij+6;7、 感想这次摄影测量作业感觉非常难！一开始审题就出现了很大问题，读不懂，不知道要求什么东西。按照书上5-2节的有关内容，采用EXCEL软件求解（因为C+虽然考了85分，但是考试成绩和实际操作水品完全不成正比，而且当初讲的时候都是些皮毛，数组用起来比较别扭，更别提用C+中的数组求解逆矩阵了，还有迭代。）老师上课没讲过，平时也不爱研究。当我用EXCEL求出第一遍的X阵的时候，感觉自己都要崩溃了，真的好烦。后来就在百度上搜求逆矩阵的C+源程序，搜到很多，但是都用不了，没办法，走投无路的我又在百度上搜起了“摄影测量后方交会”的有关内容，没想到居然出现了类似的题目。而且不止一个，我怀着激动地心情下载了一份，打开源程序我看呆了，感觉真的要是自己编出来的，该多牛，但是可悲的是我居然不大看得懂。我和孙维亚同学在下载的源程序上面看了好长时间，对里面一些地方进行修改调试，以便看得明白些。直到写报的的时候，我敢说，这份程序我也只能看懂一半的地方，而且这一半也只是理解性，要是自己编写的话，估计没辙。、下载的这份源程序，首先令人佩服的地方在于，以往我们自己编写的C+的程序都是从dos界面里面出入的，而这份源程序，通过inFile.open(控制点坐标.txt);直接从这个“控