单向空间后方交会matlab编程
基于MATLAB的摄影交会测量计算程序设计与开发
产业科技创新 Industrial Technology Innovation78Vol.2 No.17产业科技创新 2020，2（17）：78~79Industrial Technology Innovation 基于MATLAB的摄影交会测量计算程序设计与开发袁 龙，曹雅迪（四川省测绘技术服务中心，四川 成都 610036）摘要：交会测量是测量学中一种常用的测量方法，在摄影测量中单张像片空间后方交会是摄影测量基本问题之一，其基本原理是通过地面已知的若干控制点和其在摄影像片上的相应像点坐标，通过计算的方法求解其摄站参数。
本文研究通过使用单像空间后方交会的基于共线条件方程的平差解法，使用安徽省某区域的实验数据，通过matlab软件实现了基于单像空间后方交会的程序设计，计算了其左右像片的收敛情况，结果表明其计算结果的数据精度和收敛情况较好，达到了预期的实验效果。
通过单像空间的后方交会的计算结果，文章进而根据其数据结果在matlab中构建了其空间前方交会的程序设计，求解实验区待定点的地面坐标，通过前方交会坐标变换情况，分析了其结果精度和需要注意的问题。
从总的研究结果来看，单张像片的空间后方交会是前方交会计算的基础，因此，单张像片空间后方交会成为解决这两个问题以及算法程序实现的关键。
关键词：摄影交会；后方交会；前方交会；程序设计中图分类号：U495 文献标识码：A 文章编号：2096-6164（2020）17-0078-02摄影测量是指通过搭载摄影设备，获取地物空间信息的手段完成测量作业，摄影测量的主要任务是通过获取二维的像片点空间坐标得到对应目标的实际三维坐标，完成数据的测量，后期通过数据处理和软件操作，应用到空间物体的三维建模、空间分析和综合制图等领域，有着广泛的社会应用和较高的使用价值。
单像空间后方交会的数据计算从目前来说主要有三种方法：基于共线条件方程的平差解法、基于直接线性变换的解法和角锥法。
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
程序运行环境为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.单模型光束法严密平差缺少已知数据进行验证，因此如果有已知数据请代入已知数据进行验证。
摄影测量——后方交会(matlab)说课讲解
摄影测量——后方交会(m a t l a b)3.用matlab编程结果function [XS,YS,ZS,p,w,k]=hfjhx_o=[-0.08615;-0.06899;-0.05340;0.08211;-0.01478;-0.07663;0.01064;0.06443];X=[36589.641,25273.32,2195.17;37631.08,31324.51,728.69;39100.97,24934.98,2386.50;40426.54,30319.81,757.31];f=0.15324;m=40000;XS=38437;YS=27963.155;ZS=f*m+1516.9175;p=0;w=0;k=0;R=zeros(3,3);num=1;dx=zeros(6,1);whilenum==1||abs(dx(4,1)*206225)>6||abs(dx(5,1)*206225)>6||abs(dx(6,1)*206 265)>6;R(1,1)=cos(p)*cos(k)-sin(p)*sin(w)*sin(k);R(1,2)=-cos(p)*cos(k)-sin(p)*sin(w)*sin(k);R(1,3)=-sin(p)*cos(w);R(2,1)=cos(w)*sin(k);R(2,2)=cos(w)*cos(k);R(2,3)=-sin(w);R(3,1)=sin(p)*cos(k)+cos(p)*sin(w)*sin(k);R(3,2)=-sin(p)*sin(k)+cos(p)*sin(w)*cos(k);R(3,3)=cos(p)*cos(w);x_s=zeros(8,1);x_s(1,1)=-f*(R(1,1)*(X(1,1)-XS)+R(2,1)*(X(1,2)-YS)+R(3,1)*(X(1,3)-ZS))/(R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS));x_s(2,1)=-f*(R(1,2)*(X(1,1)-XS)+R(2,2)*(X(1,2)-YS)+R(3,2)*(X(1,3)-ZS))/(R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS));x_s(3,1)=-f*(R(1,1)*(X(2,1)-XS)+R(2,1)*(X(2,2)-YS)+R(3,1)*(X(2,3)-ZS))/(R(1,3)*(X(2,1)-XS)+R(2,3)*(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS));x_s(4,1)=-f*(R(1,2)*(X(2,1)-XS)+R(2,2)*(X(2,2)-YS)+R(3,2)*(X(1,3)-ZS))/(R(1,3)*(X(2,1)-XS)+R(2,3)*(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS));x_s(5,1)=-f*(R(1,1)*(X(3,1)-XS)+R(2,1)*(X(3,2)-YS)+R(3,1)*(X(1,3)-ZS))/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)-ZS));x_s(6,1)=-f*(R(1,2)*(X(3,1)-XS)+R(2,2)*(X(3,2)-YS)+R(3,2)*(X(1,3)-ZS))/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)-ZS));x_s(7,1)=-f*(R(1,1)*(X(4,1)-XS)+R(2,1)*(X(4,2)-YS)+R(3,1)*(X(1,3)-ZS))/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS));x_s(8,1)=-f*(R(1,2)*(X(4,1)-XS)+R(2,2)*(X(4,2)-YS)+R(3,2)*(X(1,3)-ZS))/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS));a=zeros(8,6);for i=1:4Z=(R(1,3)*(X(i,1)-XS)+R(2,3)*(X(i,2)-YS)+R(3,3)*(X(i,3)-ZS));a(i*2-1,1)=1/Z*(R(1,1)*f+R(1,3)*x_s(2*i-1,1));a(i*2-1,2)=1/Z*(R(2,1)*f+R(2,3)*x_s(2*i-1,1));a(i*2-1,3)=1/Z*(R(3,1)*f+R(3,3)*x_s(2*i-1,1));a(2*i,1)=1/Z*(R(1,2)*f+R(1,3)*x_s(2*i,1));a(2*i,2)=1/Z*(R(2,2)*f+R(2,3)*x_s(2*i,1));a(2*i,3)=1/Z*(R(3,2)*f+R(3,3)*x_s(2*i,1));a(i*2-1,4)=x_s(2*i,1).*sin(w)-(x_s(2*i-1,1)./f.*(x_s(2*i-1,1).*cos(k)-x_s(2*i,1).*sin(k))+f*cos(k))*cos(w);a(i*2-1,5)=-f*sin(k)-x_s(2*i-1,1)./f.*(x_s(2*i-1,1).*sin(k)+x_s(2*i,1).*cos(k));a(i*2-1,6)=x_s(2*i,1);a(2*i,4)=-x_s(2*i-1,1).*sin(w)-(x_s(2*i,1)./f.*(x_s(2*i-1,1).*cos(k)-x_s(2*i,1).*sin(k))-f*sin(k))*cos(w);a(2*i,5)=-f*sin(k)-x_s(2*i,1)./f.*(x_s(2*i-1,1).*sin(k)+x_s(2*i,1).*cos(k));a(2*i,6)=-x_s(i*2-1,1);endL=x_o-x_s;dx=((a'*a)\a')*L;num2str(a);disp('改正数');num2str(dx);size(dx);XS=XS+dx(1,1);YS=YS+dx(2,1);ZS=ZS+dx(3,1);p=p+dx(4,1);w=w+dx(5,1);k=k+dx(6,1);num=num+1;num2str(XS)num2str(YS)num2str(ZS)num2str(p)num2str(w)num2str(k)End程序中p=ψ(弧度),w=ω(弧度) k=κ(弧度)解析法相对定向是通过计算相对定向元素建立地面立体模型。
摄影测量-空间前交、后交【精选文档】
空间后交—前交程序设计(实验报告)姓名：班级：学号：时间:空间后交-前交程序设计一、实验目的用 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,Zsfunction ［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％%％%％%％％角度转度分秒function y=deg2dms(x）a=floor(x)；b=floor(（x-a)＊60);c=(x-a—b/60）＊3600；y=a+（b/100）+(c/10000);函数dms2deg%%％％%度分秒转度function y=dms2deg（x）a=floor(x）；b=floor((x-a)＊100）；c=（x-a—b/100)*10000；y=a+b/60+c/3600;函数ok%％％％％%%％%%%％％％目的是为了保证各取的值的有效值%％xy为n*1，a为1＊nfunction result=ok(xy，a)format short gi=size（xy，1)；for n=1:io=xy（n)—floor(xy(n,1));o=round（o＊（10^a（n))）/(10^a(n）);xy(n，1)=floor(xy（n,1）)+o;endformat long gresult=xy;函数rad2dmsxy%%％%求度分秒表现形式的三个外方位元素，三个角度function xydms=rad2dmsxy(xy)［a,b,c，d，e,f]=testvar(xy）;d=deg2dms（rad2deg(d））；e=deg2dms（rad2deg(e))；f=deg2dms(rad2deg(f）);xydms=[a,b，c,d,e，f］';函数spacehoujiao％％%％%%%空间后交％%% f%％输入p（2*n,1）%％像点坐标x,y,X,Y，Z,均为（n,1）function ［xy,m，R］=spacehoujiao(p，x,y，f,X，Y,Z）format long；%%％%％权的矢量化,这是等精度时的，如果非,将函数参数改为PP=diag（p)；％％求nj=size（X，2）；%％初始化Xs=0；Ys=0;Zs=0;for n=1：jXs=Xs+X（n）；Ys=Ys+Y(n）;Zs=Zs+Z(n）;endSx=sqrt((x（2）-x(1）)^2+(y（2)—y(1)）^2);％%％%两像点之间距离Sd=sqrt(（X(2）-X（1))^2+(Y（2）-Y（1))^2);%％%％两地面控制点之间距离m=Sd/Sx; ％％％％图像比例系数Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j；m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6）；l=rand（2＊j，1);％％％％for n=1:j[a（2＊n—1,:)，l(2*n—1,1）,a(2＊n,：），l(2*n,1）］=AandL(x(n)，y（n),f,q，w，k，X(n)，Y（n），Z(n），Xs,Ys,Zs);enddet=inv（a’＊P*a）*transpose（a)*P*l;％%%%%%％％％循环体while 1％％％％％%%％%％％％％%％％[dXs,dYs,dZs，dq,dw,dk]=testvar(det）;detXs=abs（dXs）；detYs=abs（dYs)；detZs=abs(dZs)；detq=abs(dq)；detw=abs（dw)；detk=abs（dk)；%%％％%%％%%if （（detXs<0。
编写单片空间后方交会程序实习一
实验一 编写单片空间后方交会程序
一、实验目的
用程序设计语言（Visual C++或者C 语言、matlab 语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。
本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
二、实验内容与要求
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。
三、实验原理
根据摄影过程的几何反转原理。
四、实验数据准备
﹡ 已知航摄仪的内方位元素0,24.15300===y x mm f （摄影比例尺为1:50000）； ﹡ 4个地面控制点的地面坐标及其对应像点的像片坐标：
1．每人必须独立上机操作；
2．实验前认真复习单片空间后方交会的有关内容。
六、实验报告
原理、实验步骤及程序框图
上机调试程序并打印结果。
后方交会MATLAB程序实习报告
《摄 影 测 量 学》单像空间后方交会实习报告班 级： XXXX姓 名： X X X学 号：XXXXXXXXXXXXX指导教师： X X X一、实习目的1、掌握空间后方交会的定义和实现算法；2、了解摄影测量平差的基本过程；3、熟练MATLAB 等程序编写。
二、实习原理利用至少三个已知地面控制点的坐标),,(A A A Z Y X A 、),,(B B B Z Y X B 、),,(C C C Z Y X C ，与其影像上对应的三个像点的影像坐标),(a a y x a 、),(b b y x b 、),(c c y x c ，根据共线方程，反求该像片的外方位元素κωϕ、、、、、S S S Z Y X 。
共线条件方程式：将共线式线性化并取一次小值项得：三、解算过程①获取已知数据。
包括影像比例尺1/m，平均摄影距离（航空摄影的航高）H，内方位元素x0、y0、f，控制点的空间坐标X、Y、Z。
②量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
③确定未知数的初始值。
单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。
或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。
④计算旋转矩阵R。
利用角元素近似值计算方向余弦值，组成R阵。
⑤逐点计算像点坐标的近似值。
利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。
⑥逐点计算误差方程式的系数和常数项，组成误差方程式。
⑦计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。
⑧解求外方位元素。
根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。
⑨检查计算是否收敛。
将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.000001弧度，当3个改正数均小于0.000001弧度时，迭代结束。
摄影测量学后方交会matlab实习报告
摄影测量原理单张影像后方交会实习目录一实习目的 (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 个方向余弦。
单像空间后方交会实验报告(c++版)
单像空间后方交会实验报告(c++版)单像空间后方交会姓名：学号：时间：目录一、作业任务..................................................... - 4 -二、计算原理..................................................... - 4 -三、算法流程..................................................... - 8 -四、源程序....................................................... - 9 -五、计算结果..................................................... - 9 -六、结果分析..................................................... - 9 -七、心得与体会................................................... - 9 -八、附页......................................................... - 9 -1.c++程序.................................................. - 10 -2.C++程序截图.............................................. - 17 -3.matlb程序 ............................................... - 17 -一、 作业任务已知条件：摄影机主距f=153.24mm ，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。
摄影测量学单像空间后方交会编程实习报告
摄影测量学单像空间后方交会编程实习报告实习背景在本次实习中，我们学习了摄影测量学单像空间后方交会的编程实现。
这是一种通过计算影像中各点的坐标来确定被摄物的三维坐标的方法，应用广泛于测绘、地理信息、建筑等领域。
本次实习采用 MATLAB 软件进行编程，目的是将理论知识应用到实际操作中，让我们更深入理解摄影测量学单像空间后方交会的原理和应用。
实习内容理论部分首先，我们在工作室进行了理论部分的学习。
老师讲解了单像空间后方交会的原理，以及如何通过影像坐标、相机外方位元素、像点坐标和像平面坐标等参数来计算被摄物的三维坐标。
在理论部分的学习过程中，我们通过公式的推导和实例分析，更加深入地理解了单像空间后方交会的原理。
实践部分实践部分是本次实习的重头戏。
我们利用 MATLAB 软件进行了单像空间后方交会的编程实现，具体步骤如下：1.输入相机外方位元素通过读取文本文件，将相机外方位元素（相机在拍摄时的姿态、位置等参数）输入到 MATLAB 中。
2.输入影像坐标通过读取文本文件，将影像中的像点坐标输入到 MATLAB 中。
3.计算像平面坐标利用相机内定标参数，将像点坐标转化为像平面坐标。
4.计算被摄物三维坐标根据单像空间后方交会的原理，利用相机外方位元素、像平面坐标和像点坐标等参数，计算被摄物的三维坐标。
5.输出结果将计算结果输出到文本文件中，以便后续的数据处理和分析。
在实际操作中，我们首先编写了 MATLAB 脚本文件，根据上述步骤逐步实现了单像空间后方交会的计算过程。
然后，我们利用自己拍摄的实际照片进行实验，将相机外方位元素和像点坐标输入到程序中，最终得到了被摄物的三维坐标结果。
实习收获通过本次实习，我从理论到实践，更深入地理解了摄影测量学单像空间后方交会的原理和应用，同时也掌握了 MATLAB 的编程技能。
在实践中，我遇到了许多问题，包括数据的输入输出、代码的调试和结果的分析等等。
通过和同学的讨论和老师的指导，我不仅解决了这些问题，还对摄影测量学的应用有了更深入的认识。
空间后交-前交程序设计实验报告
空间后交-前交程序设计(实验报告)姓名：班级：学号：时间：空间后交-前交程序设计一、实验目的用 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,Zsfunction [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%%%%%%%%角度转度分秒function y=deg2dms(x)a=floor(x);b=floor((x-a)*60);c=(x-a-b/60)*3600;y=a+(b/100)+(c/10000);函数dms2deg%%%%%度分秒转度function y=dms2deg(x)a=floor(x);b=floor((x-a)*100);c=(x-a-b/100)*10000;y=a+b/60+c/3600;函数ok%%%%%%%%%%%%%%目的是为了保证各取的值的有效值%%xy为n*1,a为1*nfunction result=ok(xy,a)format short gi=size(xy,1);for n=1:io=xy(n)-floor(xy(n,1));o=round(o*(10^a(n)))/(10^a(n));xy(n,1)=floor(xy(n,1))+o;endformat long gresult=xy;函数rad2dmsxy%%%%求度分秒表现形式的三个外方位元素，三个角度function xydms=rad2dmsxy(xy)[a,b,c,d,e,f]=testvar(xy);d=deg2dms(rad2deg(d));e=deg2dms(rad2deg(e));f=deg2dms(rad2deg(f));xydms=[a,b,c,d,e,f]';函数spacehoujiao%%%%%%%空间后交%%% f%%输入p（2*n,1）%%像点坐标x,y,X,Y,Z,均为（n,1）function [xy,m,R]=spacehoujiao(p,x,y,f,X,Y,Z)format long;%%%%%权的矢量化,这是等精度时的，如果非，将函数参数改为P P=diag(p);%%求nj=size(X,2);%%初始化Xs=0;Ys=0;Zs=0;for n=1:jXs=Xs+X(n);Ys=Ys+Y(n);Zs=Zs+Z(n);endSx=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);%%%%两像点之间距离Sd=sqrt((X(2)-X(1))^2+(Y(2)-Y(1))^2);%%%%两地面控制点之间距离m=Sd/Sx; %%%%图像比例系数Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j;m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6);l=rand(2*j,1);%%%%for n=1:j[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n ),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;%%%%%%%%%循环体while 1%%%%%%%%%%%%%%%%[dXs,dYs,dZs,dq,dw,dk]=testvar(det);detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);%%%%%%%%%if((detXs<0.01)&&(detYs<0.01)&&(detZs<0.01)&&(detq<pi/648000)&&(detw<pi/648000)&& (detk<pi/648000))break;elseV=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt((V'*P*V)/(2*j-6));%%m0需要每次的改正数算出来相加%%%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%%%for n=1:j[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n ),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;i=i+1;%%%%end%%%end[dXs,dYs,dZs,dq,dw,dk]=testvar(det);detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);V=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt((V'*P*V)/(2*n-6));%%%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%%%%%%%%%%%%%可以输出迭代次数的i%%%%%%%%%%%%Xs,Ys,Zs,q,w,k,i,dXs,dYs,dZs,dq,dw,dk,detXs,detYs,detZs %%%%%%%%%%%精度mo=m0*sqrt(Q);m=[mo(1,1),mo(2,2),mo(3,3),mo(4,4),mo(5,5),mo(6,6)]';[mXs,mYs,mZs,mq,mw,mk]=testvar(m);%%%%%%%%%输出xy=[Xs,Ys,Zs,q,w,k]';%%输出（6,1）的外方位元素m=[m0,mXs,mYs,mZs,mq,mw,mk]';%%单位误差，各元素中误差R=xyR(xy);%%旋转矩阵函数spaceqianjiao%%空间前交%输入f%输入x1,y1,x2,y2,R1,R2,xy1,xy2 （n,1）%输出X,Y,Z （n,1）function [X,Y,Z]=spaceqianjiao(x1,y1,x2,y2,f,R1,R2,xy1,xy2) i=size(x1,2);[Xs1,Ys1,Zs1,q1,w1,k1]=testvar(xy1);[Xs2,Ys2,Zs2,q2,w2,k2]=testvar(xy2);for n=1:i[X1(n),Y1(n),Z1(n)]=testvar(R1*[x1(n),y1(n),-f]');[X2(n),Y2(n),Z2(n)]=testvar(R2*[x2(n),y2(n),-f]');Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;N1=(Bx*Z2(n)-Bz*X2(n))/(X1(n)*Z2(n)-X2(n)*Z1(n));N2=(Bx*Z1(n)-Bz*X1(n))/(X1(n)*Z2(n)-X2(n)*Z1(n));X(n)=Xs1+N1*X1(n);Z(n)=Zs1+N1*Z1(n);Y(n)=0.5*((Ys1+N1*Y1(n))+(Ys2+N2*Y2(n)));end函数testvar%分割矩阵。
摄影测量实验报告(空间后方交会—前方交会)
空间后方交会—空间前方交会程序编程实验一．实验目的要求掌握运用空间后方交会-空间前方交会求解地面点的空间位置.学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。
然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程.二．仪器用具计算机、编程软件(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＊cVV(i)=V1（i)；VV(2*i+1）=V2(i)；endm0=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<=cm=input('x='）；n=input（'y='）;g1(i，1)=m;g1（i，2)=n；i=i+1；enddisp('GCP右片像框标坐标：’);global g2；g2=zeros(c,2);i=1;while i〈=cm=input（'x=’）;n=input('y=’）;g2（i，1)=m；g2（i，2）=n;i=i+1;end％%%%％%%%%%%%%％％%%％％%％％％％%％%%％％%%％％%％％%%％% function shuru（)global a；a=input（'计算总像对点数='); %完成想计算所需的像平面坐标global b1；％坐标输入,存入不同的矩阵中b1=zeros(a,2)；disp('左片像点坐标：'）i=1；while i〈=am=input（'x=’）；n=input(’y=’);b1(i，1)=m；b1（i,2)=n;i=i+1;end%%global b2;b2=zeros（a,2）；disp（’右片像点坐标：')i=1;while i〈=am=input('x=’）;n=input（'y=’）；b2(i,1）=m；b2(i,2)=n；i=i+1;end％%global c;c=input(’GCP的总数='）；disp('GCP摄影测量系坐标：’)global G;G=zeros(3,c)；i=1；while i〈=cm=input（’X='）;n=input(’Y='）;v=input（’Z=')；G(i，1）=m;G（i，2)=n;G（i,3)=v；i=i+1;end %%％％％％%%％%%%%%%％％％％％％％％%%％％%％%%%％％%%%空间前方交会和后方交会函数:function XYZ=MQZqianfangjh（）global R1 R2 a f b1 b2 Ra Rb；global X1 X2;R1=Ra；R2=Rb;R1=zeros（3,3）;R2=zeros（3,3)；global g1 g2 G V1 V2 V WF c QXX QXX1 QXX2；xs0=（G（1，1)+G(3,1))/2；ys0=（G（1，2）+G(3，2)）/2;［Xs1，Ys1，Zs1,q1,w1，k1 R]=houfangjh(g1，xs0,ys0）；％对左片调用后方交会函数R1=R；V1=V；WF1=WF;QXX1=QXX;save '左片外方位元素为。
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差).doc
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差) .程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug 下。
1. 单像空间后方交会原始数据：内方位元素x0(/mm)y0(/mm)主距f(/mm)比例尺分母00153.2450000像点坐标(/mm)物点坐标(/m)-内方位元素x0(/mm)y0(/mm)主距f(/mm)比例尺分母00153.2450000像点坐标(/mm)物点坐标(/m)：#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=6.0/206265.0||xcorrect[5][0]=6.0/206265.0||xcorrect[9][0]=6.0/ 206265.0||xcorrect[10][0]=6.0/206265.0||xcorrect[11][0]=6.0/206265.0); savedata(row,controlpoint,A[0],ATA[0],L[0],diedainumber,Xs1,Ys1,Zs1,f ai1,oumiga1,kapa1,Xs2,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");}测试结果：如果有正确的数据，请代入进行验证。
2020年MATLAB编程与应用实验报告(交会定点)
1交会定点实验报告所属课程名称 MATLAB编程与应用实验地点实验日期 21119班级学号姓名指导老师一、实验目的交会定点包括前方交会和后方交会的计算，通过编写相应函数，实现对标量和向量（或矩阵）输入参数的前方交会和后方交会计算。
二、实验内容【实验过程及成果】（程序说明、实验代码、实验数据、实验结果）程序说明前方交会x1、y1已知点A的坐标,a已知点A的交会角,x2、y2已知点B的坐标,b已知点B的交会角，dms_rad角度转弧度，通过公式计算计算P点坐标。
后方交会xA、yA已知点A的坐标,xC、yC已知点C的坐标,xB、yB已知点B的坐标,A、B 为交会角，dms_rad角度转弧度，通过公式计算计算P点坐标。
实验代码>>function [X,Y]=QFJH(x1,y1,a,x2,y2,b)%前方交会a=dms_rad(a);b=dms_rad(b);X=x*cot(b)+x*cot(a)+(y2-y1);X=Y=y*cot(b)+y*cot(a)-(x2-x1);Y=end>>function [X, Y]=HFJH(xA,yA,xC,yC,xB,yB,A,B)%后方交会A=dms_rad(A);B=dms_rad(B);a=(yA-yC).*cot(A)+(xA-xC);b=(xA-xC).*cot(A)-(yA-yC);c=(xB-xC).*cot(B)-(yB-yC);d=(yB-yC).*cot(B)-(xB-xC);k=(a+d)dx=(a-b.*k)m=a-b.*k;n=c.*k-d;X=xC+dx;Y=yC+k.*dx;实验数据>>[X1,Y1]=QFJH(367543,25614,423145,369675,273126,543121); >>[X2,Y2]=HFJH(1,1,3,2,2,3,3,4);实验结果【实验小结】（收获体会）通过此次实验了解了交会定点的基本计算公式和相应的程序编写，对交会定点的不同计算公式适应不同的程序编写有了了解，学习了用代码来实现前方交会和后方交会的计算。
摄影测量空间后交程序
实验内容:用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序程序代码:#include "stdio.h"#include "math.h"#include "iostream.h"#define N 4#define M 2*Nvoid 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 inv(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 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;}void main(){int i,m,count;double mm[M+1];double t,w,k,Xs0,Ys0,Zs0,f,S1=0.0,S2=0.0;doubleH[6]={1},a[3],b[3],c[3],Xo[N],Yo[N],Zo[N],A[6*M],B[6*M],l[M],C[36], D[6];double x[N],y[N],X[N],Y[N],Z[N];cout<<"请输入焦距(mm):f=";cin>>f;cout<<"请输入影像比例尺分母:m=";cin>>m;x[0]=-73.93; y[0]=78.706; X[0]=5083.205; Y[0]=5852.099; Z[0]=527.925; x[1]=-5.252; y[1]=78.184; X[1]=5780.02; Y[1]=5906.365; Z[1]=571.549; x[2]=-79.122; y[2]=-78.879; X[2]=5210.879; Y[2]=4258.446; Z[2]=461.81; x[3]=-9.887; y[3]=-80.089;X[3]=5909.264;Y[3]=4314.283;Z[3]=455.484;for(i=0;i<N;i++) {x[i]=x[i]/1000.0; y[i]=y[i]/1000.0; S1 += X[i];S2 += Y[i];}t=w=k=0.0;Xs0=S1/N;Ys0=S2/N;f=f/1000.0;Zs0=m*f;count=0;while(fabs(H[0])>0.00001||fabs(H[1])>0.00001||fabs(H[2])>0.00001||f abs(H[3 ])>0.00001||fabs(H[4])>0.00001||fabs(H[5])>0.00001){count+=1;a[0]=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);a[1]=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k);a[2]=-sin(t)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(t)*cos(k)+cos(t)*sin(w)*sin(k);c[1]=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k);c[2]=cos(t)*cos(w);for(i=0;i<N;i++){Xo[i]=-f*(a[0]*(X[i]-Xs0)+b[0]*(Y[i]-Ys0)+c[0]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[ 2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0));Yo[i]=-f*(a[1]*(X[i]-Xs0)+b[1]*(Y[i]-Ys0)+c[1]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y [i]-Ys0)+c[2]*(Z[i]-Zs0));Zo[i]=a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0);A[12*i+0]=(a[0]*f+a[2]*x[i])/Zo[i];A[12*i+1]=(b[0]*f+b[2]*x[i])/Zo[i];A[12*i+2]=(c[0]*f+c[2]*x[i])/Zo[i];A[12*i+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[12*i+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+5]=y[i];A[12*i+6]=(a[1]*f+a[2]*y[i])/Zo[i];A[12*i+7]=(b[1]*f+b[2]*y[i])/Zo[i];A[12*i+8]=(c[1]*f+c[2]*y[i])/Zo[i];A[12*i+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[12*i+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+11]=-x[i];l[2*i]=x[i]-Xo[i];l[2*i+1]=y[i]-Yo[i]; }transpose(A,B,M,6); mult(B,A,C,6,M,6); mult(B,l,D,6,8,1); inv(C,6);mult(C,D,H,6,6,1); Xs0+=H[0];Ys0+=H[1];Zs0+=H[2];t+=H[3];w+=H[4];k+=H[5];}for(i=0;i<M;i++) {mm[i]=0;}for(i=0;i<M;i++){mm[0]+=l[i]*l[i];}mm[0]=sqrt(mm[0]/(2*N-6)); for(i=0;i<6;i++) {mm[i+1]=sqrt(C[6*i+i])*mm[0];}cout<<"总的叠代次数为:"<<count<<endl;cout<<"像主点的空间坐标为:"<<endl;cout<<"Xs="<<Xs0<<endl;cout<<"Ys="<<Ys0<<endl;cout<<"Zs="<<Zs0<<endl;cout<<"t="<<t<<endl;cout<<"w="<<w<<endl;cout<<"k="<<k<<endl;cout<<"单位权中误差为:"<<mm[0]<<endl;for(i=1;i<7;i++){cout<<"第"<<i<<"个未知数的中误差为:"<<mm[i]<<endl; }}。
实习一 单张影像空间后方交会程序设计
实习一单张影像空间后方交会程序设计一、实习目的
学生通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，为该门课程的学习打下好的基础。
二、实习内容
根据实习要求及试验数据实现单张影像空间后方交会的计算，并将最后计算得到的旋转矩阵及外方位元素输出。
图１空间后方交会程序流程图
三、预备知识
1、熟悉一门编程语言；
2、分析已有数据的特点；
3、弄清矩阵的存储与操作方法；
4、矩阵的运算模块（函数）可直接在相关资料上查找并直接使用。
四、实习步骤
已有的控制点像点坐标及物方坐标可通过输入界面进行输入，已可定义成文本文件由程序读取文本文件；定义变量和数组来存储数据；最终计算结果的输出（可输出为文本文件）。
实习结束后提交每人提交一份实习报告。
五、思考题
1、矩阵在计算机中是怎样存储的？怎样实现对矩阵中每个元素的操作（读取与存储）？
2、单张影像空间后方交会计算时在什么条件下迭代结束？
六、实验数据
已知航摄仪的内方位元素：f=153.24mm，x0=y0=0.0mm，摄影比例尺为1:50000；4个地面控制点的地面坐标及其对应像点的像片坐标如下表。
试验一 单片空间后方交会程序设计
三、存储数据的变量及数组
针对空间后方交会计算中数据的特点，需定义相 关的变量和数组来存储数据。 存储已知数据的数组：x[]，y[]分别用来存储控制点 所对应的像点坐标；X[]，Y[]，Z[]分别用来存储地 面控制点坐标； 存储旋转矩阵数组：a[]，b[], c []分别用来存储旋 转矩阵的九个方向余弦值 误差方程系数数组： A[]数组用来存储误差方程系数 误差方程常数项数组：l []数组用来存储误差方程常 数项
解法方程，求未知数改正数 计算改正后的外方位元素 否 未知数改正数<限差否？ 整理并输出计算结果 正常结束
输出中间结果 和出错信息 非正常结束
已知航摄仪的内方位元素：f=153.24mm， x0=y0=0.0mm，摄影比例尺为1:50000； 4个地面控制 点的地面坐标及其对应像点的像片坐标如下表：
存储外方位元素改正值的数组：H[]存储某次迭 代计算出的参数改正值 相关变量：Xs0, Ys0, Zs0分别用来存储三个外 方位线元素；t，w，k分别用来存储三个外方 位角元素；f用来存储主距；m存储摄影比例尺 分母
用程序设计语言编写一个完整的单片空间后方交会程序，通 过对提供的试验数据进行计算，输出像片的外方位元素。
二、迭代计算的思想
迭代计算一般需要事先知道参数的初值，初 值选的不好，可能迭代计算收敛速度较慢，甚 至不收敛。
每计算出参数的改正值后，需要把改正值与 限差（迭代收敛的条件）进行比较，直到改正 值小于限差，迭代计算结束。
试验一 单片空间后方交会程序设计
重点： 1.理解单片空间后方交会计算流程 2.掌握迭代计算的思想 3.数据读入与存储方法
输入原始数据 归算像点坐标 x, y
一、空间后方交会 计算程序框图
是 迭 代 次 Ys0 , Z s0 , 0 , 0 , 0 组成旋转矩阵R 计算 ( x), ( y )和 l x , l y 直接生成误差方程B和L 否 所有点完否？
单向空间后方交会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);。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
%读取包含点号的像点坐标、控制点坐标分别到矩阵IP、CP% IP=load('f:');
CP=load('f:');
%删除像点坐标、控制点坐标矩阵中的点号列%
IP(:,1)=[];
CP(:,1)=[];
%像片大小%
H=4008;
W=5344;
IP1=IP(:,1)-W/2;
￥
IP2=H/2-IP(:,2);
IP=[IP1,IP2];
%焦距%
fx=;
fy=;
%内方位元素%
x0=;
y0=;
%畸变参数%
k1=;
{
k2=;
p1=;
p2=;
%外方位元素初值%
Xs=3000;
Ys=-100;
Zs=100;
Phi=0;
Omega=0;
Kappa=0;
!
m=size(IP,1);
AIP=zeros(m,2);
L=zeros(2*m,1);
A=zeros(2*m,6);
X=ones(6,1);
while X(4,1)>=||X(5,1)>=||X(6,1)>=
%旋转矩阵的系数%
a1=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);
》
b1=cos(Omega)*sin(Kappa);
b2=cos(Omega)*cos(Kappa);
b3=-sin(Omega);
c1=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa);
c2=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa);
c3=cos(Phi)*cos(Omega);
for i=1:m
%求像点坐标常数项%
r=sqrt((IP(i,1)-x0)^2+(IP(i,2)-y0)^2);
Dx=(IP(i,1)-x0)*(k1*(r^2)+k2*(r^4))+p1*(r^2+2*((IP(i,1)-x0)^2))+2*p2*(IP(i,1)-x0)*(IP(i,2)-y0);
、
Dy=(IP(i,2)-y0)*(k1*(r^2)+k2*(r^4))+p2*(r^2+2*((IP(i,2)-y0)^2))+2*p1*(IP(i,1)-x0)*(IP(i,2)-y0);
X_=a1*(CP(i,1)-Xs)+b1*(CP(i,2)-Ys)+c1*(CP(i,3)-Zs);
Y_=a2*(CP(i,1)-Xs)+b2*(CP(i,2)-Ys)+c2*(CP(i,3)-Zs);
Z_=a3*(CP(i,1)-Xs)+b3*(CP(i,2)-Ys)+c3*(CP(i,3)-Zs);
%求像点坐标的近似值%
AIP(i,1)=-fx*X_/Z_+Dx+x0;
AIP(i,2)=-fy*Y_/Z_+Dy+y0;
%求误差方程式系数%
a11=(a1*fx+a3*IP(i,1))/Z_;
a12=(b1*fx+b3*IP(i,1))/Z_;
【
a13=(c1*fx+c3*IP(i,1))/Z_;
a14=sin(Omega)*IP(i,2)-(IP(i,1)*(IP(i,1)*cos(Kappa)-IP(i,2)*sin(Kappa))/fx+fx*cos(Kappa))*cos(Om ega);
a15=-fx*sin(Kappa)-IP(i,1)*(IP(i,1)*sin(Kappa)+IP(i,2)*cos(Kappa))/fx;
a16=IP(i,2);
a21=(a2*fy+a3*IP(i,2))/Z_;
a22=(b2*fy+b3*IP(i,2))/Z_;
a23=(c2*fy+c3*IP(i,2))/Z_;
a24=-sin(Omega)*IP(i,1)-(IP(i,2)*(IP(i,1)*cos(Kappa)-IP(i,2)*sin(Kappa))/fy-fy*sin(Kappa))*cos(Om ega);
a25=-fy*cos(Kappa)-IP(i,2)*(IP(i,1)*sin(Kappa)+IP(i,2)*cos(Kappa))/fy;
a26=-IP(i,1);。
%求常数项%
L(2*i-1,1)=IP(i,1)-AIP(i,1);
L(2*i,1)=IP(i,2)-AIP(i,2);
%求误差方程式系数阵%
A(2*i-1,1)=a11;
A(2*i-1,2)=a12;
A(2*i-1,3)=a13;
A(2*i-1,4)=a14;
A(2*i-1,5)=a15;
A(2*i-1,6)=a16;
^
A(2*i,1)=a21;
A(2*i,2)=a22;
A(2*i,3)=a23;
A(2*i,4)=a24;
A(2*i,5)=a25;
A(2*i,6)=a26;
end
%求改正数%
X=pinv(A'*A)*A'*L;
%求真值%
@
Xs=Xs+X(1,1);
Ys=Ys+X(2,1);
Zs=Zs+X(3,1);
Phi=Phi+X(4,1);
Omega=Omega +X(5,1);
Kappa=Kappa +X(6,1);
end
%求误差矩阵%
V=A*X-L;
%精度评定%
>
m0=abs(sqrt(V'*V/(2*m-6)));
Q=inv(A'*A);
mi=m0*sqrt(Q);
%{fprintf('%.6f\n',A);%
%fprintf('\n');%
fprintf('%.6f\n%.6f\n%.6f\n%.6f\n%.6f\n%.6f\n',Xs,Ys,Zs,Phi,Omega,Kappa); fprintf('\n');
%fprintf('%.6f\n',L);%
%fprintf('\n');%
%fprintf('%.6f\n',V);%
%fprintf('\n');%
fprintf('%.6f\n',mi);。