单像空间摄影测量后方交会程序代码(vc++)
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug下。
1.单像空间后方交会原始数据：C语言程序：#include <stdio.h>#include <math.h>#include <iostream>double *readdata();void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa);void transpose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,double * mat2,double * result,int a,int b,int c); void inverse(double *a,int n)/*正定矩阵求逆*/{int i,j,k;for(k=0;k<n;k++){for(i=0;i<n;i++){if(i!=k)*(a+i*n+k)=-*(a+i*n+k)/(*(a+k*n+k));}*(a+k*n+k)=1/(*(a+k*n+k));for(i=0;i<n;i++){if(i!=k){for(j=0;j<n;j++){if(j!=k)*(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k);}}}for(j=0;j<n;j++){if(j!=k)*(a+k*n+j)*=*(a+k*n+k);}}}void transpose(double *m1,double *m2,int m,int n) //矩阵转置{ int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++)m2[j*m+i]=m1[i*n+j];return;}void multi(double *mat1,double *mat2,double * 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*c+j]=0;for(k=0;k<b;k++)result[i*c+j]+=mat1[i*b+k]*mat2[k*c+j];}}return;}double *readdata()FILE *fp;int i,j;int number;char datacatolog[100];//scanf("%s",datacatolog);if ((fp=fopen("控制点坐标.txt","r"))==NULL) {printf("读取数据出错!\n");return false;}fscanf(fp,"%d",&number);double *cordata=new double[number*5];for (i=0;i<number;i++){for (j=0;j<5;j++){fscanf(fp,"%lf",cordata+i*5+j);}}printf("控制点坐标数据读取成功！\n");return cordata;void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa){FILE *fp;char *file1="结算数据.txt";fp=fopen(file1,"w");fprintf(fp,"---------原始坐标数据为--------：\n");for (int i=0;i<hang;i++){for (int j=0;j<5;j++){fprintf(fp,"%7.4lf ",data[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------误差方程系数阵为：--------：\n");for (int i=0;i<hang*2;i++){for (int j=0;j<6;j++){fprintf(fp,"%7.4lf ",xishuarray[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n"); fprintf(fp,"---------法方程系数阵为：--------：\n"); for (int i=0;i<6;i++){for (int j=0;j<6;j++){fprintf(fp,"%7.5lf ",faxishu[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n"); fprintf(fp,"---------误差方程常数项为:--------：\n"); for (int i=0;i<hang*2;i++){fprintf(fp,"%lf ",l[i]);fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------迭代次数为:--------：\n");fprintf(fp,"%d\n",i);fprintf(fp,"--------------------------------\n");fprintf(fp,"-----------外方位元素为：---------\n");fprintf(fp," Xs= %lf, Ys=%lf, Zs=%lf\n",xs,ys,zs);fprintf(fp," fai= %lf, oumiga=%lf, kapa=%lf\n",fai,oumiga,kapa);fprintf(fp,"--------------------------------\n");fclose(fp);return;}void main(){int i,j;int ii,jj;int diedainumber=0;double x0=0.0,y0=0.0,f=0.0;double m=50000; //估算比例尺double fai=0,oumiga=0,kapa=0,Xs=0,Ys=0,Zs=0;double R[3][3]={0.0};double X=0.0,Y=0.0,Z=0.0,L[8][1]={0.0},A[8][6]={0.0};doublecorrect[6][1]={0.0},AT[6][8]={0.0},ATA[6][6]={0.0},ATL[6][1]={0.0};int row; //row用于存放坐标行数double *controlpoint;controlpoint=readdata();row=sizeof(controlpoint);for (i=0;i<row;i++){for (j=0;j<5;j++){printf("%3.3lf ",controlpoint[i*5+j]);}printf("\n");}/*----------内方位元素----------*/printf("请输入像片的内方位元素(mm)：\n");printf("x0=");x0/=1000.0;scanf("%lf",&x0); //double类型数据要用%lfprintf("y0=");y0/=1000.0;scanf("%lf",&y0);printf("f=");scanf("%lf",&f);f=f/1000.0;/*------------------------------*///-------确定未知数初始值------for(int i=0;i<row;i++){Xs=Xs+controlpoint[i*5+2];Ys=Ys+controlpoint[i*5+3];Zs=Zs+controlpoint[i*5+4];}Xs/=row;Ys/=row;Zs=Zs/row+m*f;//-----------------------------do{diedainumber++;//---------组成旋转矩阵--------R[0][0]=cos(fai)*cos(kapa)-sin(fai)*sin(oumiga)*sin(kapa);R[0][1]=-cos(fai)*sin(kapa)-sin(fai)*sin(oumiga)*cos(kapa);R[0][2]=-sin(fai)*cos(oumiga);R[1][0]=cos(oumiga)*sin(kapa);R[1][1]=cos(oumiga)*cos(kapa);R[1][2]=-sin(oumiga);R[2][0]=sin(fai)*cos(kapa)+cos(fai)*sin(oumiga)*sin(kapa);R[2][1]=-sin(fai)*sin(kapa)+cos(fai)*sin(oumiga)*cos(kapa);R[2][2]=cos(fai)*cos(oumiga);//-----------------------------//计算系数阵和常数项for(int i=0,k=0,j=0;i<=3;i++,k++,j++){X=R[0][0]*(controlpoint[i*5+2]-Xs)+R[1][0]*(controlpoint[i*5+3]-Y s)+R[2][0]*(controlpoint[i*5+4]-Zs);Y=R[0][1]*(controlpoint[i*5+2]-Xs)+R[1][1]*(controlpoint[i*5+3]-Y s)+R[2][1]*(controlpoint[i*5+4]-Zs);Z=R[0][2]*(controlpoint[i*5+2]-Xs)+R[1][2]*(controlpoint[i*5+3]-Y s)+R[2][2]*(controlpoint[i*5+4]-Zs);L[j][0]=controlpoint[i*5+0]-(x0-f*X/Z);L[j+1][0]=controlpoint[i*5+1]-(y0-f*Y/Z);j++;A[k][0]=(R[0][0]*f+R[0][2]*(controlpoint[i*5+0]-x0))/Z;A[k][1]=(R[1][0]*f+R[1][2]*(controlpoint[i*5+0]-x0))/Z;A[k][2]=(R[2][0]*f+R[2][2]*(controlpoint[i*5+0]-x0))/Z;A[k][3]=(controlpoint[i*5+1]-y0)*sin(oumiga)-((controlpoint[i*5+0 ]-x0)*((controlpoint[i*5+0]-x0)*cos(kapa)-(controlpoint[i*5+1]-y0)*si n(kapa))/f+f*cos(kapa))*cos(oumiga);A[k][4]=-f*sin(kapa)-(controlpoint[i*5+0]-x0)*((controlpoint[i*5+ 0]-x0)*sin(kapa)+(controlpoint[i*5+1]-y0)*cos(kapa))/f;A[k][5]=controlpoint[i*5+1]-y0;A[k+1][0]=(R[0][1]*f+R[0][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][1]=(R[1][1]*f+R[1][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][2]=(R[2][1]*f+R[2][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][3]=-(controlpoint[i*5+0]-x0)*sin(oumiga)-((controlpoint[i* 5+1]-y0)*((controlpoint[i*5+0]-x0)*cos(kapa)-(controlpoint[i*5+1]-y0) *sin(kapa))/f-f*sin(kapa))*cos(oumiga);A[k+1][4]=-f*cos(kapa)-(controlpoint[i*5+1]-y0)*((controlpoint[i* 5+0]-x0)*sin(kapa)+(controlpoint[i*5+1]-y0)*cos(kapa))/f;A[k+1][5]=-(controlpoint[i*5+0]-x0);k++;}transpose(A[0],AT[0],8,6);multi(AT[0],A[0],ATA[0],6,8,6);inverse(ATA[0],6);multi(AT[0],L[0],ATL[0],6,8,1);multi(ATA[0],ATL[0],correct[0],6,6,1);Xs=Xs+correct[0][0];Ys=Ys+correct[1][0];Zs=Zs+correct[2][0];fai=fai+correct[3][0];oumiga=oumiga+correct[4][0];kapa=kapa+correct[5][0];}while(correct[3][0]>=6.0/206265.0||correct[4][0]>=6.0/206265.0||c orrect[5][0]>=6.0/206265.0);printf("迭代次数为：%d\n",diedainumber);printf("---------误差方程系数为：--------\n");for (i=0;i<8;i++){for (j=0;j<6;j++){printf("%4.4lf ",A[i][j]);}printf("\n");}printf("--------------------------------\n");printf("求解得到的外方位元素为：\n");printf(" Xs= %lf\n",Xs);printf(" Ys= %lf\n",Ys);printf(" Zs= %lf\n",Zs);printf(" fai= %lf\n",fai);printf(" oumiga= %lf\n",oumiga);printf(" kapa= %lf\n",kapa);savedata(row,controlpoint,A[0],ATA[0],L[0],diedainumber,Xs,Ys,Zs, fai,oumiga,kapa);printf("-----------------解算结束！--------------\n");system("pause");}解算结果：2.后方交会-前方交会求解地面点坐标已知左右像片外方位元素，给出像点坐标：左像点坐标：右像点坐标：x（/m）y（/m）x（/m）y（/m）0.0053 0.0069 0.00482 0.0027C语言代码：#include <stdio.h>#include <iostream>#include <math.h>double *readdata();void savedata(int hang,double *data);double *readdata(){FILE *fp;int i,j,k;int number;char datacatolog[100];char leftdata[300];//scanf("%s",datacatolog);if ((fp=fopen("像点坐标数据.txt","r"))==NULL) {printf("读取数据出错!\n");system("pause");exit(0);}fscanf(fp,"%d",&number);double *c=new double[number*4];for (k=0;k<2;k++){fread(&leftdata,14,1,fp);for (i=0;i<number;i++){for (j=0;j<2;j++){fscanf(fp,"%lf",c+k*2+i*4+j);}}}fclose(fp);return c;}void savedata(int hang,double *data){FILE *fp;char *file1="地面点坐标数据.txt";fp=fopen(file1,"w");fprintf(fp,"---------像点对应地面点坐标为--------：\n");fprintf(fp,"\n");for (int i=0;i<hang;i++){fprintf(fp,"第%d点： ",i+1);for (int j=0;j<3;j++){fprintf(fp,"%7.4lf ",data[i*3+j]);}fprintf(fp,"\n\n");}fprintf(fp,"-----------------------------------------");fclose(fp);return;}void main(){double *imagepoint;int row;int i,j;imagepoint=readdata();row=sizeof(imagepoint);//--------------------------------------------double f=24;f/=1000;doublefai1=-0.0061,oumiga1=0.0327,kapa1=0.1711,Ys1=397367.171,Xs1=3445820.098,Zs1=148 6.212;doublefai2=0.0063,oumiga2=0.0178,kapa2=0.1489,Ys2=397367.234,Xs2=3445959.266,Zs2=1490 .096;// printf("请输入左像片的外方位元素：\n");//printf("Xs1= ");//scanf("%lf",&Xs1);//printf("Ys1= ");//scanf("%lf",&Ys1);//printf("Zs1= ");//scanf("%lf",&Zs1);//printf("fai1= ");//scanf("%lf",&fai1);//printf("oumiga1= ");//scanf("%lf",&oumiga1);//printf("kapa1= ");//scanf("%lf",&kapa1);//printf("请输入右像片的外方位元素：\n");//printf("Xs2= ");//scanf("%lf",&Xs2);//printf("Ys2= ");//scanf("%lf",&Ys2);//printf("Zs2= ");//scanf("%lf",&Zs2);//printf("fai2= ");//scanf("%lf",&fai2);//printf("oumiga2= ");//scanf("%lf",&oumiga2);//printf("kapa2= ");//scanf("%lf",&kapa2);double Bx=Xs2-Xs1,By=Ys2-Ys1,Bz=Zs2-Zs1;double N1=0,N2=0;double X1=0,Y1=0,Z1=0,X2=0,Y2=0,Z2=0;double R1[3][3]={0.0};double R2[3][3]={0.0};double GEOdata[4][3]={0.0};for (i=0;i<row;i++){//---------组成左影像旋转矩阵--------R1[0][0]=cos(fai1)*cos(kapa1)-sin(fai1)*sin(oumiga1)*sin(kapa1);R1[0][1]=-cos(fai1)*sin(kapa1)-sin(fai1)*sin(oumiga1)*cos(kapa1);R1[0][2]=-sin(fai1)*cos(oumiga1);R1[1][0]=cos(oumiga1)*sin(kapa1);R1[1][1]=cos(oumiga1)*cos(kapa1);R1[1][2]=-sin(oumiga1);R1[2][0]=sin(fai1)*cos(kapa1)+cos(fai1)*sin(oumiga1)*sin(kapa1);R1[2][1]=-sin(fai1)*sin(kapa1)+cos(fai1)*sin(oumiga1)*cos(kapa1);R1[2][2]=cos(fai1)*cos(oumiga1);//-----------------------------------//---------组成右影像旋转矩阵--------R2[0][0]=cos(fai2)*cos(kapa2)-sin(fai2)*sin(oumiga2)*sin(kapa2);R2[0][1]=-cos(fai2)*sin(kapa2)-sin(fai2)*sin(oumiga2)*cos(kapa2);R2[0][2]=-sin(fai2)*cos(oumiga2);R2[1][0]=cos(oumiga2)*sin(kapa2);R2[1][1]=cos(oumiga2)*cos(kapa2);R2[1][2]=-sin(oumiga2);R2[2][0]=sin(fai2)*cos(kapa2)+cos(fai2)*sin(oumiga2)*sin(kapa2);R2[2][1]=-sin(fai2)*sin(kapa2)+cos(fai2)*sin(oumiga2)*cos(kapa2);R2[2][2]=cos(fai2)*cos(oumiga2);//-------------像空辅系坐标-------------X1=R1[0][0]*imagepoint[i*4+0]+R1[0][1]*imagepoint[i*4+1]-R1[0][2]*f;Y1=R1[1][0]*imagepoint[i*4+0]+R1[1][1]*imagepoint[i*4+1]-R1[1][2]*f;Z1=R1[2][0]*imagepoint[i*4+0]+R1[2][1]*imagepoint[i*4+1]-R1[2][2]*f;X2=R2[0][0]*imagepoint[i*4+2]+R2[0][1]*imagepoint[i*4+3]-R2[0][2]*f;Y2=R2[1][0]*imagepoint[i*4+2]+R2[1][1]*imagepoint[i*4+3]-R2[1][2]*f;Z2=R2[2][0]*imagepoint[i*4+2]+R2[2][1]*imagepoint[i*4+3]-R2[2][2]*f;//--------------------------------------//------------点投影系数-------------N1=(Bx*Z2-Bz*X2)/(X1*Z2-Z1*X2);N2=(Bx*Z1-Bz*X1)/(X1*Z2-Z1*X2);//-----------------------------------//------------计算地面点坐标------------GEOdata[i][0]=Xs1+N1*X1;GEOdata[i][1]=Ys1+By+N2*Y2;GEOdata[i][2]=Zs1+N1*Z1;//--------------------------------------}//--------------------------------------------for (i=0;i<4;i++){printf("第%d个地面点坐标： ",i+1);for (j=0;j<3;j++){printf("%lf ",GEOdata[i][j]);}printf("\n\n");}savedata(row,GEOdata[0]);system("pause");}测试结果：3.单模型光束法严密平差缺少已知数据进行验证，因此如果有已知数据请代入已知数据进行验证。
摄影测量后方交会程序
摄影测量后方交会程序（c/c++）输入数据截图：结果截图：程序源代码（其中的矩阵求逆在前面已经有了，链接）：#include <stdio.h>#include <stdlib.h>#include <math.h>const double PRECISION=1e-5;typedef double DOUBLE[5];int InputData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f); int InverseMatrix(double *matrix,const int &row);int main(int argc, char* argv[]){DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if(Resection(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if (Data!=NULL){delete []Data;}printf("解算完毕...\n");do{printf("计算结果保存于\"结果.txt\"文件中\n""请选择操作（输入P打开结果数据，R打开原始数据，其它退出程序）："); fflush(stdin); //刷新输入流char order=getchar();if ('P'==order || 'p'==order){system("结果.txt");}else if ('R'==order || 'r'==order){system("data.txt");}elsebreak;system("cls");}while(1);system("PAUSE");return 0;}/***********************************************函数名：InputData*函数介绍：从文件（data.txt）中读取数据,*文件格式如下：*点数 m(未知写作0)* 内方位元素(f x0 y0)*编号 x y X Y Z*下面是一个实例：4 0153.24 0 01 -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*参数：(in/out)Num(点数)，*(in/out)Data(存放数据),m,f,x0,y0*返回值：int ，0成功，1文件打开失败，2控制点个*数不足，3文件格式错误*作者：vcrs*完成时间：09-10-4**********************************************/int InputData(int &Num, DOUBLE *&Data,double &m,double &f) {double x0,y0;FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fscanf(fp_input,"%d%lf",&Num,&m);if (Num<4){return 2;}fscanf(fp_input,"%lf%lf%lf",&f,&x0,&y0);f/=1000;if (m<0 || f<0){return 3;}Data=new DOUBLE[Num];double *temp= new double[Num-1];double scale=0;int i;for (i=0;i<Num;i++){//读取数据，忽略编号if(fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&Data[i][0],&Data[i][1],&Data[i][2],&Data[i][3],&Data[i][4])!=5){return 3;}//单位换算成mData[i][0]/=1000.0;Data[i][1]/=1000.0;}//如果m未知则归算其值if (0==m){for (i=0;i<Num-1;i++){temp[i]=(Data[i][2]-Data[i+1][2])/(Data[i][0]-Data[i+1][0])+ (Data[i][3]-Data[i+1][3])/(Data[i][1]-Data[i+1][1]);scale+=temp[i]/2.0;}m=scale/(Num-1);}fclose(fp_input);delete []temp;return 0;}/***********************************************函数名：MatrixMul*函数介绍：求两个矩阵的积,*参数：Jz1(第一个矩阵),row(第一个矩阵行数),*Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个*矩阵列数),(out)JgJz(存放结果矩阵)*返回值：void*作者：vcrs*完成时间：09-10-4**********************************************/void MatrixMul(double *Jz1,const int &row,double *Jz2,const int &line,const int &com,double *JgJz){for (int i=0;i<row;i++){for (int j=0;j<line;j++){double temp=0;for (int k=0;k<com;k++){temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j));}*(JgJz+i*line+j)=temp;}}}/***********************************************函数名：OutPut*函数介绍：向结果.txt文件输出数据*参数：Q协因数阵，m精度，m0单位权中误差，6个外*方位元素，旋转矩阵*返回值：int，0成功，1失败*作者：vcrs*完成时间：09-10-4**********************************************/int OutPut(const double *&Q,const double *&m,const double &m0, const double &Xs,const double &Ys,const double &Zs,const double &Phi,const double &Omega,const double &Kappa,const double *R){FILE *fp_out;if (!(fp_out=fopen("结果.txt","w"))){return 1;}FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"\n空间后方交会程序(C\\C++)\n遥感信息工程学院\n班级：" "00000\n学号：0000000\n姓名：vcrs\n\n");fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"已知数据：\n\n已知点数：");int num;double temp,x,y;fscanf(fp_input,"%d%lf",&num,&temp);fprintf(fp_out,"%d\n",num);fprintf(fp_out,"摄影比例尺（0表示其值位置）：");fprintf(fp_out,"%10.0lf\n",temp);fprintf(fp_out,"内方位元素(f x0 y0)：");fscanf(fp_input,"%lf%lf%lf",&temp,&x,&y);fprintf(fp_out,"%10lf\t%10lf\t%10lf\n",temp,x,y);for (int i=0;i<num;i++){double temp[5];fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&temp[0],&temp[1],&temp[2],&temp[3],&temp[4]);fprintf(fp_out,"%3d\t%10lf\t%10lf\t%10lf\t%10lf\t%10lf\n", i+1,temp[0],temp[1],temp[2],temp[3],temp[4]);}fclose(fp_input);fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"计算结果如下：\n\n外方位元素：\n");fprintf(fp_out,"\tXs=%10lf\n",Xs);fprintf(fp_out,"\tYs=%10lf\n",Ys);fprintf(fp_out,"\tZs=%10lf\n",Zs);fprintf(fp_out,"\tPhi=%10lf\n",Phi);fprintf(fp_out,"\tOmega=%10lf\n",Omega);fprintf(fp_out,"\tKappa=%10lf\n\n",Kappa);fprintf(fp_out,"旋转矩阵：\n");for (i=0;i<3;i++){fprintf(fp_out,"\t");for (int j=0;j<3;j++){fprintf(fp_out,"%10lf\t",*(R+i*3+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n单位权中误差：%10lf\n\n",m0);fprintf(fp_out,"协因数阵：\n");for (i=0;i<6;i++){fprintf(fp_out,"\t");for (int j=0;j<6;j++){fprintf(fp_out,"%20lf\t",*(Q+i*6+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n外方位元素精度：");for (i=0;i<6;i++){fprintf(fp_out,"%10lf\t",m[i]);}fprintf(fp_out,"\n");fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fclose(fp_out);return 0;}/***********************************************函数名：Resection*函数介绍：计算*参数：Num(点数)，Data(数据)，m,,f(焦距),x0,y0*返回值：int，0成功，其它失败*作者：vcrs*完成时间：09-10-4**********************************************/int Resection(const int &Num,const DOUBLE *&Data,const double &m, const double &f){double Xs=0,Ys=0,Zs=0;int i,j;//设置初始值for (i=0;i<Num;i++){Xs+=Data[i][2];Ys+=Data[i][3];}Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);double R[3][3]={0.0};double *L=new double[2*Num];typedef double Double6[6];Double6 *A=new Double6[2*Num];double *AT=new double[2*Num*6];double *ATA=new double[6*6];double *ATL=new double[6];double *Xg=new double[6];//迭代计算do{//旋转矩阵R[0][0]=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa); R[0][1]=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa); R[0][2]=-sin(Phi)*cos(Omega);R[1][0]=cos(Omega)*sin(Kappa);R[1][1]=cos(Omega)*cos(Kappa);R[1][2]=-sin(Omega);R[2][0]=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa); R[2][1]=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa); R[2][2]=cos(Phi)*cos(Omega);for (i=0;i<Num;i++){double X=R[0][0]*(Data[i][2]-Xs)+R[1][0]*(Data[i][3]-Ys)+R[2][0]*(Data[i][4]-Zs);double Y=R[0][1]*(Data[i][2]-Xs)+R[1][1]*(Data[i][3]-Ys)+ R[2][1]*(Data[i][4]-Zs);double Z=R[0][2]*(Data[i][2]-Xs)+R[1][2]*(Data[i][3]-Ys)+ R[2][2]*(Data[i][4]-Zs);double xxx,yyy;xxx=-f*X/Z;yyy=-f*Y/Z;//常数项L[2*i]=Data[i][0]-(-f*X/Z);L[2*i+1]=Data[i][1]-(-f*Y/Z);A[2*i][0]=(R[0][0]*f+R[0][2]*(xxx))/Z;A[2*i][1]=(R[1][0]*f+R[1][2]*(xxx))/Z;A[2*i][2]=(R[2][0]*f+R[2][2]*(xxx))/Z;A[2*i][3]=(yyy)*sin(Omega)-(((xxx)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))+f*cos(Kappa))*cos(Omega);A[2*i][4]=-f*sin(Kappa)-((xxx)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i][5]=(yyy);A[2*i+1][0]=(R[0][1]*f+R[0][2]*(yyy))/Z;A[2*i+1][1]=(R[1][1]*f+R[1][2]*(yyy))/Z;A[2*i+1][2]=(R[2][1]*f+R[2][2]*(yyy))/Z;A[2*i+1][3]=-(xxx)*sin(Omega)-(((yyy)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))-f*sin(Kappa))*cos(Omega);A[2*i+1][4]=-f*cos(Kappa)-((yyy)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i+1][5]=-(xxx);}//求矩阵A的转置矩阵ATfor (i=0;i<2*Num;i++){for (j=0;j<6;j++){*(AT+j*2*Num+i)=A[i][j];}}//求ATAMatrixMul(AT,6,&A[0][0],6,2*Num,ATA);if(InverseMatrix(ATA,6))return 1;MatrixMul(AT,6,L,1,2*Num,ATL);MatrixMul(ATA,6,ATL,1,6,Xg);Xs+=Xg[0];Ys+=Xg[1];Zs+=Xg[2];Phi+=Xg[3];Omega+=Xg[4];Kappa+=Xg[5];} while(fabs(Xg[0])>=PRECISION ||fabs(Xg[1])>=PRECISION || fabs(Xg[2])>=PRECISION ||fabs(Xg[3])>=PRECISION ||fabs(Xg[4])>=PRECISION || (Xg[5])>=PRECISION);//注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值，//由于变换很小忽略double *Q=ATA;double *V=new double[2*Num];MatrixMul(&A[0][0],2*Num,Xg,1,6,V);double VTV=0;for(i=0;i<2*Num;i++){V[i]-=L[i];VTV+=V[i]*V[i];}double m0=sqrt(VTV/(2*Num-6));double *mm=new double[6];for (i=0;i<6;i++){mm[i]=sqrt(*(Q+i*6+i))*m0;}OutPut(Q,mm,m0,Xs,Ys,Zs,Phi,Omega,Kappa,&R[0][0]);delete []L;delete []A;delete []AT;delete []ATA;delete []ATL;delete []Xg;delete []mm;delete []V;return 0;}void swap(double &a,double &b){double temp=a;a=b;b=temp;}/********************************************** *函数名：InverseMatrix*函数介绍：求矩阵的逆（高斯-约当法）*输入参数：(in/out)matrix（矩阵首地址），*（in）row（矩阵阶数）*输出参数：matrix（原矩阵的逆矩阵）*返回值：int ，0成功，1失败*调用函数：swap(double&,double&)*作者：vcrs*完成时间：09-10-4**********************************************/ int InverseMatrix(double *matrix,const int &row) {double *m=new double[row*row];double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;i<row;i++){for (j=0;j<row;j++){*pt=*ptemp;ptemp++;pt++;}}int k;int *is=new int[row],*js=new int[row];for (k=0;k<row;k++){double max=0;//全选主元//寻找最大元素for (i=k;i<row;i++){for (j=k;j<row;j++){if (fabs(*(m+i*row+j))>max){max=*(m+i*row+j);is[k]=i;js[k]=j;}}}if (0 == max){return 1;}//行交换if (is[k]!=k){for (i=0;i<row;i++){swap(*(m+k*row+i),*(m+is[k]*row+i)); }}//列交换if (js[k]!=k){for (i=0;i<row;i++){swap(*(m+i*row+k),*(m+i*row+js[k])); }}*(m+k*row+k)=1/(*(m+k*row+k));for (j=0;j<row;j++){if (j!=k){*(m+k*row+j)*=*((m+k*row+k));}for (i=0;i<row;i++){if (i!=k){for (j=0;j<row;j++){if(j!=k){*(m+i*row+j)-=*(m+i*row+k)**(m+k*row+j); }}}}for (i=0;i<row;i++){if(i!=k){*(m+i*row+k)*=-(*(m+k*row+k));}}}int r;//恢复行列for (r=row-1;r>=0;r--){if (js[r]!=r){for (j=0;j<row;j++){swap(*(m+r*row+j),*(m+js[r]*row+j));}if (is[r]!=r){for (i=0;i<row;i++){swap(*(m+i*row+r),*(m+i*row+is[r])); }}}ptemp=matrix;pt=m;for (i=0;i<row;i++){for (j=0;j<row;j++){*ptemp=*pt;ptemp++;pt++;}}delete []is;delete []js;delete []m;return 0;}。
摄影测量学单像空间后方交会编程实习报告(精品资料).doc
【最新整理，下载后即可编辑】摄影测量学单像空间后方交会编程实习报告班级：130x姓名：xx学号：2013302590xxx指导老师：李欣一、实习目的通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。
了解摄影测量平差的基本过程，掌握空间后方交会的定义和实现算法。
二、实习内容根据学习的单像空间后方交会的知识，用程序设计语言（C++或C语言）编写一个完整的单像空间后方交会程序，通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
三、实习数据已知航摄仪的内方位元素：fk =153.24mm，x=y=0，摄影比例尺为1:15000；4个地面控制点的地面坐标及其对应像点的像片坐标：四、实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。
因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。
可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。
单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。
五、 实习流程1. 获取已知数据。
从摄影资料中查取影像比例尺1/m ，平均摄影距离（航空摄影的航高、内方位元素x 0，y 0，f ；获取控制点的空间坐标X t ，Y t ，Z t 。
2. 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
摄影测量后方交会VC实现代码及实习报告
摄影测量课间实习单片空间后方交会班级：09031姓名：吴煜晖学号：20093025901232011-10-8一、实习原理单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素XS ，YS，ZS，ψ，ω，κ。
由于空间后方交会所采用的数学模型共线方程是非线性函数，为了便于外方位元素的解求，须首先对共线方程进行线性化，然后进行后方交会，最后在精度评定。
二、实习过程1、实习所用数据本次实习数据采用课本P.44 27题所给数据。
如下图：2、程序流程图及界面设计本程序程序框图如下：本人采用Visual C++6.0编此程序，利用MFC来设计程序。
主程序页面设计如下：子窗口（即进行计算后所得结果页面）设计如下：3、程序代码本程序代码较多，在此讲部分重要代码列出，其余代码参见程序源代码。
对话框类头文件内声明的类成员及函数（后来增加的）：public:void houfangjiaohui();double fi,w,k; //影像外方位角double Xs0,Ys0,Zs0; //后方交会所求解double a[3],b[3],c[3];double mx[6],m0;void inv(double *a,int n); /*正定矩阵求逆*/void transpose(double *m1, double *m2, int m, int n); //矩阵转置void mult(double *m1, double *m2, double *result, int i_1, int j_12, int j_2); //矩阵相乘部分核心代码，在此讲后方交会计算代码给出：void CFinalworkDlg::houfangjiaohui(){double t_fk;fi=w=k=0.0;int i; //中间变量double t_x[4],t_y[4],t_X[4],t_Y[4],t_Z[4];//将单位换成米t_x[0]=m_x1/1000.0; t_x[1]=m_x2/1000.0; t_x[2]=m_x3/1000.0;t_x[3]=m_x4/1000.0;t_y[0]=m_y1/1000.0; t_y[1]=m_y2/1000.0; t_y[2]=m_y3/1000.0;t_y[3]=m_y4/1000.0;t_X[0]=m_X1; t_X[1]=m_X2; t_X[2]=m_X3; t_X[3]=m_X4;t_Y[0]=m_Y1; t_Y[1]=m_Y2; t_Y[2]=m_Y3; t_Y[3]=m_Y4;t_Z[0]=m_Z1; t_Z[1]=m_Z2; t_Z[2]=m_Z3; t_Z[3]=m_Z4;Xs0=Ys0=Zs0=0.0;for(i=0;i<4;i++){Xs0+=t_X[i];Ys0+=t_Y[i];}//确定未知数初始值Xs0/=4;Ys0/=4;t_fk=m_fk/1000.0;Zs0=t_fk*m_mk;double x00=m_x00;double y00=m_y00;double A[8*6];double AT[6*8];double ATA[6*6];double L[8];double ATL[6*1];double Xo[4],Yo[4],Zo[4],Xom,Yom,Zom;double Delta[6];while(1){a[0]=cos(fi)*cos(k)-sin(fi)*sin(w)*sin(k);a[1]=-cos(fi)*sin(k)-sin(fi)*sin(w)*cos(k);a[2]=-sin(fi)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(fi)*cos(k)+cos(fi)*sin(w)*sin(k);c[1]=-sin(fi)*sin(k)+cos(fi)*sin(w)*cos(k);c[2]=cos(fi)*cos(w);for(i=0;i<4;i++){Xom=a[0]*(t_X[i]-Xs0)+b[0]*(t_Y[i]-Ys0)+c[0]*(t_Z[i]-Zs0);Yom=a[1]*(t_X[i]-Xs0)+b[1]*(t_Y[i]-Ys0)+c[1]*(t_Z[i]-Zs0);Zom=a[2]*(t_X[i]-Xs0)+b[2]*(t_Y[i]-Ys0)+c[2]*(t_Z[i]-Zs0);Xo[i]=-t_fk*Xom/Zom;Yo[i]=-t_fk*Yom/Zom;Zo[i]=Zom;A[i*12+0]=(a[0]*t_fk+a[2]*(t_x[i]-x00))/Zo[i];A[i*12+1]=(b[0]*t_fk+b[2]*(t_x[i]-x00))/Zo[i];A[i*12+2]=(c[0]*t_fk+c[2]*(t_x[i]-x00))/Zo[i];A[i*12+3]=(t_y[i]-y00)*sin(w)-cos(w)*(t_fk*cos(k)+((t_x[i]-x00)*cos(k)-(t _y[i]-y00)*sin(k))*(t_x[i]-x00)/t_fk);A[i*12+4]=-t_fk*sin(k)-(t_x[i]-x00)*((t_x[i]-x00)*sin(k)+(t_y[i]-y00)*cos(k ))/t_fk;A[i*12+5]=t_y[i]-y00;A[i*12+6]=(a[1]*t_fk+a[2]*(t_y[i]-y00))/Zo[i];A[i*12+7]=(b[1]*t_fk+b[2]*(t_y[i]-y00))/Zo[i];A[i*12+8]=(c[1]*t_fk+c[2]*(t_y[i]-y00))/Zo[i];A[i*12+9]=-(t_x[i]-x00)*sin(w)-((t_y[i]-y00)*((t_x[i]-x00)*cos(k)-(t_y[i]-y0 0)*sin(k))/t_fk-t_fk*sin(k))*cos(w);A[i*12+10]=-t_fk*cos(k)-(t_y[i]-y00)*((t_x[i]-x00)*sin(k)+(t_y[i]-y00)*cos(k))/t_fk;A[i*12+11]=-(t_x[i]-x00);L[i*2+0]=(t_x[i]-x00-Xo[i]);L[i*2+1]=(t_y[i]-y00-Yo[i]);}transpose(A,AT,8,6);mult(AT,A,ATA,6,8,6);mult(AT,L,ATL,6,8,1);inv(ATA,6);mult(ATA,ATL,Delta,6,6,1);Xs0+=Delta[0];Ys0+=Delta[1];Zs0+=Delta[2];fi+=Delta[3];w+=Delta[4];k+=Delta[5];if((fabs(Delta[0])>0.00001||fabs(Delta[1])>0.00001||fabs(Delta[2])>0.00001|| fabs(Delta[3])>0.00001||fabs(Delta[4])>0.00001||fabs(Delta[5])>0.00001)==0) {double V[8];double s=0.0;mult(A,Delta,V,8,6,1);for(i=0;i<8;i++){V[i]-=L[i];s+=V[i]*V[i];}m0=sqrt(s/2.0); //单位权中误差for(i=0;i<6;i++) //精度评定{mx[i]=sqrt(ATA[7*i])*m0;}break;}}}三、实习结果程序主界面如下：在该界面中，我们可以手动输入各项数据，也可以清空所有数据以及将所有空恢复默认数据（即题目所给数据）。
摄影测量学单像空间后方交会编程实习报告材料
摄影测量学单像空间后方交会编程实习报告班级： 130x姓名： xx学号： 2013302590xxx指导老师：李欣一、实习目的通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。
了解摄影测量平差的基本过程，掌握空间后方交会的定义和实现算法。
二、实习内容根据学习的单像空间后方交会的知识，用程序设计语言（C++或C语言）编写一个完整的单像空间后方交会程序，通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
三、实习数据已知航摄仪的内方位元素：fk =153.24mm，x=y=0，摄影比例尺为1:15000；4个地面控制点的地面坐标及其对应像点的像片坐标：四、实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。
因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。
可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。
单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。
五、实习流程1.获取已知数据。
从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y，；获取控制点的空间坐标Xt，Yt，Zt。
2.量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
单像空间后方交会和双像解析空间后方 前方交会的算法程序实现
（2）空间前方交会数学模型：
空间前方交会的计算步骤：由已知的外方位角元素及像点的坐标，根据（a）式计算像空间 辅助坐标；
X1 x 1 X 1 x1 Y R y Y R y 1 1 1 1 1 1 Z1 f ， Z1 f
1.单像空间后方交会的算法程序实现：
（1）空间后方交会的基本原理：对于遥感影像,如何获取像片的外方位元素,一直
是摄影测量工作者探讨的问题,其方法有:利用雷达(Radar)、 全球定位系统(GPS)、 惯性导航 系统(I N S)以及星像摄影机来获取像片的外方位元素;也可以利用一定数量的地面控制点, 根据共线方程,反求像片的外方位元素,这种方法称为单像空间后方交会(如图 1 所示)。
(4)空间后方交会的程序框图：
（5）单像空间后方交会的算法程序实现：
已知条件 摄影机主距 f=153.24mm，x0=0，y0=0, 像片比例尺为 1:40000，有四对点的像点坐标与相应 的地面坐标如下表。
点号
像点坐标
地面坐标
x(mm) 1 2 3 4 -86.15 -53.40 -14.78 10.46
（1）
由式 2 计算地面点的地面摄影测量坐标。
X A X S 1 N1 X 1 Y Y N Y A S1 1 1 ZA Z S1 N1Z1 X S 2 N2 X 2 N Y Y S 2 2 2 ZS 2 N2Z2 X S 1 BX N 2 X 2 B N Y Y S1 Y 2 2 Z S1 BZ N2Z2 （2）
空间后方交会程序
一． 实验目的： 掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。
二． 仪器用具及已知数据文件： 计算机windows xp 系统，编程软件（VISUAL C++6.0），地面控制点在摄影测量坐标系中的坐标及其像点坐标文件shuju.txt 。
三． 实验内容：单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标根据共线方程反求影像的外方位元素。
数学模型：共线条件方程式： )(3)(3)(3)(1)(1)(1Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f x -+-+--+-+--= )(3)(3)(3)(2)(2)(2Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f y -+-+--+-+--= 求解过程： （1）获取已知数据。
从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影测量坐标。
（2）量测控制点的像点坐标并做系统改正。
（3）确定未知数的初始值。
在竖直摄影且地面控制点大致分布均匀的情况下，按如下方法确定初始值，即： n X X S ∑=0,n Y Y S ∑=0,n Z mf Z S ∑=0 φ =ω=κ=0 式中；m 为摄影比例尺分母；n 为控制点个数。
（4）用三个角元素的初始值，计算个方向余弦，组成旋转矩阵R 。
（5）逐点计算像点坐标的近似值。
利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x ）、（y ）。
（6）逐点计算误差方程式的系数和常数项，组成误差方程式。
（7）计算法方程的系数矩阵A A T 和常数项l A T ，组成法方程式。
（8）解法方程，求得外方位元素的改正数dXs ,S dY ，s dZ ，d φ,d ω,d κ。
（9）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。
摄影测量学单像空间后方交会编程实习报告
摄影测量学单像空间后方交会编程实习报告班级:130x姓名:xx学号:2013302590xxx指导老师:李欣一、实习目得通过对提供得数据进行计算,输出像片得外方位元素并评定精度。
深入理解单像空间后方交会得思想,体会在有多余观测情况下,用最小二乘平差方法编程实现解求影像外方位元素得过程。
通过尝试编程实现加强编程处理问题得能力与对实习内容得理解,通过对实验结果得分析,增强综合运用所学知识解决实际问题得能力。
了解摄影测量平差得基本过程,掌握空间后方交会得定义与实现算法。
二、实习内容根据学习得单像空间后方交会得知识,用程序设计语言(C++或C语言)编写一个完整得单像空间后方交会程序,通过对提供得数据进行计算,输出像片得外方位元素并评定精度。
三、实习数据已知航摄仪得内方位元素:fk =153、24mm,x=y=0,摄影比例尺为1:15000;4个地面控制点得地面坐标及其对应像点得像片坐标:2 -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四、实习原理如果我们知道每幅影像得6个外方位元素,就能确定被摄物体与航摄影像得关系。
因此,如何获取影像得外方位元素,一直就是摄影测量工作者所探讨得问题。
可采取得方法有:利用雷达、全球定位系统(GPS)、惯性导航系统(INS)以及星相摄影机来获取影像得外方位元素;也可以利用影像覆盖范围内一定数量得控制点得空间坐标与摄影坐标,根据共线条件方程,反求该影像得外方位元素,这种方法称为单幅影像得空间后方交会。
单像空间后方交会得基本思想就是:以单幅影像为基础,从该影像所覆盖地面范围内若干控制点得已知地面坐标与相应点得像坐标量测值出发,根据共线条件方程,解求该影像在航空摄影时刻得外方位元素Xs,Ys,Zs,ϕ,ω,κ。
摄影测量后方交会近景摄影测量DLT直接线性变换代码见介绍
近景摄影测量实习报告班级： 07021班学号： 0062姓名：方毅日期： 2010年12月1日第一部分摄影、像点量测以及DLT 和单像后方空间交会解算1. 实习的目的和意义数码相机摄影：为后续的摄影测量解析处理提供质量合格的数字影像，了解所用数码相机的特点及使用，学习获取数字影像的方法。
像点量测：量测所拍摄的高精度室内三维控制场中控制点的像点坐标，为后续摄影测量解析处理准备计算数据。
直接线性变换（DLT ）与单像空间后方交会解算：加深理解近景摄影测量直接线性变换与单像空间后方交会的理论，学习准备数据和调试程序的方法。
2.实习原理 DLT 直接线性变换直接线性变换解法是建立像点坐标仪坐标和相应物点物方空间坐标直接的线性关系的算法。
它的基本关系式如下：1234910115678910110101l X l Y l Z l x l X l Y l Z l X l Y l Z l y l X l Y l Z +++⎧+=⎪+++⎪⎨+++⎪+=⎪+++⎩（）展开可得到以i l 为未知数的方程：1234910115678910110000000000l X l Y l Z l xl X xl Y xl Z x l X l Y l Z l yl X yl Y yl Z y +++++++++++=⎧⎨+++++++++++=⎩（） 当有n 个控制点时，即可列出2n 个方程式，写为矩阵的形式如下：111111111111111111111211100000000110000001n n n n n n n n n n n nnn nn nn n n X Y Z x X x Y x Z x l X Y Z y X y Y y Z y l X Y Z x X x Y x Z x l X Y Z y X y Y y Z y -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦⎣⎦（） 即AX L =。
摄影测量学单像空间后方交会程序设计作业
using System;using ;using ;namespace单像空间后方交会{class Program{static void Main string args{int x0, y0, i, j; double f, m;"请输入像片比例尺：";m = ;"请输入像片的内方位元素x0:";//均以毫米为单位x0 = ;"请输入像片的内方位元素y0:";y0 = ;"请输入摄影机主距f:";f = ;;//输入坐标数据double, zuobiao = new double4, 5;for i = 0; i < 4; i++{for j = 0; j < 5; j++{if j < 3{"请输入第{0}个点的第{1}个地面坐标：", i + 1, j + 1; zuobiaoi, j = ;}else{"请输入第{0}个点的第{1}个像点坐标：", i + 1, j - 2; zuobiaoi, j = ;}} ;}//归算像点坐标for i = 0; i < 4; i++{for j = 3; j < 5; j++{if j == 3zuobiaoi, j = zuobiaoi, j - x0;elsezuobiaoi, j = zuobiaoi, j - y0;}}//计算和确定初值double zs0 = m f, xs0 = 0, ys0 = 0;for i = 0; i < 4; i++{xs0 = xs0 + zuobiaoi, 0;ys0 = ys0 + zuobiaoi, 1;}xs0 = xs0 / 4;ys0 = ys0 / 4;//逐点计算误差方程系数double, xishu = new double8, 6;for i = 0; i < 8; i += 2{double x, y;x = zuobiaoi / 2, 3; y = zuobiaoi / 2, 4;xishui, 0 = xishui + 1, 1 = -1 / m; xishui, 1 = xishui + 1, 0 = 0; xishui, 2 = -x / m f; xishui, 3 = -f 1 + x x / f f;xishui, 4 = xishui + 1, 3 = -x y / f; xishui, 5 = y; xishui + 1, 2 = -y / m f; xishui + 1, 4 = -f 1 + y y / f f; xishui + 1, 5 = -x;}//计算逆阵double, dMatrix =matrixChematrixTransxishu, xishu;double, dReturn = ReverseMatrixdMatrix, 6;"逆矩阵为：";if dReturn = null{matrixOutdReturn;}//求解过程double phi0 = 0, omega0 = 0, kappa0 = 0; int q = 0;double, r = new double3, 3;double, jinsi = new double4, 2;double chazhi = new double8;double jieguo = new double6;double, zhong = matrixChedReturn, matrixTransxishu;do{ //计算旋转矩阵rr0, 0 = phi0 kappa0 - phi0 omega0 kappa0;r0, 1 = phi0 kappa0 - phi0 omega0 kappa0;r0, 2 = phi0 omega0;r1, 0 = omega0 kappa0;r1, 1 = omega0 kappa0;r1, 2 = omega0;r2, 0 = phi0 kappa0 + phi0 omega0 kappa0;r2, 1 = phi0 kappa0 + phi0 omega0 kappa0;r2, 2 = phi0 omega0;//计算x,y的近似值for i = 0; i < 4; i++{jinsii, 0 = -f r0, 0 zuobiaoi, 0 - xs0 + r1, 0 zuobiaoi, 1 - ys0 + r2, 0 zuobiaoi, 2 - zs0 / r0, 2 zuobiaoi, 0 - xs0 + r1, 2 zuobiaoi, 1 - ys0 + r2, 2 zuobiaoi, 2 - zs0; jinsii, 1 = -f r0, 1 zuobiaoi, 0 - xs0 + r1, 1 zuobiaoi, 1 - ys0 + r2, 1 zuobiaoi, 2 - zs0 / r0, 2 zuobiaoi, 0 - xs0 + r1, 2 zuobiaoi, 1 - ys0 + r2, 2 zuobiaoi, 2 - zs0; }for i = 0; i < 8; i += 2{chazhii = zuobiaoi / 2, 3 - jinsii / 2, 0;chazhii + 1 = zuobiaoi / 2, 4 - jinsii / 2, 1;}for i = 0; i < 0; i++{double k = 0;for j = 0; j < 1; j++{k = k + zhongi, j chazhij;}jieguoi = k;}//求新的近似值xs0 += jieguo0; ys0 += jieguo1; zs0 += jieguo2;phi0 += jieguo3; omega0 += jieguo4; kappa0 += jieguo5; q++;if q > 1000break;} while jieguo0 > || jieguo1 > || jieguo2 > ;"共进行了{0}次运算", q;"旋转矩阵为";matrixOutr;for i = 0; i < 0; i++{"第{0}个外方位元素为：{1}", i + 1, jieguoi;}}//矩阵转置public static double, matrixTrans double, X{double, A = X;double, C = new double1, 0;for int i = 0; i < 1; i++for int j = 0; j < 0; j++{Ci, j = Aj, i;}return C;}//矩阵输出public static void matrixOut double, X{double, C = X;for int i = 0; i < 0; i++{for int j = 0; j < 1; j++{" {0}", Ci, j;}"\n";}}//二维矩阵相乘public static double, matrixChe double, X, double, Y{int i, j, n; double m;double, C = X; double, D = Y;double, E = new double0, 0;for i = 0; i < 0; i++{for n = 0; n < 0; n++{m = 0;for j = 0; j < 1; j++{m = m + Ci, j Dj, n;}Ei, n = m;}}return E;}//计算行列式的值public static double MatrixValue double, MatrixList, int Level {double, dMatrix = new double Level, Level;for int i = 0; i < Level; i++for int j = 0; j < Level; j++dMatrixi, j = MatrixListi, j;double c, x;int k = 1;for int i = 0, j = 0; i < Level && j < Level; i++, j++{if dMatrixi, j == 0{int m = i;for ; dMatrixm, j == 0; m++ ;if m == Levelreturn 0;else{for int n = j; n < Level; n++{c = dMatrixi, n;dMatrixi, n = dMatrixm, n;dMatrixm, n = c;}k = -1;}}for int s = Level - 1; s > i; s--{x = dMatrixs, j;for int t = j; t < Level; t++dMatrixs, t -= dMatrixi, t x / dMatrixi, j;}}double sn = 1;for int i = 0; i < Level; i++{if dMatrixi, i = 0sn = dMatrixi, i;return 0;}return k sn;}//计算逆阵public static double, ReverseMatrix double, dMatrix, int Level {double dMatrixValue = MatrixValuedMatrix, Level;if dMatrixValue == 0 return null;double, dReverseMatrix = new double Level, 2 Level;double x, c;for int i = 0; i < Level; i++{for int j = 0; j < 2 Level; j++{if j < LeveldReverseMatrixi, j = dMatrixi, j;elsedReverseMatrixi, j = 0;}dReverseMatrixi, Level + i = 1;}for int i = 0, j = 0; i < Level && j < Level; i++, j++{if dReverseMatrixi, j == 0{int m = i;for ; dMatrixm, j == 0; m++ ;if m == Levelreturn null;else{for int n = j; n < 2 Level; n++dReverseMatrixi, n += dReverseMatrixm, n; }}x = dReverseMatrixi, j;if x = 1{for int n = j; n < 2 Level; n++if dReverseMatrixi, n = 0dReverseMatrixi, n /= x;}for int s = Level - 1; s > i; s--{x = dReverseMatrixs, j;for int t = j; t < 2 Level; t++dReverseMatrixs, t -= dReverseMatrixi, t x; }}for int i = Level - 2; i >= 0; i--{for int j = i + 1; j < Level; j++if dReverseMatrixi, j = 0c = dReverseMatrixi, j;for int n = j; n < 2 Level; n++dReverseMatrixi, n -= c dReverseMatrixj, n; }}double, dReturn = new double Level, Level;for int i = 0; i < Level; i++for int j = 0; j < Level; j++dReturni, j = dReverseMatrixi, j + Level;return dReturn;}}}。
单片空间后方交会C#源代码
单片空间后方交会C#源代码主方法：private void Cal_Click(object sender, EventArgs e) {string[] lines = RichText.Text.Split('\n');long m = lines.Length;m = m - 1;//真实数据行数double[] Coor_x = new double[m];//已知点x坐标double[] Coor_y = new double[m];//已知点x坐标double[] Coor_X = new double[m];//已知点X坐标double[] Coor_Y = new double[m];//已知点Y坐标double[] Coor_Z = new double[m];//已知点Z坐标///赋值for (int i = 0; i < m; i++){string[] FJstring = Regex.Split(lines[i+1], ",");Coor_x[i] = 0.001*(Convert.T oDouble(FJstring[0])); Coor_y[i] = 0.001 *( Convert.T oDouble(FJstring[1])); Coor_X[i] = Convert.ToDouble(FJstring[2]);Coor_Y[i] = Convert.ToDouble(FJstring[3]);Coor_Z[i] = Convert.ToDouble(FJstring[4]);}if (textBox_m.Text == ""){MessageBox.Show("请输入参数！");}if (textBox_m.Text != ""){double M = double.Parse(textBox_m.Text);//比例尺double f = 0.001 * (double.Parse(textBox_f.Text));//焦距double x0 = 0.001 * double.Parse(textBox_x0.Text);//内方位元素x0double y0 = 0.001 * double.Parse(textBox_y0.Text);//内方位元素y0double X0 = 0, Y0 = 0, Z0 = 0;//外方位坐标元素初始值double min = (double.Parse(textBox_k.Text));//焦距double angle1 = 0, angle2 = 0, angle3 = 0;//外方位角元素初始值for (int i = 0; i < m; i++){//累加X0 = Coor_X[i] + X0;Y0 = Coor_Y[i] + Y0;Z0 = Coor_Z[i] + Z0;}X0 = X0 / m;Y0 = Y0 / m;Z0 = Z0 / m + M * f;double[,] A = new double[2 * m, 6];//A，二维：（2 * m）*6 double[,] l = new double[2 * m, 1];//l，二维：（2 * m）*1double[,] V = new double[2 * m, 6];//V，二维：（2 * m）*6 double[,] XXX = new double[6, 1];//XXX，6*1 所求参数值，即外方位元素改正数///循环解求误差方程，以限差为标准结束循环Matrix ad = new Matrix();Cls_AL ae = new Cls_AL();int num = 0;//记录循环次数for (int i = 0; ; i++){//求A，lae.Cal_AL(X0, Y0, Z0, angle1, angle2, angle3, m, f, x0, y0, Coor_x, Coor_y, Coor_X, Coor_Y, Coor_Z, A, l);//求XXXXXX = ad.Cal_Chengfa(ad.Cal_Chengfa(ad.Cal_Qiuni(ad.Cal_Chengfa(ad .Cal_Zhuanzhi(A), A)), ad.Cal_Zhuanzhi(A)), l);///改正外方位元素X0 = X0 + XXX[0, 0];Y0 = Y0 + XXX[1, 0];Z0 = Z0 + XXX[2, 0];angle1 = angle1 + XXX[3, 0];angle2 = angle2 + XXX[4, 0];angle3 = angle3 + XXX[5, 0];num++;if ((Math.Abs(XXX[3, 0]) < min) && (Math.Abs(XXX[4, 0]) < min) && (Math.Abs(XXX[5, 0]) < min)){ break; }}///精度评定;V = ad.Cal_Jianfa(ad.Cal_Chengfa(A, XXX), l);double[,] c = new double[1, 1];//单位权中误差c = ad.Cal_Chengfa(ad.Cal_Zhuanzhi(V), V);double cc = Math.Sqrt(c[0, 0] / (2 * m - 6));double[,] AA = new double[6, 6];//外方位元素改正数的协因数阵QQ AA = ad.Cal_Qiuni(ad.Cal_Chengfa(ad.Cal_Zhuanzhi(A), A));RichText.Text = ("外方位元素为：" + '\n' + "Xs = " + Convert.T oString(X0) + '\n' + "Ys = " + Convert.ToString(Y0) + '\n' + "Zs = " + Convert.T oString(Z0) + '\n' + "angel_1 = " + Convert.T oString(angle1) + '\n' + "angel_2 = " + Convert.T oString(angle2) + '\n' + "angel_3 = " +Convert.T oString(angle3) + '\n' + "单位权中误差为：" + Convert.T oString(cc) + '\n' + "Xs精度为： " + Convert.T oString(cc * Math.Sqrt(AA[0, 0])) + '\n' + "Ys精度为：" + Convert.ToString(cc * Math.Sqrt(AA[1, 1])) + '\n' + "Zs精度为：" + Convert.ToString(cc * Math.Sqrt(AA[2, 2])) + '\n' + "迭代次数为：" + Convert.T oString(num));}}动态类库：using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace dll{////// Matrix类包含三个方法：矩阵乘法，矩阵转置，矩阵求逆///。
单片空间后方交会编程
矩阵相乘2_*1_2_*12_12_*1_2*1j i j j j i result m m ＝void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2){int i,j,k;for(i=0;i<i_1;i++)for(j=0;j<j_2;j++){result[i*j_2+j]=0.0;for(k=0;k<j_12;k++)result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];}return;}矩阵求逆n n n n m m *1*11＝ int invers_matrix (double *m1,int n ){int *is ,*js ;int i ,j ,k ,l ,u ,v ;double temp ,max_v ;is =(int *)malloc (n *sizeof (int ));js =(int *)malloc (n *sizeof (int ));if (is ==NULL ||js ==NULL ){printf ("out of memory!\n");return (0);}for (k =0;k <n ;k ++){max_v =0.0;for (i =k ;i <n ;i ++)for (j =k ;j <n ;j ++){temp =fabs (m1[i *n +j ]);if (temp >max_v ){max_v =temp ; is [k ]=i ; js [k ]=j ;}}if (max_v ==0.0){free (is ); free (js );printf("invers is not availble!\n");return(0);}if(is[k]!=k)for(j=0;j<n;j++){u=k*n+j; v=is[k]*n+j;temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}if(js[k]!=k)for(i=0;i<n;i++){u=i*n+k; v=i*n+js[k];temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}l=k*n+k;m1[l]=1.0/m1[l];for(j=0;j<n;j++)if(j!=k){u=k*n+j;m1[u]*=m1[l];}for(i=0;i<n;i++)if(i!=k)for(j=0;j<n;j++)if(j!=k){u=i*n+j;m1[u]-=m1[i*n+k]*m1[k*n+j];}for(i=0;i<n;i++)if(i!=k){u=i*n+k;m1[u]*=-m1[l];}}for(k=n-1;k>=0;k--){if(js[k]!=k)for(j=0;j<n;j++){u=k*n+j; v=js[k]*n+j;temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}if(is[k]!=k)for(i=0;i<n;i++){u=i*n+k; v=i*n+is[k];temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}}free (is ); free (js );return (1);}矩阵转置m n Tn m m m **21＝ void transpose (double *m1,double *m2,int m ,int n ){int i ,j ;for (i =0;i <m ;i ++)for (j =0;j <n ;j ++)m2[j *m +i ]=m1[i *n +j ];return ;}提交验收成果要求打印：1．程序实习报告（理解后方交会原理）；2 . 程序流程图；3．打印源程序代码（C 和C++独立完成，附带程序注释说明）；4．计算结果输出（包括： 原始计算数据， 计算的外方位元素并对其评定精度；。
单向空间后方交会matlab编程
%读取包含点号的像点坐标、控制点坐标分别到矩阵IP、CP%IP=load(,f:1)；CP=load(,f:,)；%删除像点坐标、控制点坐标矩阵中的点号列％IP(:,1 冃];cp(：,i)=n；%像片大小％H 二4008;W=5344;IP1=IP(:J卜W/2;¥IP2=H/2.|P(:,2);IP=[IP1JP2];%焦距％fx=；%内方位元素％x0=;yo=；%畸变参数％kl=;k2=;pl=；p2=;%外方位元素初值％Xs=3000;Ys=-100;Zs=100;Phi=0;Omega=0;Kappa=0;Im=size(IF>l);2);AIP=zeros(m/L=zeros(2*m,l)；A=zeros(2*m,6);X=ones(6l);zl)>=l ^(64)>=while X(4,l)>=| |X(5f%旋转矩阵的系数％al=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa); a2=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa); a3=-sin(Phi)*cos(Omega);bl=cos(Omega)*sin(Kappa);b2=cos(Omega)*cos(Kappa);b3=・sin(Omega);cl=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa); c2=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(0mega)*cos(Kappa); c3=cos(Phi)*cos(Omega);for i=l:m%求像点坐标常数项％r=sqrt((IP(iz l)-xO)A2+(IP(i/2)-yO)A2);Dx=(IP(i/l)-xO)*(kl*(r A2)+k2*(r A4))+pl*(r A2+2*((IP(i/l)-xO)A2))+2*p2*(IP(i/l)-xO)*(IP(i/2)-yO);Dy=(IP(i/2)-yO)*(kl*(r A2)+k2*(r A4))+p2*(r A2+2*((IP(i/2)-yO)A2))+2*pl*(IP(i/l)-xO)*(IP(i/2)-yO);X-=al*(CP(i/l)-Xs)+bl*(CP(i/2)-Ys)+cl*(CP(i/3)-Zs);Y-=a2*(CP(i/l)-Xs)+b2*(CP(i/2)-Ys)+c2*(CP(i/3)-Zs);Z-=a3*(CP(i/l)-Xs)+b3*(CP(i/2)-Ys)+c3*(CP(i/3)-Zs);%求像点坐标的近似值％AIP(izl)=-fx*X JZ_+Dx+xO;AIP(i,2)=-fy*YVZ_+Dy+yO;%求误差方程式系数％all=(al*fx+a3#IP(i/l))/Z_;al2=(bl*fx+b3*IP(i/l))/Z」(al3=(cl*fx+c3*IP(i/l))/Z_;al4=sin(Omega)*IP(i/2)-(IP(i/l)*(IP(i/l)*cos(Kappa)-IP(i/2)*sin(Kappa))/fx+fx*cos(Kappa))*cos(Om ega)；al5=-fx*sin(Kappa)-IP(i/l)*(IP(i/l)*sin(Kappa)+IP(i/2)*cos(Kappa))/fx;al6=IP(i，2);a21=(a2*fy+a3*IP(i/2))/ZJa22=(b2*fy+b3*IP(iz2))/Zja23=(c2*fy+c3*IP(iz2))/Za24=-sin(Omega)*IP(i/l)-(IP(i/2)*(IP(i/l)*cos(Kappa)-IP(i/2)*sin(Kappa))/fy-fy*sin(Kappa))*cos{Om ega);a25=-fy*cos(Kappa)-IP(i/2)*(IP(i/l)*sin(Kappa)+IP(i/2)*cos(Kappa))/fy;a26=-IP(U); %求常数项％L(2*i-l/l)=IP(i,l)-AIP(i/l)；L(2*i,l)=IP(i/2)-AIP(i,2); %求误差方程式系数阵％A(2*i-l,l)=all;A(2*i-l,2)=al2;A(2*i-13)=al3;A(2*i-l,4)=al4;A(2*i-l,5)=al5;A(2*i-l,6)=al6; AA(2*i4)=a21;A(2*iz2)=a22;A(2*b3)=a23;A(2*i,4)=a24;A(2*b5)=a25;A(2*i,6)=a26;end%求改正数％X=pinv(A,*A)*A,*L;%求真值％Xs=Xs+X(l,l)；Ys=Ys+X(2,l);Zs=Zs+X(3,l);Phi 二Phi+X(4J);Omega=Omega +X(5,1);Kappa=Kappa +X(6,1)；end%求误差矩阵％V=A*X-L;%精度评泄％>mO=abs(sqrt(V,*V/(2*m-6))); Q二inv(A 忖A);mi=mO*sqrt(Q);%{fpri ntf('%.6f\n\A);%%fprintf(,\n,)；%fprintf('%.6f\n%.6f\n%.6f\n%.6f\ n%.6f\n%.6f\n,/Xs/Ys/Zs/Phi,Omega/Kappa);fprintfCXn1)；%fprintf(,%.6f\n\L);% %fprintf(,\n,);%%fprintf(,%.6f\n\V);% %fprintf('\n1);%fprintf('%・6f\n：mi);。
空间后方交会程序设计
单片空间后方交会程序设计专业：测绘工程班级：测绘101姓名：张雯学号：2010220200512 用C 或VC++语言实现单片后方交汇的计算。
以单幅影像为基础，从该影象所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素Xs ，Ys ，Zs,ф,ω,κ.共线条件方程如下：x-x0=-f*[a1(X-Xs)+b1(Y-Ys)+c1(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] y-y0=-f*[a2(X-Xs)+b2(Y-Ys)+c2(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] x,y 为像点的像平面坐标；x0,y0,f 为影像的外方位元素；Xs,Ys,Zs 为摄站点的物方空间坐标；X,Y,Z 为物方点的物方空间坐标；输入原始数据归算像点坐标x,,y 计算和确定初值X s 0, Y s 0, Z s 0, φ0,ω0,κ0 组成旋转矩阵R 计算（x ），（y ）和l x ，l y 逐点组成误差方程并法化 所有点完否？ 解法方程，求未知数改正数 计算改正后外方位元素未知数改正数<限差否？整理并输出计算结果正常输出迭代次数小于限差否？ 输出中间结果和出错信息4个地面控制点的地面坐标及其对应像点的像片坐标：结果输出：已知条件:像点坐标x,y:-53.4 82.21-14.78 -76.6310.46 64.43153.24 -86.15地点坐标Xa,Ya,Za:37631.1 31324.5 728.6939101 24935 2386.540426.5 30319.8 757.31-68.99 36589.4 25273.3摄影比例尺：1:50000内方位元素:x0=y0=0 f=153.24mm计算结果:旋转矩阵：0.997709 0.0675334 0.00398399-0.0675254 0.997715 -0.00211178-0.0041175 0.00183792 0.99999像点坐标位:(单位：mm)-86.15 -68.99-53.41 82.21-14.78 -76.6310.47 64.43外方位元素：Xs=39795.435 精度为:1.1254813Ys=27476.479 精度为:1.2437625Zs=7572.6929 精度为:0.48380521q=-0.0039840098 精度为:0.00018182003 w=0.0021117837 精度为:0.00015959235 k=-0.067576934 精度为:7.2440432e-005 迭代次数：4Press any key to continue。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
单像空间后方交会程序西南交通大学土木工程学院测绘工程张慧鑫学号：20030593 输入文件形式如下：C++源程序如下：#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=39689; //估算比例尺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); //从文件中读取控制点的影像坐标和地面坐标，存入数组B double 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;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 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[j+1][0]=B[i][1]-(y0-fk*Y/Z);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))/ fk;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];}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/2062 65.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);outFile.open("aim.txt",ios::app); //打开并添加aim.txt文件outFile.precision(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;outFile.precision(5);for(i=0;i<6;i++){for(int j=0;j<6;j++)outFile<<setw(14)<<D[i][j]<<setw(14);outFile<<endl;}outFile.close();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;inFile.open("控制点坐标.txt");while(!inFile.eof()){for (int i=0;i<a;i++)for(int j=0;j<b;j++)inFile>>mat[i][j];}inFile.close();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];}程序的结果输出如下：（包括文本输出结果和荧屏输出中间数据）。