任务已知f=153.24mm，m=10000，限差0.1’各点坐标点号像点坐标地面坐标x（mm）y（mm）X（m）Y（m）Z（m）1 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31求近似垂直摄影情况下后方交会解设计任务1、确定未知数的初始值：Φ0 =ω0 =К0 = 0 , 内方位元素，，f=153.24mm。

； = 38437m ；= 27963.16m2、计算旋转矩阵R利用角元素的近似值计算方向余弦值，组成R阵根据《摄影测量学》P32中的公式（3-9），初步计算R阵R[0][0]=cos(Φ)*cos(K)-sin(Φ)*sin(W)*sin(K);R[0][1]=-cos(Φ)*sin(K)-sin(Φ)*sin(W)*cos(K);R[0][2]=-sin(Φ)*cos(W);R[1][0]=cos(W)*sin(K);R[1][1]=cos(W)*cos(K);R[1][2]=-sin(W);R[2][0]=sin(Φ)*cos(K)+cos(Φ)*sin(W)*sin(K);R[2][1]=-sin(Φ)*sin(K)+cos(Φ)*sin(W)*cos(K);R[2][2]=cos(Φ)*cos(W);得初始R阵3、逐点计算近似值（x），（y）：带入《摄影测量学》P61的公式（5-1）；得4、组成误差方程式：按（5-8）；（5-9b）、（5-4）式逐点计算误差方程式的系数和常数项根据Lx=x-(x)；Ly=y-(y)得解得A阵为根据公式解除X阵，检验Φ1 、ω1 、К1 并不在限差0.1’之内;利用解得的Φ1 、ω1 、К1，带入2中循环；知道解出合格结果！5、结果，通过C++程序解得结果6、源程序代码#include <iostream>#include <fstream>#include <cmath>#include <string>#include <iomanip>using namespace std;const int n=6;void inverse (double c[n][n]);template<typename T1,typename T2>void transpose (T1*mat1,T2*mat2,int a,int b);template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);template<typename T>void input (T*mat,int a,int b);template<typename T>void 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 B[4][5]={0.0},R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];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,L[8][1],A[8][6];//确定未知数的出始值for(int i=0;i<4;i++){Xs=Xs+B[i][2];Ys=Ys+B[i][3];Zs=Zs+B[i][4];}Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk;//cout<<Xs<<" "<<Ys<<" "<<Zs<<endl;//int f=0;do//迭代计算{f++;//组成旋转矩阵R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K);R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K);R[0][2]=-sin(Q)*cos(W);R[1][0]=cos(W)*sin(K);R[1][1]=cos(W)*cos(K);R[1][2]=-sin(W);R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K);R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);R[2][2]=cos(Q)*cos(W);for(int t=0;t<3;t++){for(int j=0;j<3;j++)cout<<" "<<R[t][j];cout<<endl;}cout<<endl;//计算系数阵和常数项for(int i=0,k=0,j=0;i<=3;i++,k++,j++){X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs);Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs);Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs);//L[j][0]=B[i][0]-(x0-fk*X/Z); //l=x-(x)//插入cout<<L[j][0]<<" ";//L[j+1][0]=B[i][1]-(y0-fk*Y/Z); //l=y-(y)//插入cout<<L[j+1][0]<<endl;//j++;A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z;A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z;A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0) *sin(K))/fk+fk*cos(K))*cos(W);A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/f k;A[k][5]=B[i][1]-y0;A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z;A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z;A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][1]-y0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin(K))*cos(W);A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K)) /fk;A[k+1][5]=-(B[i][0]-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+XG[0][0]; Ys=Ys+XG[1][0]; Zs=Zs+XG[2][0];Q=Q+XG[3][0]; W=W+XG[4][0]; K=K+XG[5][0];//cout<<Q<<"**"<<W<<"**"<<K<<endl;//}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265 .0);cout<<"迭代次数为："<<f<<endl;//精度评定double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];multi(A,XG,AXG,8,6,1);for( i=0;i<8;i++) //计算改正数V[i][0]=AXG[i][0]-L[i][0];transpose (V,VT,1,8);multi(VT,V,VTV,1,8,1);m0=VTV[0][0]/2;for(i=0;i<6;i++)for(int j=0;j<6;j++)D[i][j]=m0*ATA[i][j];//屏幕输出误差方程系数阵、常数项、改正数output(A,"误差方程系数阵A为：",8,6);output(L,"常数项L为：",8,1);output(XG,"改正数为：",6,1);out("aim.txt",ios::app); //打开并添加aim.txt文件out(10);//以文件的形式输出像片外方位元素、旋转矩阵、方差阵outFile<<"一、像片的外方位元素为："<<endl<<endl;outFile<<setw(10)<<"Xs="<<Xs<<setw(10)<<"Ys="<<Ys<<setw(10)<<"Zs="<<Zs<<endl;outFile<<setw(20)<<"航向倾角为："<<Q<<setw(10)<<"旁向倾角为："<<W<<setw(10)<<"像片旋角为："<<K<<endl;outFile<<endl<<"二、旋转矩阵R为:"<<endl<<endl;for( i=0;i<3;i++){for(int j=0;j<3;j++)outFile<<setw(25)<<R[i][j]<<setw(25);outFile<<endl;}outFile<<endl;outFile<<setw(0)<<"三、精度评定结果为："<<endl;out(5);for(i=0;i<6;i++){for(int j=0;j<6;j++)outFile<<setw(14)<<D[i][j]<<setw(14);outFile<<endl;}out();return 0;}template<typename T1,typename T2>void transpose(T1*mat1,T2*mat2,int a,int b){ int i,j;for(i=0;i<b;i++)for(j=0;j<a;j++)mat2[j][i]=mat1[i][j];return;}template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c){ int i,j,k;for(i=0;i<a;i++){for(j=0;j<c;j++){result[i][j]=0;for(k=0;k<b;k++)result[i][j]+=mat1[i][k]*mat2[k][j];}}return;}template <typename T>void input (T*mat,int a,int b){ ifstream inFile;in("控制点坐标.txt");while(!in()){for (int i=0;i<a;i++)for(int j=0;j<b;j++)inFile>>mat[i][j];}in();return;}template<typename T>void output(T*mat,char*s,int a,int b){ cout<<setw(15)<<s<<endl;for(int i=0;i<a;i++){for(int j=0;j<b;j++)cout<<setw(13)<<mat[i][j];cout<<endl;}return;}void inverse(double c[n][n]){ int i,j,h,k;double p;double q[n][12];for(i=0;i<n;i++)//构造高斯矩阵for(j=0;j<n;j++)q[i][j]=c[i][j];for(i=0;i<n;i++)for(j=n;j<12;j++){if(i+6==j)q[i][j]=1;elseq[i][j]=0;}for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据for(i=k+1;i<n;i++){if(q[i][h]==0)continue;p=q[k][h]/q[i][h];for(j=0;j<12;j++){ q[i][j]*=p;q[i][j]-=q[k][j];}}for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据for(i=k-1;i>=0;i--){if(q[i][h]==0)continue;p=q[k][h]/q[i][h];for(j=0;j<12;j++){q[i][j]*=p;q[i][j]-=q[k][j];}}for(i=0;i<n;i++)//将对角线上数据化为1{ p=1.0/q[i][i];for(j=0;j<12;j++)q[i][j]*=p;}for(i=0;i<n;i++) //提取逆矩阵for(j=0;j<n;j++)c[i][j]=q[i][j+6];}7、感想这次摄影测量作业感觉非常难！一开始审题就出现了很大问题，读不懂，不知道要求什么东西。