近景单张像片空间后方交会

实验一近景单张像片空间后方交会一、实验目的通过单张像片空间后方交会算法编程，掌握空间后方交会的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。

二、实验内容利用C++编程平台，通过给定的单片像点的像点坐标、内方位元素和控制点物方空间坐标，计算出像片的外方位元素。

三、实验步骤：1、编程流程图：2、编程代码：#include "iostream"#include"stdio.h"#include "stdlib.h"#include<math.h>#define N 10using namespace std;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;}/* ######################矩阵求逆##################*/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);}/*------------------------------矩阵转置----------------------------------*/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 main(){double Xs,Ys,Zs,q,w,k;//外方位元素//double a[3],b[3],c[3];//旋转矩阵//double x0,y0,f;//内方位元素//double x[N],y[N];double X[N],Y[N],Z[N];double x1[N],y1[N];double m;//中误差//double L[2*N];double XX[6];double A[2*N][6];double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];//外方位元素初始置//int i,n=0;//迭代次数//double sum=0,m0;/*for(i=0;i<N;i++){printf("请输入第%d个点的地面坐标：",i+1); //----------输入点地面坐标-scanf_s("%lf%lf%lf",&X[i],&Y[i],&Z[i]);} */X[0]=1617.712 ;Y[0]=1883.870 ;Z[0]=-4917.308;X[1]=1615.572 ;Y[1]= 1687.514 ; Z[1]= -4913.750;X[2]=1614.640 ;Y[2]= 1315.066 ;Z[2]=-4916.724;X[3]=1615.784 ;Y[3]= 1001.358 ;Z[3]= -4916.691;X[4]=1617.197 ;Y[4]= 690.040 ;Z[4]= -4917.070;X[5]=1615.933 ;Y[5]= 387.492;Z[5]= -4917.234;X[6]=1618.946 ; Y[6]= 145.491;Z[6]= -4917.807;X[7]=1616.753 ;Y[7]= -186.640 ;Z[7]= -4918.531;X[8]=1617.714 ;Y[8]= -482.264;Z[8]= -4918.561;X[9]=1617.859 ;Y[9]= -680.767;Z[9]= -4918.799;/*for(i=0;i<N;i++){printf("请输入第%d个点的像片坐标：",i+1); //---------输入点像片坐标scanf_s("%lf%lf",&x[i],&y[i]);}cout<<endl;*/x[0]=1105.312313;y[0]=3367.979242;x[1]=1111.483438;y[1]=3203.998009;x[2]=1128.936009;y[2]=2886.146374;x[3]=1145.729428;y[3]=2615.385301;x[4]=1163.579521;y[4]=2343.603584;x[5]=1179.483899;y[5]=2076.672738;x[6]=1196.664025;y[6]=1861.971885;x[7]=1216.048311;y[7]=1565.046087;x[8]=1236.196347;y[8]=1299.275413;x[9]=1249.998721;y[9]=1120.52425;/*-----------------设定外方位元素初始值--------------*/x0=48.6401;y0=11.4765;f=4548.5949;Xs=3391.658;Ys=-135.645;Zs=97.250;q=-0.0102;w=0.0620;k=-0.0887;XX[3]=1;/*------------------迭代计算--------------------------*/while((XX[3]>0.00001 || XX[4]>0.00001 || XX[5]>0.00001)&&n<100) {/*----------------旋转矩阵R-----------------------*/a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);/*-----------------像点坐标计算值------------------*/for(i=0;i<N;i++){X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs);Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs);Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);x1[i]=x0-f*X0[i]/Z0[i];y1[i]=y0-f*Y0[i]/Z0[i];}/*-------------误差方程中各偏导数的值--------------*/for(i=0;i<N;i++){A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i];A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i];A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin( k))/f+f*cos(k))*cos(w);A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;A[2*i][5]=y[i]-y0;L[2*i]=x[i]-x1[i];A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];A[1+2*i][1]=((b[1]*f+b[2]*(y[i]-y0)))/Z0[i];A[1+2*i][2]=((c[1]*f+c[2]*(y[i]-y0)))/Z0[i];A[1+2*i][3]=-(x[i]-x0)*sin(w)-((y[i]-y0)*((x[i]-x0)*cos(k)-(y[i]-y0)*sin(k))/f-f*sin(k))*cos(w);A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k ))/f;A[1+2*i][5]=-x[i]+x0;L[1+2*i]=y[i]-y1[i];}/*-------------------解法方程--------------------*/transpose(&A[0][0],&At[0][0],2*N,6);mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);invers_matrix(&result1[0][0],6);mult(&At[0][0],L,&result2[0][0],6,2*N,1);mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);Xs+=XX[0];Ys+=XX[1];Zs+=XX[2];q+=XX[3];w+=XX[4];k+=XX[5];n++;}/*----------------旋转矩阵R-----------------------*/a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);cout<<"迭代次数为:"<<n<<endl;printf("\n像片外方位元素的解\n");cout<<"航向顷角q:"<<q<<endl;cout<<"旁向倾角w:"<<w<<endl;cout<<"像片旋角k:"<<k<<endl;printf("Xs=%lf\t\tYs=%lf\t\tZs=%lf\n",Xs,Ys,Zs);cout<<endl;printf("旋转矩阵R:\n");for(i=0;i<3;i++)printf("%lf\t",a[i]);printf("\n");for(i=0;i<3;i++)printf("%lf\t",b[i]);printf("\n");for(i=0;i<3;i++)printf("%lf\t",c[i]);printf("\n");/*-------------------计算单位权中误差---------------*/ for(i=0;i<2*N;i++)sum+=L[i]*L[i];m0=sqrt(sum/(2*N-6));cout<<"单位权中误差m0="<<m0<<endl;cout<<"测量精度:"<<endl;cout<<"Xs的精度为:"<<m0*sqrt(result1[0][0])<<endl;cout<<"Ys的精度为:"<<m0*sqrt(result1[1][1])<<endl;cout<<"Zs的精度为:"<<m0*sqrt(result1[2][2])<<endl;cout<<"q的精度为:"<<m0*sqrt(result1[3][3])<<endl;cout<<"w的精度为:"<<m0*sqrt(result1[4][4])<<endl;cout<<"k的精度为:"<<m0*sqrt(result1[5][5])<<endl;system("pause");}四、实验结果分析五、实验心得通过这次实验，使我掌握了在C++及Microsoft visual studio 环境下单张像片空间后方交会算法编程及实现，掌握空间后方交会的基本原理和基本步骤，并且锻炼了我编程的能力，看到自己编出的程序成功实现后，感到很有成就感，这也是我有了更大的自信，也为以后的工作打下了更好的基础，感谢老师的这次实习。