摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug 下。
1.单像空间后方交会
C语言程序：
#include
#include
#include
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 { for(i=0;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 { if(i!=k) { for(j=0;j { if(j!=k) *(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k); } } } for(j=0;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 for(j=0;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 {for(j=0;j {result[i*c+j]=0; for(k=0;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 { 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 { 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 { 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 { 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}; double correct[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 { 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类型数据要用%lf printf("y0="); y0/=1000.0; scanf("%lf",&y0); printf("f="); scanf("%lf",&f); f=f/1000.0; /*------------------------------*/ //-------确定未知数初始值------ for(int i=0;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]-Ys)+R[2][0]*(co ntrolpoint[i*5+4]-Zs); Y=R[0][1]*(controlpoint[i*5+2]-Xs)+R[1][1]*(controlpoint[i*5+3]-Ys)+R[2][1]*(con trolpoint[i*5+4]-Zs); Z=R[0][2]*(controlpoint[i*5+2]-Xs)+R[1][2]*(controlpoint[i*5+3]-Ys)+R[2][2]*(con trolpoint[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)*((control point[i*5+0]-x0)*cos(kapa)-(controlpoint[i*5+1]-y0)*sin(kapa))/f+f*cos(kapa))*cos(o umiga); 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)*((cont rolpoint[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(ka pa)+(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||correct[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,kap a); printf("-----------------解算结束！--------------\n"); system("pause"); } 解算结果： 2.后方交会-前方交会求解地面点坐标 已知左右像片外方位元素，给出像点坐标： 左像点坐标：右像点坐标： x（/m）y（/m）x（/m）y（/m） 0.0053 0.0069 0.00482 0.0027 0.00556 0.00229 0.00521 -0.00191 0.00815 0.00499 0.00773 0.00089 -0.00247 0.00156 -0.00279 -0.00292 C语言代码： #include #include #include 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 { 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 { 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; double fai1=-0.0061,oumiga1=0.0327,kapa1=0.1711,Ys1=397367.171,Xs1=3445820.098,Zs1=1486.212; double fai2=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 { //---------组成左影像旋转矩阵-------- 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语言代码： #include #include 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,double xs2,double ys2,double zs2,double fai2,double oumiga2,double kapa2); 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 { for(i=0;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 { if(i!=k) { for(j=0;j { if(j!=k) *(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k); } } } for(j=0;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 for(j=0;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 {for(j=0;j {result[i*c+j]=0; for(k=0;k result[i*c+j]+=mat1[i*b+k]*mat2[k*c+j]; } } return; } 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"); return false; } fscanf(fp,"%d",&number); double *c=new double[number*7]; for (k=0;k<3;k++) { fread(&leftdata,14,1,fp); for (i=0;i<4;i++) { if (k==2) { for (j=0;j<3;j++) { fscanf(fp,"%lf",c+k*2+i*7+j); } } else { for (j=0;j<2;j++) { fscanf(fp,"%lf",c+k*2+i*7+j); } } } } fclose(fp); return c; } 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,double xs2,double ys2,double zs2,double fai2,double oumiga2,double kapa2) { FILE *fp; char *file1="结算数据.txt"; fp=fopen(file1,"w"); fprintf(fp,"---------原始坐标数据为--------：\n"); for (int i=0;i { for (int j=0;j<7;j++) { fprintf(fp,"%7.4lf ",data[i*7+j]); } fprintf(fp,"\n"); } fprintf(fp,"--------------------------------\n"); fprintf(fp,"---------误差方程系数阵为：--------：\n"); for (int i=0;i { for (int j=0;j<12;j++) { fprintf(fp,"%4.3e ",xishuarray[i*12+j]); } fprintf(fp,"\n"); } fprintf(fp,"--------------------------------\n"); fprintf(fp,"---------法方程系数阵为：--------：\n"); for (int i=0;i<12;i++) { for (int j=0;j<12;j++) { fprintf(fp,"%4.3e ",faxishu[i*12+j]); } fprintf(fp,"\n"); } fprintf(fp,"--------------------------------\n"); fprintf(fp,"---------误差方程常数项为:--------：\n"); for (int i=0;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," Xs1= %lf, Ys1=%lf, Zs1=%lf\n",xs,ys,zs); fprintf(fp," fai1= %lf, oumiga1=%lf, kapa1=%lf\n",fai,oumiga,kapa); fprintf(fp,"----------------------------------------\n"); fprintf(fp,"-----------右像的外方位元素为：---------\n"); fprintf(fp," Xs2= %lf, Ys2=%lf, Zs2=%lf\n",xs2,ys2,zs2); fprintf(fp," fai2= %lf, oumiga2=%lf, kapa2=%lf\n",fai2,oumiga2,kapa2); fprintf(fp,"----------------------------------------\n"); fclose(fp); return; } void main() { double *controlpoint; int row; int i,j; controlpoint=readdata(); row=sizeof(controlpoint); double f=0.024; double x0left=0,y0left=0; double x0right=0,y0right=0; double fai1=0,oumiga1=0,kapa1=0,Ys1=397510,Xs1=3445853,Zs1=1455; double fai2=0,oumiga2=0,kapa2=0,Ys2=397513,Xs2=3445979,Zs2=1453.685; double H=830; printf("请输入摄影机焦距："); scanf("%f=lf",&f); f=f/1000; printf("请输入摄影时相对航高："); scanf("%f=lf",&H); 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); int diedainumber=0; //double Xunkonwen=0.0,Yunkonwen=0.0,Zunkonwen=0.0; //待定点的坐标初始化 double X=0.0,Y=0.0,Z=0.0,L[16][1]={0.0},A[16][12]={0.0}; //先求两张相片的外方位元素double xcorrect[12][1]={0.0},AT[12][16]={0.0},ATA[12][12]={0.0},ATL[12][1]={0.0}; double R1[3][3]={0.0}; double R2[3][3]={0.0}; //----------------仅计算左右像片的外方位元素，没有考虑位置点---------------- do { diedainumber++; //---------组成左影像旋转矩阵-------- 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); //----------------------------------- //计算系数阵和常数项 for(int i=0,k=0,j=0;i { X=R1[0][0]*(controlpoint[i*7+4]-Xs1)+R1[1][0]*(controlpoint[i*7+5]-Ys1)+R1[2][0]*(controlp oint[i*7+6]-Zs1); Y=R1[0][1]*(controlpoint[i*7+4]-Xs1)+R1[1][1]*(controlpoint[i*7+5]-Ys1)+R1[2][1]*(controlp oint[i*7+6]-Zs1); Z=R1[0][2]*(controlpoint[i*7+4]-Xs1)+R1[1][2]*(controlpoint[i*7+5]-Ys1)+R1[2][2]*(controlp oint[i*7+6]-Zs1); L[j][0]=controlpoint[i*7+0]+f*X/Z; L[j+1][0]=controlpoint[i*7+1]+f*Y/Z; j++; //----------左像的误差方程---------- A[k][0]=-(f/H)*cos(kapa1); A[k][1]=-(f/H)*sin(kapa1); A[k][2]=-(controlpoint[i*7+0])/H; A[k][3]=-(f+(controlpoint[i*7+0]*controlpoint[i*7+0])/f)*cos(kapa1)+(controlpoint[i*7+0]*co ntrolpoint[i*7+1])/f*sin(kapa1); A[k][4]=-(f+(controlpoint[i*7+0]*controlpoint[i*7+0])/f)*sin(kapa1)-(controlpoint[i*7+0]*co ntrolpoint[i*7+1])/f*cos(kapa1); A[k][5]=controlpoint[i*7+1]; //A[k][12]=f/H; //A[k][13]=0; //A[k][14]=(controlpoint[i*7+0])/H; A[k+1][0]=(f/H)*sin(kapa1); A[k+1][1]=-(f/H)*cos(kapa1); A[k+1][2]=-(controlpoint[i*7+1])/H; A[k+1][3]=(f+(controlpoint[i*7+0]*controlpoint[i*7+0])/f)*sin(kapa1)+(controlpoint[i*7+0]*c ontrolpoint[i*7+1])/f*cos(kapa1); A[k+1][4]=-(f+(controlpoint[i*7+0]*controlpoint[i*7+0])/f)*cos(kapa1)-(controlpoint[i*7+0]* controlpoint[i*7+1])/f*sin(kapa1); A[k+1][5]=-controlpoint[i*7+0]; //A[k+1][12]=0; //A[k+1][13]=f/H; //A[k+1][14]=(controlpoint[i*7+1])/H; //------------------------------------ k=k+1; } //---------组成右影像旋转矩阵-------- 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); //----------------------------------- for(int i=0,k=row*2,j=row*2;i { X=R2[0][0]*(controlpoint[i*7+4]-Xs2)+R1[1][0]*(controlpoint[i*7+5]-Ys2)+R1[2][0]*(controlp oint[i*7+6]-Zs2); Y=R2[0][1]*(controlpoint[i*7+4]-Xs2)+R1[1][1]*(controlpoint[i*7+5]-Ys2)+R1[2][1]*(controlp oint[i*7+6]-Zs2); Z=R2[0][2]*(controlpoint[i*7+4]-Xs2)+R1[1][2]*(controlpoint[i*7+5]-Ys2)+R1[2][2]*(controlp oint[i*7+6]-Zs2); L[j][0]=controlpoint[i*7+2]+f*X/Z; L[j+1][0]=controlpoint[i*7+3]+f*Y/Z; j++; //----------右像的误差方程---------- A[k][6]=-(f/H)*cos(kapa2); A[k][7]=-(f/H)*sin(kapa2); A[k][8]=-(controlpoint[i*7+2])/H; A[k][9]=-(f+(controlpoint[i*7+2]*controlpoint[i*7+2])/f)*cos(kapa2)+(controlpoint[i*7+2]*co ntrolpoint[i*7+3])/f*sin(kapa2); A[k][10]=-(f+(controlpoint[i*7+2]*controlpoint[i*7+2])/f)*sin(kapa2)-(controlpoint[i*7+2]*c ontrolpoint[i*7+3])/f*cos(kapa2); A[k][11]=controlpoint[i*7+3]; //A[k][12]=f/H; //A[k][13]=0; //A[k][14]=(controlpoint[i*7+2])/H; A[k+1][6]=(f/H)*sin(kapa2); A[k+1][7]=-(f/H)*cos(kapa2); A[k+1][8]=-(controlpoint[i*7+3])/H; A[k+1][9]=(f+(controlpoint[i*7+2]*controlpoint[i*7+2])/f)*sin(kapa2)+(controlpoint[i*7+2]*c ontrolpoint[i*7+3])/f*cos(kapa2); A[k+1][10]=-(f+(controlpoint[i*7+2]*controlpoint[i*7+2])/f)*cos(kapa2)-(controlpoint[i*7+2] *controlpoint[i*7+3])/f*sin(kapa2); A[k+1][11]=-controlpoint[i*7+2]; //A[k+1][12]=0; //A[k+1][13]=f/H; //A[k+1][14]=(controlpoint[i*7+3])/H; //------------------------------------ k=k+1; } transpose(A[0],AT[0],16,12); multi(AT[0],A[0],ATA[0],12,16,12); inverse(ATA[0],12); multi(AT[0],L[0],ATL[0],12,16,1); multi(ATA[0],ATL[0],xcorrect[0],12,12,1); Xs1=Xs1+xcorrect[0][0]; Ys1=Ys1+xcorrect[1][0]; Zs1=Zs1+xcorrect[2][0]; fai1=fai1+xcorrect[3][0]; oumiga1=oumiga1+xcorrect[4][0]; kapa1=kapa1+xcorrect[5][0]; Xs2=Xs2+xcorrect[6][0]; Ys2=Ys2+xcorrect[7][0]; Zs2=Zs2+xcorrect[8][0]; fai2=fai2+xcorrect[9][0]; oumiga2=oumiga2+xcorrect[10][0]; kapa2=kapa2+xcorrect[11][0]; //Xunkonwen=Xunkonwen+XG[12][0]; //Yunkonwen=Yunkonwen+XG[13][0]; //Zunkonwen=Zunkonwen+XG[14][0]; }while(xcorrect[3][0]>=6.0/206265.0||xcorrect[4][0]>=6.0/206265.0||xcorrect[5][0]>=6.0/20 6265.0||xcorrect[9][0]>=6.0/206265.0||xcorrect[10][0]>=6.0/206265.0||xcorrect[11][0]>=6.0/20 6265.0); savedata(row,controlpoint,A[0],ATA[0],L[0],diedainumber,Xs1,Ys1,Zs1,fai1,oumiga1,kapa1,Xs 2,Ys2,Zs2,fai2,oumiga2,kapa2); printf("迭代次数为：%d\n\n",diedainumber); printf("左像片的外方位元素为：\n"); printf("xs1=%lf,ys1=%lf,zs1=%lf\n\n",Xs1,Ys1,Zs1); printf("fai1=%lf,oumiga1=%lf,kapa1=%lf\n\n",fai1,oumiga1,kapa1); printf("右像片的外方位元素为：\n"); printf("xs2=%lf,ys2=%lf,zs2=%lf\n\n",Xs2,Ys2,Zs2); printf("fai2=%lf,oumiga2=%lf,kapa2=%lf\n",fai2,oumiga2,kapa2); system("pause"); } 测试结果： 空间后方交会的解算 一． 空间后方交会的目的 摄影测量主要利用摄影的方法获取地面的信息，主要是是点位信息，属性信息，因此要对此进行空间定位和建模，并首先确定模型的参数，这就是空间后方交会的目的，用以求出模型外方位元素。 二． 空间后方交会的原理 空间后方交会的原理是共线方程。 共线方程是依据相似三角形原理给出的，其形式如下 111333222333()()() ()()() ()()()()()()A S A S A S A S A S A S A S A S A S A S A S A S a X X b Y Y c Z Z x f a X X a Y Y a Z Z a X X b Y Y c Z Z y f a X X a Y Y a Z Z -+-+-=--+-+--+-+-=--+-+- 上式成为中心投影的构线方程， 我们可以根据几个已知点，来计算方程的参数，一般需要六个方程，或者要三个点，为提高精度，可存在多余观测，然后利用最小二乘求其最小二乘解。 将公式利用泰勒公式线性化，取至一次项，得到其系数矩阵A ；引入改正数（残差）V ，则可将其写成矩阵形式： V AX L =- 其中 111333222333[,]()()()()()()()()()()()()()()T x y A S A S A S x A S A S A S A S A S A S y A S A S A S L l l a X X b Y Y c Z Z l x x x f a X X a Y Y a Z Z a X X b Y Y c Z Z l y y y f a X X a Y Y a Z Z =-+-+-=-=+-+-+--+-+-=-=+-+-+- 则1()T T X A A A L -= X 为外方位元素的近似改正数， 由于采用泰勒展开取至一次项，为减少误差，要将的出的值作为近似值进行迭代，知道小于规定的误差 三． 空间后方交会解算过程 1. 已知条件 近似垂直摄影 摄影测量原理 单张影像后方交会实习 目录 一实习目的 (3) 二实习原理 (3) 1. 间接平差 (3) 2. 共线方程 (3) 3. 单向空间后方交会 (4) 三计算流程 (4) 1. 求解步骤 (4) 2.计算机框图 (4) 四程序实现 (5) 五结果分析 (6) 1.外方位元素 (6) 2.误差 (6) 3.旋转矩阵R (7) 六实习体会 (7) 1. 平台的选择 (7) 2.问题的解决 (7) 3.心得体会 (8) 七代码展示 (8) 一实习目的 为了增强同学们对后方交会公式的理解，培养同学们对迭代循环编程的熟悉感，本次摄影测量课间实习内容定为用C语言或其他程序编写单片空间后方交会程序，最终输出像点坐标、地面坐标、单位权中误差、外方位元素及其精度。 已知四对点的影像坐标和地面坐标如下。内方位元素fk=153.24mm，x0=y0=0。 本次实习，我使用了matlab2014进行后方交会程序实现。结果与参考答案一致，精度良好。 二实习原理 题干中有四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标，由此可列出8个误差方程，存在2个多余观测（n=2）。故可利用间接平差的最小二乘法则求解。 由于共线方程是非线性函数模型，为了方便计算，需要将其“线性化”。但如果仅取泰勒级数展开式的一次项，未知数的近似值改正是不精确的。因此必须采用迭代趋近法计算，直到外方位元素的改正值小于限差。 1.间接平差 间接平差为平差计算最常用的方法。在确定多个未知量的最或然值时，选择它们之间不存在任何条件关系的独立量作为未知量组成用未知量表达测量的函数关系、列出误差方程式，按最小二乘法原理求得未知量的最或然值的平差方法。 在一个间接平差问题中，当所选的独立参数X个数与必要观测值t个数相等时，可将每个观测值表达成这t个参数的函数，组成观测方程。 函数模型为：L = BX + d。 2.共线方程 共线方程是中心投影构像的数学基础，也是各种摄影测量处理方法的重要理论基础。 式中： x，y 为像点的像平面坐标； x0，y0，f 为影像的内方位元素； XS，YS，ZS 为摄站点的物方空间坐标； XA，YA，ZA 为物方点的物方空间坐标； ai，bi，ci （i = 1,2,3）为影像的3 个外方位角元素组成的9 个方向余弦。 单张相片后方交会 目录 ●作业任务 (3) ●解算原理 (3) ●具体过程 (4) ●算法描述及程序流程 (4) ●计算结果 (7) ●结果分析 (8) ●心得体会及建议 (8) ●参考文献 (9) 一，作业任务 已知摄影机主距f=153.24mm，四对点的像点坐标与相应地面坐标列入下表： 表1-1 计算近似垂直摄影情况下后方交会解。 二，解算原理 【关键词1】中心投影构像方程 在摄影测量学中，最重要的方程就是中心投影构像方程（图2-1）。这个方程 将地面点在地面摄影测量坐标系中的坐标（物方坐标）和地面点对应像点的像平 面坐标联系起来。在解析摄影测量与数字摄影测量中是极其有用的。在以后将要 学习到的双像摄影测量光束法、解析测图仪原理及数字影像纠正等都要用到该 式。 图2-1 在上述公式中：x和y分别为以像主点为原点的像点坐标，相应地面点坐标 为X,Y,Z，相片主距f以及外方位元素Xs，Ys，Zs，ψ，ω，κ。 而在此次作业中，就是已知四个地面控制点的坐标以及其对应的像点坐标， 通过间接平差原理来求解此张航片的外方位元素。 【关键词2】间接平差 在一个平差问题中，当所选的独立参数X的个数等于必要观测值t时，可将 每个观测值表达成这t个参数的函数，组成观测方程，然后依据最小二乘原理求 解，这种以观测方程为函数模型的平差方法，就是间接平差方法 间接平差的函数模型为： 随机模型为： 平差准则为：VtPV=min 【关键词3】单像空间后方交会 利用至少三个已知地面控制点的坐标A（Xa，Ya，Za）、B（Xb，Yb，Zb）、Z（Xc， Yc，Zc），与其影像上对应的三个像点的影像坐标a（xa，ya）、b（xb，yb）、c（xc，yc），根据共线方程，反求该像点的外方位元素Xs，Ys，Zs，ψ，ω，κ。这种解算方法是以单张像片为基础，亦称单像空间后方交会。 在此次作业中，就是已知四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标。由此可以列出8个误差方程，存在两个多余观测数，则n=2。故可利用间接平差里，依据最小二乘法则，进行求解。由于共线条件方程是非线性函数模型，为了便于计算，需把非线性函数表达式用泰勒公式展开成现行形式，即“线性化”。而又因为仅取泰勒级数展开式的一次项，未知数的近似值改正是粗略的，所以必须计算采用逐渐趋近法，解求过程需要反复趋近，直至改正值小于限差为止。 三，具体过程 1，获取已知点数据：从摄影资料中查取像片比例尺1/m，平均航高，内方元素x0，y0，f；从外业测量成果中，获取控制点的地面测量坐标Xt，Yt，Zt，并转换成摄影测量坐标X,Y,Z。 2，量测控制点的像点坐标：将控制点标刺在像片上，利用立体坐标量测仪量测控制点的像框坐标，并经像点坐标改正，得到像点坐标x，y。 3，确定未知数的初始值：在竖直摄影测量情况下，角元素的初始值为0，及ψ=ω=κ=0； 线元素中，Zso =m*f+(Z[0]+Z[1]+Z[2]+Z[3])/4，Xso，Yso的取值可用四个角点上制点坐标的平均值，即：Xso=(X[0]+X[1]+X[2]+X[3])/4;Yso=(Y[0]+Y[1]+Y[2]+Y[3])/4；4，计算旋转矩阵R：利用角元素的近似值计算方向余弦，组成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); 5，逐点计算像点坐标的近似值：利用未知数的近似值按共线方程计算控制点像点坐标的近似值（x）、（y）； 6，组成误差方程式：参照教材（5-8）式、（5-9b）式、（5-4）式逐点计算误差方程的系数阵和常数项。 7，组成法方程：计算法方程的系数矩阵与常数项。 8，解求外方位元素：根据法方程，按间接平差原理解求外方位元素改正值，并与相应的近似值求和，得到外方位元素的新的近似值。 9，检查计算是否收敛：将求得的外方位元素的改正值与规定的限差比较，小于限差则计算终止，否则用新的近似值重复第4至第8步骤计算，直至满足要求为止。 四，算法描述及程序流程。 算法描述（图4-1）： 单像空间后方交会和双像解析空间后方-前 方交会的算法程序实现 遥感科学与技术 摘要：如果已知每张像片的6个外方位元素，就能确定被摄物体与航摄像片的关系。因此，利用单像空间后方交会的方法，可以迅速的算出每张像片的6个外方位元素。而前方交会的计算，可以算出像片上点对应于地面点的三维坐标。基于这两点，利用计算机强大的运算能力，可以代替人脑快速的完成复杂的计算过程。 关键词：后方交会，前方交会，外方位元素，C++编程 0.引言： 单张像片空间后方交会是摄影测量基本问题之一,是由若干控制点及其相应像点坐标求解摄站参数(X S,Y S,ZS,ψ、ω、κ)。单像空间后方交会主要有三种方法:基于共线条件方程的平差解法、角锥法、基于直接线性变换的解法。而本文将介绍第一种方法，基于共线条件方程反求象片的外方位元素。 而空间前方交会先以单张像片为单位进行空间后方交会，分别求出两张像片的外方位元素，再根据待定点的一对像点坐标，用空间前方交会的方法求解待定点的地面坐标。可以说，这种求解地面点的坐标的方法是以单张像片空间后方交会为基础的，因此，单张像片空间后方交会成为解决这两个问题以及算法程序实现的关键。 1.单像空间后方交会的算法程序实现： （1）空间后方交会的基本原理：对于遥感影像,如何获取像片的外方位元素,一直是摄影测量工作者探讨的问题,其方法有:利用雷达(Radar)、全球定位系统(GPS)、惯性导航系统(I N S)以及星像摄影机来获取像片的外方位元素;也可以利用一定数量的地面控制点,根据共线方程,反求像片的外方位元素,这种方法称为单像空间后方交会(如图1所示)。 图中,地面坐标X i、Yi、Zi和对应的像点坐标x i、yi是已知的,外方位元素XS、Y S、ZS(摄站点坐标),ψ、ω、κ(像片姿态角)是待求的。 （2）空间后方交会数学模型：空间后方交会的数学模型是共线方程, 即中心投影的构像方程： 式中X、Y、Z是地面某点在地面摄影测量坐标系中的坐标,x,y是该地面点在像片上的构像点的像片坐标,对 于空间后方交会而言它们是已知的,还有主距f是已知的。而9个方向余弦a 1,a 2,a3;b1,b 2,b 3;c 1,c2,c 3是未知的，具体表达式可以取 空间后方交会编程实习报告 一实习目的 用程序设计语言（Visual C++或者C语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。 二实习内容 利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。 三实习数据 已知航摄仪的内方位元素：f k =153.24mm，x =y =0.0mm，摄影比例尺为1:50000； 4个地面控制点的地面坐标及其对应像点的像片坐标： 四实习原理 如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。 单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，t，w，k。 五实习流程 （1）获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。 （2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。 （3）确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值： 摄影测量学 单像空间后方交会实习报告 、管路敷设技术通过管线不仅可以解决吊顶层配置不规范高中资料试卷问题，而且可保障各类管路习题到位。在管路敷设过程中，要加强看护关于管路高中资料试卷连接管口处理高中资料试卷弯扁度固定盒位置保护层防腐跨接地线弯曲半径标等，要求技术交底。管线敷设技术中包含线槽、管架等多项方式，为解决高中语文电气课件中管壁薄、接口不严等问题，合理利用管线敷设技术。线缆敷设原则：在分线盒处，当不同电压回路交叉时，应采用金属隔板进行隔开处理；同一线槽内强电回路须同时切断习题电源，线缆敷设完毕，要进行检查和检测处理。、电气课件中调试对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行 高中资料试卷调整试验；通电检查所有设备高中资料试卷相互作用与相互关系，根据生产工艺高中资料试卷要求，对电气设备进行空载与带负荷下高中资料试卷调控试验；对设备进行调整使其在正常工况下与过度工作下都可以正常工作；对于继电保护进行整核对定值，审核与校对图纸，编写复杂设备与装置高中资料试卷调试方案，编写重要设备高中资料试卷试验方案以及系统启动方案；对整套启动过程中高中资料试卷电气设备进行调试工作并且进行过关运行高中资料试卷技术指导。对于调试过程中高中资料试卷技术问题，作为调试人员，需要在事前掌握图纸资料、设备制造厂家出具高中资料试卷试验报告与相关技术资料，并且了解现场设备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况 ，然后根据规范与规程规定，制定设备调试高中资料试卷方案。 、电气设备调试高中资料试卷技术电力保护装置调试技术，电力保护高中资料试卷配置技术是指机组在进行继电保护高中资料试卷总体配置时，需要在最大限度内来确保机组高中资料试卷安全，并且尽可能地缩小故障高中资料试卷破坏范围，或者对某些异常高中资料试卷工况进行自动处理，尤其要避免错误高中资料试卷保护装置动作，并且拒绝动作，来避免不必要高中资料试卷突然停机。因此，电力高中资料试卷保护装置调试技术，要求电力保护装置做到准确灵活。对于差动保护装置高中资料试卷调试技术是指发电机一变压器组在发生内部故障时，需要进行外部电源高中资料试卷切除从而采用高中资料试卷主要保护装置。 摄影测量后方交会程序（c/c++） 输入数据截图： 结果截图： 程序源代码（其中的矩阵求逆在前面已经有了，链接）： #include } 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"); } else break; system("cls"); }while(1); system("PAUSE"); return 0; } /********************************************** *函数名：InputData *函数介绍：从文件（data.txt）中读取数据, *文件格式如下： *点数 m(未知写作0) * 内方位元素(f x0 y0) *编号 x y X Y Z *下面是一个实例： 4 0 153.24 0 0 1 -86.15 -68.99 36589.41 25273.3 2 2195.17 程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug 下。 1.单像空间后方交会 C语言程序： #include 摄影测量学空间后方交会实验报告测绘101徐斌摄影测量学实验报告 实验一、单像空间后方交会 学院: 建测学院 班级: 测绘101 姓名: 徐斌 学号: 26 一( 实验目的 1.深入了解单像空间后方交会的计算过程; 2.加强空间后方交会基本公式和误差方程式，法线方程式的记忆; 3.通过上机调试程序加强动手能力的培养。 二(实验原理 以单幅影像为基础，从该影像所覆盖地面范围内若干控制点和相应点的像坐标量测值出发，根据共线条件方程，求解该影像在航空摄影时刻的相片外方位元素。 三(实验内容 1.程序图框图 2.实验数据 (1)已知航摄仪内方位元素f,153.24mm，Xo,Yo,0。限差0.1秒 (2)已知4对点的影像坐标和地面坐标: 影像坐标地面坐标 x(mm) y(mm) X(m) Y(m) Z(m) 1 -86.15 -68.99 36589.41 25273.32 2195.17 2 -53.40 82.21 37631.08 31324.51 728.69 3 -14.78 -76.63 39100.97 24934.98 2386.50 4 10.46 64.43 40426.54 30319.81 757.31 3.实验程序 Form1.cs 程序 using System; using System.Collections.Generic; using https://www.sodocs.net/doc/2b11774748.html,ponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; using System.IO; namespace 后方交会1 { public partial class Form1 : Form { public Form1() { InitializeComponent(); } public double f,m, Xs, Ys, Zs, a1, a2, a3, b1, b2, b3, c1, c2, c3, q, w, k; public static int N,s; public double[] x = new double[4]; 空间后方交会程序 ————————————————————————————————作者：————————————————————————————————日期： ? 一． 实验目的: 掌握摄影测量空间后方交会的原理,利用计算机编程语言实现空间 后方交会外方位元素的解算。 二． 仪器用具及已知数据文件： 计算机wind ｏws xp 系统,编程软件(VI ＳＵA Ｌ C ＋+6．0)，地面控 制点在摄影测量坐标系中的坐标及其像点坐标文件shu ｊu.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)量测控制点的像点坐标并做系统改正。 (３)确定未知数的初始值。在竖直摄影且地面控制点大致分布均匀 的情况下,按如下方法确定初始值,即: n X X S ∑=0,n Y Y S ∑=0,n Z mf Z S ∑=0 φ =ω=κ=0 式中;ｍ为摄影比例尺分母;ｎ为控制点个数。 （４)用三个角元素的初始值,计算个方向余弦，组成旋转矩阵R 。 (5)逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标代入共 线方程式，逐点计算像点坐标的近似值（x ）、（y ）。 （6)逐点计算误差方程式的系数和常数项,组成误差方程式。 (7）计算法方程的系数矩阵A A T 和常数项l A T ,组成法方程式。 （8）解法方程，求得外方位元素的改正数dXs ,S dY ，s dZ ,d φ,ｄω，d κ。 （９）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。 { System; System.Collections.Generic; System.Linq; System.Text; namespace 单像空间后方交会 { class Program { static void Main( string [] args) for (j = 0; j < 5; j++) if (j < 3) "请输入第 {0} 个点的第 {1} 个地面坐标： ", i + 1, j + 1); double .Parse( Console .ReadLine()); "请输入第 {0} 个点的第 {1} 个像点 坐标： ", i + 1, j - 2); double .Parse( Console .ReadLine()); Console .WriteLine(); // 归算像点坐标 (i = 0; i < 4; i++) for (j = 3; j < 5; j++) if (j == 3) zuobiao[i, j] = zuobiao[i, j] - x0; else zuobiao[i, j] = zuobiao[i, j] - y0; // 计算和确定初值 double zs0 = m * f, xs0 = 0, ys0 = 0; for (i = 0; i < 4; i++) else using using using using x0 = y0 = int x0, y0, i, j; double f, m; Console .Write( " 请输入像片比例尺： "); double .Parse( Console .ReadLine()); Console .Write( " 请输入像片的内方位元素 x0:" ); // 均以毫米为单 位 int .Parse( Console .ReadLine()); Console .Write( " 请输入像片的内方位元素 y0:" ); int .Parse( Console .ReadLine()); Console .Write( " 请输入摄影机主距 f:" ); double .Parse( Console .ReadLine()); Console .WriteLine(); // 输入坐标数据 double [,] zuobiao = new double [4, 5]; (i = 0; i < 4; i++) for Console .Write( zuobiao[i, j] = Console .Write( zuobiao[i, j] = for 作业报告 空间后方交会 专业：测绘工程 班级：2008级(1)班姓名：陈闻亚 指导教师：陈强 2010 年 4 月16 日 1 作业任务------------------------------------------------------------------------------------ 3 2 作业思想 --------------------------------------------------------------------------------------- 3 3 作业条件及数据 -------------------------------------------------------------------- 3 4 作业过程--------------------------------------------------------------------------- 3 5 源程序----------------------------------------------------------------------------- 4 6 计算结果--------------------------------------------------------------------------- 17 7心得体会与建议----------------------------------------------------------------------------- 17 1 作业任务 计算近似垂直摄影情况下后方交会解。即利用摄影测量空间后方交会的方法，获取相片的6个外方位元素。限差为0.1。 2作业思想 利用摄影测量空间后方交会的方法求解。该方法的基本思想是利用至少三个一直地面控制点的坐标A（X A，Y A，Z A）、B（X B，Y B，Z B）C（X C，Y C，Z C），与其影像上对应的三个像点的影像坐标a（x a，y a）、b（x b，y b）、c（x c，y c），根据共线方程，反求该相片的外方位元素X S、Y S、Z S、φ、ω、κ。 3作业条件及数据 已知摄影机主距f=153.24mm，四对点的像点坐标与相应的地面坐标列入下表： 4作业过程 4．1 获取已知数据 相片比例尺1/m=1：10000，内方位元素f=153.24mm，x0，y0；获取控制点的地面测量坐标X t、Y t、Z t。 4．2 量测控制点的像点坐标： 本次作业中为已知。见表1。 一． 实验目的： 掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。 二． 仪器用具及已知数据文件： 计算机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）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素 的新值。 空间后方交会-空间前方交会程序编程实验一．实验目的要求 掌握运用空间后方交会-空间前方交会求解地面点的空间位置。学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程。 二．仪器用具 计算机、编程软件（MATLAB） 三．实验数据 实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。此四对坐标运用最小二乘法求解左右像片的外方位元素，即完成了空间后方的过程。另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。实验数据如下： 内方位元素：f=152.000mm，x0=0，y0=0 四．实验框图 此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0.00003，相当于0.1’的角度值）为止。在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。 在空间后方交会中运用的数学模型为共线方程 确定Xs，Ys，Zs的初始值时，对于左片可取地面左边两个GCP的坐标的平均值作为左片Xs 和Ys的初始值，取右边两个GCP的坐标平均值作为右片Xs 和Ys的初始值。Zs可取地面所有GCP的Z坐标的平均值再加上航高。 空间前方交会的数学模型为： 五．实验源代码 function Main_KJQHFJH() global R g1 g2 m G a c b1 b2; m=10000;a=5;c=4; feval(@shuru); %调用shuru（）shurujcp（）函数完成像点及feval(@shurujcp); %CCP有关数据的输入 XYZ=feval(@MQZqianfangjh); %调用MQZqianfangjh（）函数完成空间前方、%%%%%% 单位权中误差%%%% %后方交会计算解得外方位元素 global V1 V2; %由于以上三个函数定义在外部文件中故需VV=[]; %用feval（）完成调用过程 for i=1:2*c VV(i)=V1(i);VV(2*i+1)=V2(i); end m0=sqrt(VV*(VV')/(2*c-6)); disp('单位权中误差m0为正负:');disp(m0); %计算单位权中误差并将其输出显示 输入GCP像点坐标及地面摄影测量坐标系坐标的函数和输入所求点像点坐标函数： function shurujcp() global c m; m=input('摄影比例尺：'); %输入GCP像点坐标数据函数并分别将其c=input('GCP的总数='); % 存入到不同的矩阵之中 disp('GCP左片像框标坐标：'); global g1;g1=zeros(c,2); i=1; while i<=c m=input('x='); n=input('y='); g1(i,1)=m;g1(i,2)=n; i=i+1; end disp('GCP右片像框标坐标：'); global g2;g2=zeros(c,2); i=1; while i<=c m=input('x='); n=input('y='); g2(i,1)=m;g2(i,2)=n; i=i+1; end 空间后交-前交程序设计 (实验报告) 姓名： 班级： 学号： 时间： 空间后交-前交程序设计 一、实验目的 用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序 ⑴提交实习报告：程序框图、程序源代码、计算结果、体会 ⑵计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度 二、实验数据 f=150.000mm，x0=0，y0=0 三、实验思路 1.利用空间后方交会求左右像片的外方位元素 (1).获取m(于像片中选取两点，于地面摄影测量坐标系中选取同点，分别计算距离，距离比值即为m)，x,y,f,X,Y,Z (2).确定未知数初始值Xs,Ys,Zs,q,w,k (3).计算旋转矩阵R (4).逐点计算像点坐标的近似值（x），（y） (5).组成误差方程式 (6).组成法方程式 (7).解求外方位元素 (8).检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤（3）-（7） 2.利用求出的外方位元素进行空间前交，求出待定点地面坐标(1).用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1,R2 (2).根据左、右像片的外方位元素，计算摄影基线分量Bx，By，Bz (3).计算像点的像空间辅助坐标（X1,Y1,Z1）和（X2,Y2,Z2） (4).计算点投影系数N1和N2 (5).计算未知点的地面摄影测量坐标 四、实验过程 ⑴程序框图 函数AandL %求间接平差时需要的系数 %%%已知 %a=像点坐标x,b=像点坐标y，f内方位元素主距 %φ=q,ψ=w,κ=k %像空间坐标系X,Y,Z %地面摄影测量坐标系Xs,Ys,Zs function [A1,L1,A2,L2]=AandL(a,b,f,q,w,k,X,Y,Z,Xs,Ys,Zs) %%%%%%%%%%%选择矩阵元素 a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a3=-sin(q)*cos(w); b1=cos(w)*sin(k); b2=cos(w)*cos(k); b3=-sin(w); c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c3=cos(q)*cos(w); %%%%%%%共线方程的分子分母 X_=a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs); Y_=a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs); Z_=a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs); %%%%%%%近似值 x=-f*X_/Z_; y=-f*Y_/Z_; %%%%%%%A组成L组成 a11=1/Z_*(a1*f+a3*x); a12=1/Z_*(b1*f+b3*x); a13=1/Z_*(c1*f+c3*x); a21=1/Z_*(a2*f+a3*y); a22=1/Z_*(b2*f+b3*y); a23=1/Z_*(c2*f+c3*y); a14=y*sin(w)-(x/f*(x*cos(k)-y*sin(k))+f*cos(k))*cos(w); a15=-f*sin(k)-x/f*(x*sin(k)+y*cos(k)); a16=y; a24=-x*sin(w)-(y/f*(x*cos(k)-y*sin(k))-f*sin(k))*cos(w); a25=-f*cos(k)-y/f*(x*sin(k)+y*cos(k)); a26=-x; lx=a-x; ly=b-y; %%%%%%%%%组成一个矩阵，并返回 A1=[a11,a12,a13,a14,a15,a16]; A2=[a21,a22,a23,a24,a25,a26]; L1=lx; L2=ly; 函数deg2dms 空间后方交会的直接解 空间后方交会，即由物方已知若干个控制点以及相应的像点坐标，解求摄站的坐标与影像的方位，这是一个摄影测量的基本问题。通常采用最小二乘解算，由于原始的观测值方程是非线性的，因此，一般空间后方交会必须已知方位元素的初值，且解算过程是个迭代解算过程。但是，在实时摄影测量的某些情况下，影像相对于物方坐标系的方位是任意的，且没有任何初值可供参考。这时常规的空间后方交会最小二乘算法就无法处理，而必须建立新的空间后方交会的直接解法。 直接解法的基本思想是将它分成两步：先求出三个已知点i P 到摄站S 的距离i S ；然后求出摄站S 的坐标和影像方位。 物方一已知点()i i i i ,Z ,Y X P 在影像上的成像()i i i ,y x p ，根据影像已知的内方位元素()0 ,y f,x 可求得从摄站()S S S S ,Z ,Y X 到已知点i P 的观测方向i ,βαi 。 () ??? ????-+-= -=2 020 tan tan x x f y y βf x x αi i i i i (1) 距离方程组可以写成如下形式： ?? ??? =+++=+++=+++020202312 1133123232 3322322122 2211221b x x x a x b x x x a x b x x x a x (2) 其中()j ;i ,,i,j S ,b a ij ij ij ij ≠===321cos ?。因此，解算摄站S 到三个 控制点的距离问题，被归结为解算一个三元二次联立方程组的问题。这个方程组的解算方法选用迭代法。 迭代计算公式可写成： 中南大学 本科生课程设计(实践)任务书、设计报告 (摄影测量与遥感概论) 题目空间后方-前方交会 学生姓名 指导教师邹峥嵘 学院地球科学与信息物理学院 专业班级测绘0902班 学生学号 一、实验目的 通过对数字影像空间后交前交的程序设计实验，要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位元素及其精度，然后通过求得的外方位元素求解未知点的地面摄影测量坐标，达到通过摄影测量量测地面地理数据的目的。 二、实验要求 用C、VB或者Matlab编写空间后方交会-前方交会计算机程序。 提交实验报告：程序框图，程序源代码、计算结果及体会。 计算结果：地面点坐标、外方位元素及精度。 完成时间：2011年11月17日。 三、实验数据 f=150.000mm，x0=0，y0=0 四、实验思路 利用后方交会得出两张像片各自的外方位元素 1)获取已知数据：从摄影资料中插曲像片比例尺、平均航高、内 方位元素以及控制点的地面摄影测量坐标及对应的像点坐标。 2)确定未知数的初始值：在竖直摄影的情况下，胶原素的初始值 为0，线元素其中Zs=m*f+，Xs=，Ys=。 3)计算旋转矩阵R。 4)逐点计算像点坐标的近似值：利用共线方程。 5)组成误差方程并法化。 6)解求外方位元素。 7)检查计算是否收敛。 利用解求出的外方位元素进行前方交会 1)用各自像片的角元素计算出左右像片的旋转矩阵R1和R2。 2)根据左右像片的外方位元素计算摄影基线分量Bx,By，Bz。 3)逐点计算像点的空间辅助坐标。 单像空间后方交会程序报告 指导老师：刘老师 班级：测绘101 姓名：尚锋 学号： 19号 1、应用程序的主入口部分的代码： using System; using System.Collections.Generic; using System.Linq; using System.Windows.Forms; namespace单像空间后方交会 { static class Program { /// public double Z; } private double f; //主距 private double u; //u为外方位元素，下面5个相同 private double w; private double k; private double Xs; private double Ys; private double Zs; private image_point[] p = new image_point[4]; //四个控制点 private double[] R = new double[9]; //旋转矩阵 private double[] a = new double[8]; //像点坐标近似值 private double[,] A = new double[8, 6]; //误差方程式系数 private double[] L = new double[8]; //误差方程式常数项 private int count = 0; //统计代次数 public Tongyong(double g, double[] q) //构造函数，初始化各变量,单位m { f = g; for (int i = 0; i < 4; i++) { int j = i * 5; p[i].x = q[j]; p[i].y = q[j + 1]; p[i].X = q[j + 2]; p[i].Y = q[j + 3]; p[i].Z = q[j + 4]; } double ave = 0, sum = 0; //求比例尺分母 for (int i = 0; i < 3; i++) { for (int j = i + 1; j < 4; j++) { sum += Math.Sqrt(Math.Pow(p[i].X - p[j].X, 2) + Math.Pow(p[i].Y - p[j].Y, 2)) / Math.Sqrt(Math.Pow(p[i].x - p[j].x, 2) + Math.Pow(p[i].y - p[j].y, 2)); } } ave = sum / 6; u = 0; //给定外方位元素的初始值,角度均设置为0 w = 0; k = 0; Xs = (p[0].X + p[1].X + p[2].X + p[3].X) / 4; //Xs为四个控制点X的平均值，Ys类似空间后方交会的解算
摄影测量学后方交会matlab实习报告
摄影测量后方交会
单像空间后方交会和双像解析空间后方-前方交会的算法程序实现
空间后方交会编程实习报告
单像空间前方交会实习报告
摄影测量后方交会程序
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
摄影测量学空间后方交会实验报告测绘101徐斌
空间后方交会程序
摄影测量学单像空间后方交会程序设计作业
作业4--空间后方交会
空间后方交会程序
摄影测量实验报告(空间后方交会—前方交会)
空间后交-前交程序设计实验报告
空间后方交会的直接解
前方后方空间交会实验报告
单像空间后方交会程序报告
相关文档
- 摄影测量-空间前交、后交
- 前方交会实验报告(含VB程序代码)
- 单像空间前方交会实习报告
- 空间后交-前交程序设计实验报告
- 摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
- 前方交会实验报告
- 摄影测量前方交会后方交会matlab源程序
- 摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
- 空间前方交会及其应用
- 空间前方交会
- 空间前方交会及其应用
- 单像空间后方交会和双像解析空间后方-前方交会的算法程序实现
- 摄影测量实验报告(空间后方交会—前方交会)
- 摄影测量学 第三章 立体像对空间前方交会
- 工程摄影测量学--第13讲 空间前方交会
- 立体像对的空间前方交会
- 第3章 3节立体像对空间前方交会
- 前方后方空间交会实验报告
- 空间前方交会及其应用
- 中国大学MOOC慕课爱课程(18)--空间前方交会网课刷课