单向后方交会实验报告C++

实验报告班级：测绘一班学号:日期：2016.5.5目录一、计算原理.........................................3二、算法流程.........................................4三、源程序...........................................5四、计算结果........................................13五、结果分析........................................13六、心得体会........................................13一、计算原理已知条件摄影机主距f=153.24mm，x0=0.01mm，y0=0.02mm,像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。以单像空间后方交会方法，求解该像片的外方位元素。二、算法流程（1）获取已知数据。从航摄资料中差取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影坐标。（2）量测控制点的像点坐标并作系统误差改正。（3）确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即nXXS0，nYYS0，ZnmfZS100000（4）用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵R点号像点坐标地面坐标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coscoscossincossinsinsinsincoscossinsincoscossincoscossincossinsinsincossinsinsincoscos321321321cccbbbaaa（5）逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标；带入共线方程式，逐点近似像点坐标的近似值（x）、（y）。（6）逐点计算误差方程式的系数和常数项，组成误差方程式。（7）计算法方程的系数矩阵AAT和常数项LAL，组成法方程式。（8）解法方程，求得外方位元素的改正数dddddd、、、、、SSSZYX。（9）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。KSKSKSKSKSKSKSKSKSdZZZdYYYdXXX111,,KKKKKKKKKddd111,,（10）将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。负责用新的近似值重复（4）-（9），直到满足要求为止。用共线方程进行空间后方交会的程序框如图所示。输入原始数据像点坐标计算，系统误差正确定外方位因素初始值组成旋转矩阵R逐点组成误差方程式并法化所有像点完否是解法方程，求外方位元素改正数是否否二．源程序//2014113214徐福辉#includeiostream#includefstream#includeiomanipusingnamespacestd;constintn=6;voidinverse(doublec[n][n]);templatetypenameT1,typenameT2voidtranspose(T1*mat1,T2*mat2,inta,intb);templatetypenameT1,typenameT2voidmulti(T1*mat1,T2*mat2,T2*result,inta,intb,intc);intmain(){doublex[4][2]={-0.08616,-0.06897,-0.05341,0.08223,-0.01479,-0.07661,0.01045,0.06445};doubleX[4][3]={36589.41,25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98,2386.50,40426.54,30319.81,757.31};inti,j,m=1;//n为迭代次数doubleX0[6]={0};//设定未知数(XS,YS,ZS,ψ，ω，κ)初始值doublef=0.15324;//摄影机主距f=153.24mmdoublea=1/40000.0;//像片比例尺为1:40000doubleR[3][3]={0};//初始化旋转矩阵Rdoubledayue_x[8]={0};//用于存放像点估计值计算改正后的外方位元素外方位元素改正数是否小于限差是输出计算成果，计算并结束结束并显示错误信息迭代次数小于ndoubleA[8][6]={0};//系数阵doubleAT[6][8]={0};//A的转置矩阵doubleL[8]={0};//存放常数项constdoublepi=3.1415926535897932;doubleAsum[6][6]={0};doublejieguo2[6]={0};doublejieguo1[6][8]={0};doublesumXYZ[3]={0};cout.precision(5);cout已知像点坐标为：\n;for(i=0;i4;i++)for(j=0;j2;j++){coutfixed;if(j==0){coutxi+1=setw(10)x[i][j];}else{coutyi+1=setw(10)x[i][j]endl;}}cout已知地面四个点的坐标为：\n;for(i=0;i4;i++)for(j=0;j3;j++){if(j==0){coutXi+1=X[i][j];}elseif(j==1){coutYi+1=X[i][j];}else{coutZi+1;cout=;coutX[i][j]endl;}}coutendl;for(j=0;j3;j++)for(i=0;i4;i++)sumXYZ[j]+=X[i][j];for(i=0;i2;i++)X0[i]=sumXYZ[i]/4;//X0,Y0初始化X0[i]=1/a*f+sumXYZ[2]/4.0;//Z0初始化do{R[0][0]=cos(X0[3])*cos(X0[5])-sin(X0[3])*sin(X0[4])*sin(X0[5]);R[0][1]=-cos(X0[3])*sin(X0[5])-sin(X0[3])*sin(X0[4])*cos(X0[5]);R[0][2]=-sin(X0[3])*cos(X0[4]);R[1][0]=cos(X0[4])*sin(X0[5]);R[1][1]=cos(X0[4])*cos(X0[5]);R[1][2]=-sin(X0[4]);R[2][0]=sin(X0[3])*cos(X0[5])+cos(X0[3])*sin(X0[4])*sin(X0[5]);R[2][1]=-sin(X0[3])*sin(X0[5])+cos(X0[3])*sin(X0[4])*cos(X0[5]);R[2][2]=cos(X0[3])*cos(X0[4]);//第一个像点的估计值，其坐标位于X[0][0],X[0][1],X[0][2]dayue_x[0]=-f*(R[0][0]*(X[0][0]-X0[0])+R[1][0]*(X[0][1]-X0[1])+R[2][0]*(X[0][2]-X0[2]))/(R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2]));dayue_x[1]=-f*(R[0][1]*(X[0][0]-X0[0])+R[1][1]*(X[0][1]-X0[1])+R[2][1]*(X[0][2]-X0[2]))/(R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2]));//第二个像点的估计值,其坐标位于X[1][0],X[1][1],X[1][2]dayue_x[2]=-f*(R[0][0]*(X[1][0]-X0[0])+R[1][0]*(X[1][1]-X0[1])+R[2][0]*(X[1][2]-X0[2]))/(R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2]));dayue_x[3]=-f*(R[0][1]*(X[1][0]-X0[0])+R[1][1]*(X[1][1]-X0[1])+R[2][1]*(X[1][2]-X0[2]))/(R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2]));//第三个像点的估计值,其坐标位于X[2][0],X[2][1],X[2][2]dayue_x[4]=-f*(R[0][0]*(X[2][0]-X0[0])+R[1][0]*(X[2][1]-X0[1])+R[2][0]*(X[2][2]-X0[2]))/(R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2]));dayue_x[5]=-f*(R[0][1]*(X[2][0]-X0[0])+R[1][1]*(X[2][1]-X0[1])+R[2][1]*(X[2][2]-X0[2]))/(R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2]));//第四个像点的估计值,其坐标位于X[3][0],X[3][1],X[3][2]dayue_x[6]=-f*(R[0][0]*(X[3][0]-X0[0])+R[1][0]*(X[3][1]-X0[1])+R[2][0]*(X[3][2]-X0[2]))/(R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2]));dayue_x[7]=-f*(R[0][1]*(X[3][0]-X0[0])+R[1][1]*(X[3][1]-X0[1])+R[2][1]*(X[3][2]-X0[2]))/(R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2]));for(i=0;i4;i++){//第i个像点估计值放在dayue_x[2*(i-1)]A[2*i][0]=(R[0][0]*f+R[0][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));A[2*i][1]=(R[1][0]*f+R[1][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));A[2*i][2]=(R[2][0]*f+R[2][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R