第一篇：摄影测量空间后前交会VC++
#include
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 for(j=0;j result[i*j_2+j]=0.0; for(k=0;k result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2]; } return;} /*---------------矩阵求逆---------------------*/ void inverse(double c[n][n]){ int i,j,h,k;double p;double q[n][12];for(i=0;i for(j=0;j q[i][j]=c[i][j]; for(i=0;i for(j=n;j<12;j++){ if(i+6==j) q[i][j]=1; else q[i][j]=0; } for(h=k=0;k for(i=k+1;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 p=1.0/q[i][i]; for(j=0;j<12;j++) q[i][j]*=p;} for(i=0;i for(j=0;j c[i][j]=q[i][j+6];} /*---------------矩阵转置---------------------*/ 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 main(){ double Xs,Ys,Zs,q,w,k;double Xsr,Ysr,Zsr,qr,wr,kr;// double a[3],b[3],c[3];double ar[3],br[3],cr[3];// double x0,y0,f;double x[N],y[N];double X[N],Y[N],Z[N];double x1[N],y1[N];double xr[N],yr[N];// double x1r[N],y1r[N];// double m;double L[2*N];double Lr[2*N];// double XX[6];double XXr[6];// double A[2*N][6];double B[2*N][6];// double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];double X0r[N],Y0r[N],Z0r[N],Bt[6][2*N],result1r[6][6],result2r[6][1];// double R1[3][3],R2[3][3],XI1[3],XII1[3],XI2[3],XII2[3],XI3[3],XII3[3],XI4[3],XII4[3],XI5[3],XII5[3];double N1,N2,N3,N4,N5,BX,BY,BZ,G[5][3];double RS1[3][1],RSR1[3][1],RS2[3][1],RSR2[3][1],RS3[3][1],RSR3[3][1],RS4[3][1],RSR4[3][1],RS5[3][1],RSR5[3][1];int i,j,n=0,nr=0;double sum=0;double m1,n1,m2,n2,m3,n3,m4,n4,m5,n5,p1,q1,p2,q2,p3,q3,p4,q4,p5,q5; /*---------------输入点地面坐标---------------------*/ X[0]=5083.205;X[1]=5780.02;X[2]=5210.879;X[3]=5909.264;Y[0]=5852.099;Y[1]=5906.365;Y[2]=4258.446;Y[3]=4314.283;Z[0]=527.925;Z[1]=571.549;Z[2]=461.81;Z[3]=455.48;m1=0.051758;n1=0.081555;p1=-0.039953;q1=0.078463;m2=0.014618;n2=-0.000231;p2=-0.076016;q2=0.000036;m3=0.04988;n3=-0.000792;p3=-0.042201;q3=-0.001022;m4=0.086243;n4=-0.001346;p4=-0.007706;q4=-0.002112;m5=0.048135;n5=-0.079962;p5=-0.044438;q5=-0.079736;/*---------------输入点像片坐标---------------------*/ x[0]=16.012;x[1]=88.56;x[2]=13.362;x[3]=82.24;y[0]=79.963;y[1]=81.134;y[2]=-79.37;y[3]=-80.027; xr[0]=-73.93;xr[1]=-5.252;xr[2]=-79.122;xr[3]=-9.887;yr[0]=78.706;yr[1]=78.184;yr[2]=-78.879;yr[3]=-80.089; /*-----------------设定外方位元素初始值--------------*/ x0=0;y0=0;f=152.00;m=10000;Xs=0;Ys=0;Zs=f*m/1000;Xsr=0;Ysr=0;Zsr=f*m/1000;// q=0;w=0;k=0;qr=0;wr=0;kr=0;// XX[3]=1;XXr[3]=1;// /*------------------迭代计算左片Xs,Ys,Zs,q,w,k--------------------------*/ while((XX[3]>6/206265 || XX[4]>6/206265 || XX[5]>6/206265)&&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 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);//求A转置矩阵 mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);//A矩阵 X A转置矩阵-> result1 inverse(result1);// result1 求逆 mult(&At[0][0],L,&result2[0][0],6,2*N,1);// A转置矩阵 X L-> result2 mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);// result1 X result2-> XX Xs+=XX[0]; } Ys+=XX[1];Zs+=XX[2];q+=XX[3];w+=XX[4];k+=XX[5];n++;/*------------------迭代计算有片Xsr,Ysr,Zsr,qr,wr,kr--------------------------*/ while((XXr[3]>6/206265 || XXr[4]>6/206265 || XXr[5]>6/206265)&&nr<100){//右点算Xs,Ys,Zs,q,w,k /*----------------旋转矩阵R-----------------------*/ ar[0]=cos(qr)*cos(kr)-sin(qr)*sin(wr)*sin(kr); ar[1]=-cos(qr)*sin(kr)-sin(qr)*sin(wr)*cos(kr); ar[2]=-sin(qr)*cos(wr); br[0]=cos(wr)*sin(kr); br[1]=cos(wr)*cos(kr); br[2]=-sin(wr); cr[0]=sin(qr)*cos(kr)+cos(qr)*sin(wr)*sin(kr); cr[1]=-sin(qr)*sin(kr)+cos(qr)*sin(wr)*cos(kr); cr[2]=cos(qr)*cos(wr); /*-----------------像点坐标计算值------------------*/ for(i=0;i X0r[i]=ar[0]*(X[i]-Xsr)+br[0]*(Y[i]-Ysr)+cr[0]*(Z[i]-Zsr); Y0r[i]=ar[1]*(X[i]-Xsr)+br[1]*(Y[i]-Ysr)+cr[1]*(Z[i]-Zsr); Z0r[i]=ar[2]*(X[i]-Xsr)+br[2]*(Y[i]-Ysr)+cr[2]*(Z[i]-Zsr); x1r[i]=x0-f*X0r[i]/Z0r[i]; y1r[i]=y0-f*Y0r[i]/Z0r[i]; } /*-------------误差方程中各偏导数的值--------------*/ for(i=0;i B[2*i][0]=((ar[0]*f+ar[2]*(xr[i]-x0)))/Z0r[i]; B[2*i][1]=((br[0]*f+br[2]*(xr[i]-x0)))/Z0r[i]; B[2*i][2]=((cr[0]*f+cr[2]*(xr[i]-x0)))/Z0r[i]; B[2*i][3]=(yr[i]-y0)*sin(wr)-((xr[i]-x0)*((xr[i]-x0)*cos(kr)-yr[i]*sin(kr))/f+f*cos(kr))*cos(wr); B[2*i][4]=-f*sin(kr)-(xr[i]-x0)*((xr[i]-x0)*sin(kr)+(yr[i]-y0)*cos(kr))/f; B[2*i][5]=yr[i]-y0; Lr[2*i]=xr[i]-x1r[i]; B[1+2*i][0]=((ar[1]*f+ar[2]*(yr[i]-y0)))/Z0r[i]; B[1+2*i][1]=((br[1]*f+br[2]*(yr[i]-y0)))/Z0r[i]; B[1+2*i][2]=((cr[1]*f+cr[2]*(yr[i]-y0)))/Z0r[i]; B[1+2*i][3]=-(xr[i]-x0)*sin(wr)-((yr[i]-y0)*((xr[i]-x0)*cos(kr)-(yr[i]-y0)*sin(kr))/f-f*sin(kr))*cos(wr); } } B[1+2*i][4]=-f*cos(kr)-(yr[i]-y0)*((xr[i]-x0)*sin(kr)+(yr[i]-y0)*cos(kr))/f;B[1+2*i][5]=-xr[i]+x0;Lr[1+2*i]=yr[i]-y1r[i];/*-------------------解法方程--------------------*/ transpose(&B[0][0],&Bt[0][0],2*N,6);//求A转置矩阵 mult(&Bt[0][0],&B[0][0],&result1r[0][0],6,2*N,6);//A矩阵 X A转置矩阵-> result1 inverse(result1r);// result1 求逆 mult(&Bt[0][0],Lr,&result2r[0][0],6,2*N,1);// A转置矩阵 X L-> result2 mult(&result1r[0][0],&result2r[0][0],&XXr[0],6,6,1);// result1 X result2-> XX Xsr+=XXr[0];Ysr+=XXr[1];Zsr+=XXr[2];qr+=XXr[3];wr+=XXr[4];kr+=XXr[5];nr++;BX=Xsr-Xs;//BX,BY,BZ BY=Ysr-Ys;BZ=Zsr-Zs;for(i=0;i<3;i++){//R1,R2赋值 R1[0][i]=a[i]; R1[1][i]=b[i]; R1[2][i]=c[i]; R2[0][i]=ar[i]; R2[1][i]=br[i]; R2[2][i]=cr[i];} RS1[0][0]=m1;RS1[1][0]=n1;RS1[2][0]=-f/1000;// 给 将要与R1，R2相乘的 RS矩阵和RSR矩阵 赋值 RS2[0][0]=m2;RS2[1][0]=n2;RS2[2][0]=-f/1000;RS3[0][0]=m3;RS3[1][0]=n3;RS3[2][0]=-f/1000;RS4[0][0]=m4;RS4[1][0]=n4;RS4[2][0]=-f/1000;RS5[0][0]=m5;RS5[1][0]=n5;RS5[2][0]=-f/1000; RSR1[0][0]=p1;RSR1[1][0]=q1;RSR1[2][0]=-f/1000;RSR2[0][0]=p2;RSR2[1][0]=q2;RSR2[2][0]=-f/1000;RSR3[0][0]=p3;RSR3[1][0]=q3;RSR3[2][0]=-f/1000;RSR4[0][0]=p4;RSR4[1][0]=q4;RSR4[2][0]=-f/1000;RSR5[0][0]=p5;RSR5[1][0]=q5;RSR5[2][0]=-f/1000; mult(&R1[0][0],&RS1[0][0],&XI1[0],3,3,1);//R1 X RS mult(&R1[0][0],&RS2[0][0],&XI2[0],3,3,1);mult(&R1[0][0],&RS3[0][0],&XI3[0],3,3,1);mult(&R1[0][0],&RS4[0][0],&XI4[0],3,3,1);mult(&R1[0][0],&RS5[0][0],&XI5[0],3,3,1);mult(&R2[0][0],&RSR1[0][0],&XII1[0],3,3,1);//R2 X RSR mult(&R2[0][0],&RSR2[0][0],&XII2[0],3,3,1);mult(&R2[0][0],&RSR3[0][0],&XII3[0],3,3,1);mult(&R2[0][0],&RSR4[0][0],&XII4[0],3,3,1);mult(&R2[0][0],&RSR5[0][0],&XII5[0],3,3,1);N1=(BX*XII1[2]-BZ*XII1[0])/(XI1[0]*XII1[2]-XI1[2]*XII1[0]);//给N1 赋值 有五个点 N2=(BX*XII2[2]-BZ*XII2[0])/(XI2[0]*XII2[2]-XI2[2]*XII2[0]);N3=(BX*XII3[2]-BZ*XII3[0])/(XI3[0]*XII3[2]-XI3[2]*XII3[0]);N4=(BX*XII4[2]-BZ*XII4[0])/(XI4[0]*XII4[2]-XI4[2]*XII4[0]);N5=(BX*XII5[2]-BZ*XII5[0])/(XI5[0]*XII5[2]-XI5[2]*XII5[0]);G[0][0]=Xs+N1*XI1[0];//五个点XA,YA,ZA计算G[i][j] i:第几个点，j：是x，y还是z G[0][1]=Ys+N1*XI1[1];G[0][2]=Zs+N1*XI1[2];G[1][0]=Xs+N2*XI1[0];G[1][1]=Ys+N2*XI1[1];G[1][2]=Zs+N2*XI1[2];G[2][0]=Xs+N3*XI1[0];G[2][1]=Ys+N3*XI1[1];G[2][2]=Zs+N3*XI1[2];G[3][0]=Xs+N4*XI1[0];G[3][1]=Ys+N4*XI1[1];G[3][2]=Zs+N4*XI1[2];G[4][0]=Xs+N5*XI1[0];G[4][1]=Ys+N5*XI1[1];G[4][2]=Zs+N5*XI1[2];/*----------------旋转矩阵R-----------------------*/ cout<<“--------左--------”< ”< cout<<“Ys ”< cout<<“Zs ”< cout<<“迭代次数为:”< ”< cout<<“Ysr ”< cout<<“Zsr ”< 5”< %-.5f %.5f %-.5f %.5f %-.5f n”,XI1[0],XI2[0],XI3[0],XI4[0],XI5[0]);printf(“ Y1 %-.5f %.5f %-.5f %.5f %-.5f n”,XI1[1],XI2[1],XI3[1],XI4[1],XI5[1]);printf(“ Z1 %-.5f %.5f %-.5f %.5f %-.5f n”,XI1[2],XI2[2],XI3[2],XI4[2],XI5[2]);printf(“ X2 %-.5f %.5f %-.5f %.5f %-.5f n”,XII1[0],XII2[0],XII3[0],XII4[0],XII5[0]);printf(“ Y2 %-.5f %.5f %-.5f %.5f %-.5f n”,XII1[1],XII2[1],XII3[1],XII4[1],XII5[1]);printf(“ Z2 %-.5f %.5f %-.5f %.5f %-.5f n”,XII1[2],XII2[2],XII3[2],XII4[2],XII5[2]); printf(“ N1 %-.5f %.5f %-.5f %.5f %-.5f nn”,N1,N2,N3,N4,N5);printf(“ XA %-.5f %.5f %-.5f %.5f %-.5f n”,G[0][0],G[1][0],G[2][0],G[3][0],G[4][0]);printf(“ YA %-.5f %.5f %-.5f %.5f %-.5f n”,G[0][1],G[1][1],G[2][1],G[3][1],G[4][1]);printf(“ ZA %-.5f %.5f %-.5f %.5f %-.5f n”,G[0][2],G[1][2],G[2][2],G[3][2],G[4][2]); } 摄影测量实习报告 实习内容:单片空间后方交会编程 实习者:李友兵 学号:0810050121 指导老师:张金平老师 实习时间:2011.05.30——2011.06.03 一、实习目的与任务 此次摄影测量实习主要是要自主编程实现单像空间后方交会，通过已知的内方位元素和控制点像点坐标和地面坐标求解六个外方位元素，在此过程中深入理解单像空间后方交会的原理和对编程的熟悉和理解（我用的是C语言编程），和对时间的合理运用，对知识的综合运用，培养理论的实际运用能力，任务是在一个星期内自主完成。 二、单片空间后方交会理论基础 单像空间后方交会：是通过以像点平面坐标为观测值，以控制点坐标为已知值，利用共线条件方程和最小二乘原理，运用间接平差方法，通过迭代求解6个外方位元素。程序设计的思路与流程 三、程序设计的思路与流程（编程框架和步骤） 1、根据内方位元素和已知的数据将控制点的框标坐标转换为像平面坐标系坐标：x=x′-x0，y=y′-y0 2、确定未知数的初始值：角元素初始值κ0=ω0=Φ0=0；线元素初始值Zs0=H=mf，Xs0=(X1+X2+X3+X4)/4,Ys0=（Y1+Y2+Y3+Y4）/4 3、利用角元素初始值计算方向余弦值组成旋转矩阵R a1=cosΦ*cosκ-sinΦ*sinω*sinκ,a2=-cosΦ*sinκ-sinΦ*sinω*cosκ，a3=-sinΦ*cosω b1=cosω*sinκ,b2=cosω*cosκ,b3=-sinω c1=sinΦ*cosκ+cosΦ*sinω*sinκ,c2=-sinΦ*sinκ+cosΦ*sinω*cosκ,c3=cosΦ*cosω 4、逐点计算控制点的像点坐标的近似值，共线方程为: x=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)) y=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)) (x)=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs))(y)=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)) 5、组成误差方程：Vx=a11*dXs+a12*dYs+a13*dZs+a14*dΦ+a15*dω+a16*dκ+(x)-x Vy=a21*dXs+a22*dYs+a23*dZs+a24*dΦ+a25*dω+a26*dκ+(y)-y 系数求解(是共线方程分别外方位元素求导，是共线方程线性化的系数)：变量代换 A= a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)B= a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)C= a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)a11=(a1*f+a3*x)/C,a12=(b1*f+b3*x)/C,a13=(c1*f+c3*f)/C,a14=y*sinω-((x*(x*cosκ-y*sinκ))/f+f*cosκ)cosω,a15=-f*sinκ-(x*(x*sinκ+y*cosκ)/f),a16=y a21=(a2*f+a3*x)/C,a22=(b2*f+b3*x)/C,a23=（c2*f+c3*x）/C,a24=-xsinω-((x*(x*cosκ-y*sinκ))/f-f*sinκ)cosω,a25=-f*cosκ-(y*(x*sinκ+y*cosκ)/f),a26 =-x 误差方程的常系数（是像点坐标观测值与计算的近似值的差值）： lx=x-(x),ly=y-(y) 6、组成法方程，解求外方位元素改正数X=(ATA)-1ATL(A为误差方程的系数矩阵，L为误差方程的常系数矩阵，通过步骤5求得，此处先求ATA再求矩阵的逆矩阵，解得的改正数加上相应的近似值得到外方位元素新的近似值) 7、检查计算是否收敛：将求得的外方位元素改正数与规定的限差比较，大于限差继续迭代，小于限差则终止。 四、各子函数详细设计的关键技术参数 子函数（输入函数、Input，矩阵求积Matrixmultiply，计算函数Resection，矩阵转置Matrixtranspose，矩阵求逆Matrixinverse，输出函数Output）主要用了Matrixtranspose，矩阵的行变列，列变行,参数为要转置的矩阵，转置后的矩阵，要转置矩阵的行列数，Matrixinverse,求矩阵的代数余子式，参数有要求逆的矩阵和，逆矩阵的行数，Matrixmultiply，一矩阵的行乘以二矩阵的列，参数为一矩阵，二矩阵，所求矩阵，一的行，一的列，二的列 五、像片外方位元素解算结 六、实习体会 从周一开始进行整个实习框架进行构建，编写了实习的思路既步骤；然后是梳理框架，要用到的子函数，子函数参数，怎么编写，有不懂的东西网上参考资料；再是自己动手编程，发现编程要求很细腻，出不得一点差错，程序编写出来后，有很多的错误，要求逐一进行改正，修改；最后才得以运行。通过本次实习，我深刻理解了单片空间后方交会原理，进一步熟悉理解了C语言编程，认识到了时间的搭配的重要性和参考资料的必要性，当遇到难题是要迎难而上，达到突破，当完成是能感到一丝丝欣慰和成就感。 附录:源代码 #include “stdio.h” #include “math.h” #include “Matrixmultiply.c” #include “Matrixtranspose.c” #include “Matrixinverse.c” void main(){ int i,j,k,f=0;double x0=0.00018, y0=0.00026,fk=0.15324; //内方位元素 double m=40000;//估算比例尺 double R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];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]; double B[4][5]={-0.08615,-0.06899,36589.41,25273.32,2195.17,-0.05340,0.08221,37631.08,31324.51,728.69,-0.01478,-0.07663,39100.97,24934.98,2386.80,0.01046,0.06443,40426.54,30319.81,757.31}; for(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=m*fk;//求得外方位线元素的初始值 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(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++; } Matrixtranspose(A,AT,8,6);//此为转置函数的调用,求AT Matrixmultiply(AT,A,ATA,6,8,6);//此为矩阵相乘函数的调用，求ATA Matrixinverse(ATA,6);//此为求矩阵逆函数的调用，求ATA的逆 Matrixmultiply(AT,L,ATL,6,8,1);//此为矩阵相乘函数的调用，求ATL Matrixmultiply(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]>=0.00000001||XG[4][0]>=0.00000001||XG[5][0]>=0.00000001);//当限差满足要求时要再一次进行旋转矩阵的求解 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);printf(“迭代次数:%d”,f); //屏幕输出误差方程系数阵、常数项、改正数 printf(“nn误差方程系数矩阵A为:nn”);for(i=0;i<6;i++){ for(j=0;j<6;j++) printf(“%13.5e ”,A[i][j]); printf(“n”);} printf(“n常数项L为:nn”); for(i=0;i<8;i++){ for(j=0;j<1;j++) printf(“%13.5e ”,L[i][j]); printf(“n”);} printf(“n改正数XG为:nn”); for(i=0;i<6;i++){ for(j=0;j<1;j++) printf(“%13.5e ”,XG[i][j]); printf(“n”);} printf(“n相片的外方位元素为:nn”); printf(“ Xs=%13.7e, Ys=%13.7e, Zs=%13.7e nn”,Xs,Ys,Zs); printf(“ Q=%13.7e, W=%13.7e, K=%13.7e n”,Q,W,K);printf(“n旋转矩阵R为:nn”); for(i=0;i<3;i++){ for(j=0;j<3;j++) printf(“%13.5e ”,R[i][j]); printf(“n”);} } //子函数 #include #include #include int Matrixinverse(a,n) int n; double a[]; { int *is,*js,i,j,k,l,u,v; double d,p; is=malloc(n*sizeof(int)); js=malloc(n*sizeof(int)); for(k=0;k<=n-1;k++) { d=0.0; for(i=k;i<=n-1;i++) for(j=k;j<=n-1;j++) { l=i*n+j;p=fabs(a[l]); if(p>d){ d=p;is[k]=i;js[k]=j;} } if(d+1.0==1.0) { free(is);free(js);printf(“err**not invn”); return(0); } if(is[k]!=k) for(j=0;j<=n-1;j++) { u=k*n+j;v=is[k]*n+j; p=a[u];a[u]=a[v];a[v]=p; } if(js[k]!=k) for(i=0;i<=n-1;i++) { u=i*n+k;v=i*n+js[k]; p=a[u];a[u]=a[v];a[v]=p; } l=k*n+k; a[l]=1.0/a[l]; for(j=0;j<=n-1;j++) if(j!=k) { u=k*n+j;a[u]=a[u]*a[l];} for(i=0;i<=n-1;i++) if(i!=k) for(j=0;j<=n-1;j++) if(j!=k) { u=i*n+j; a[u]=a[u]-a[i*n+k]*a[k*n+j]; } for(i=0;i<=n-1;i++) if(i!=k) { u=i*n+k;a[u]=-a[u]*a[l];} } for(k=n-1;k>=0;k--) { if(js[k]!=k) for(j=0;j<=n-1;j++) { u=k*n+j;v=js[k]*n+j; p=a[u];a[u]=a[v];a[v]=p; } if(is[k]!=k) for(i=0;i<=n-1;i++) { u=i*n+k;v=i*n+is[k]; p=a[u];a[u]=a[v];a[v]=p; } } free(is);free(js); return(1); } //子函数 void Matrixmultiply(a,b,c,m,n,k)int m,n,k;double a[],b[],c[];{ int i,j,l,u;for(i=0;i for(j=0;j { u=i*k+j;c[u]=0.0; for(l=0;l c[u]+=a[i*n+l]*b[l*k+j]; } return;} //子函数 void Matrixtranspose(a,b,m,n)int m,n;double a[],b[];{ int i,j,u;for(i=0;i for(j=0;j { u=j*m+i;b[u]=0.0;b[j*m+i]=a[i*n+j]; } } return;} 4D指的是DEM、DOM、DLG、DRG。意义如下： 数字高程模型（Digital Elevation Model 简称DEM）是在高斯投影平面上规则格网点平面坐标（x,y）及其高程（z）的数据集。Dem的水平间隔可随地貌类型不同而改变。根据不同的高程精度，可分为不同等级产品。 数字正射影像图（Digital Orthophoto Map简称DOM）是利用数字高程模型对扫描处理的数字化的航空相片 / 遥感相片（单色 / 彩色），经逐象元进行纠正，再按影像镶嵌，根据图幅范围剪裁生成的影像数据。一般带有公里格网、图廓内 / 外整饰和注记的平面图。 数字线划地图（Digital Line Graphic简称DLG）是现有地形图上基础地理要素的矢量数据集，且保存要素间空间关系和相关的属性信息。 数字栅格地图（Digital Raster Graphic简称DRG）是纸质地形图的数字化产品。每幅图经扫描、纠正、图幅处理及数据压缩处理后，形成在内容、几何精度和色彩上与地形图保持一致的栅格文件。 像点位移：当像片倾斜、地面起伏时，地面点在航摄像片上构像相对于理想情况下的构像所产生的位置差异称像点位移 引起原因：1.像片倾斜引起的像点位移 2.地形起伏引起的像点位移 像片的内方位元素：摄影物镜后节点与像片之间相互位置的参数 像片外方位元素：已建立的摄影光束，确定像片摄影瞬间在地面直角坐标系中空间位置和姿态的参数 像平面坐标系:用以表示像点在像平面上的位置，通常采用右手坐标系x,y轴的选择按需要而定，在解析和数字摄影测量中，常根据框标来确定像平面坐标系，称为像框标坐标系。像空间坐标系:以摄影中心S为坐标原点，x，y轴与像平面坐标系的x，y轴平行，z轴与主光轴重合，形成像空间右手直角坐标系S-xyz。 像空间辅助坐标系:此坐标系的原点仍选在摄影中心S，坐标轴系的选择视需要而定，通常有三种选取方法。其一是取铅垂方向为z轴，航向为X轴，构成右手直角坐标系，见图(a)。其二是以每条航线内第一张像片的像空间坐标系作为像空间辅助坐标系，见图(b)。其三是以每个像片对的左片摄影中心为坐标原点，摄影基线方向为X轴，以摄影基线及左片主光轴构成的面作为XZ平面，构成右手直角坐标系，如图(c)。用S-XYZ表示。 地面测量坐标系：空间大地坐标基准下的高斯-克吕克6°带或3°带投影的平面直角坐标与定义的从某一基准面量起的高程，两者结合而成的空间左手直角坐标系。 摄影测量坐标系：是指描述摄影测量模型的空间直角坐标系。其原点选在某摄站或某一已知点，X轴大体与航线方向一致，Z轴与铅垂线方向一致且向上为正的右旋空间直角坐标系。像点位移：航空像片是地面的中心投影，根据中心投影的原理，无论是带有起伏状态的地形，还是高出地面的任何物体，反映到航空像片上的像点与其平面位置相比，一般都会产生位置的移动，这种像点位置的移动，叫做像点位移。主要是由像片倾斜、地面点相对于基准面的高差和物理因素（如摄影材料变形、压平误差、摄影物镜畸变、大气折光和地球曲率等）产生。 为什么研究像点位移的规律可以清楚知道透镜成像的大小，虚实6-=等基本情况，对于应用透镜解决生活和实际问题是必不可少的像点位移的规律： 单向空间后方交会：已知至少3个地面控制点的坐标A，B，C，与其影像上对应的三个像点的影像坐标a，b，c，根据共线方程，反求该像片的6个外方位元素。 立体像对双像前方交会：现已知这两张像片的内外方位元素，设想将该像片的内外方位元素值置于摄影时的位置，显然同名射线S1a1和S2a2必然交于地面点A。这种由立体像对中两张像片的内外方位元素和像点坐标来确定相应地面点的地面坐标的方法。 立体像对双像前方交会的目的：利用立体像对上的同名像点，才能得到两条同名射线在空间相交的点。 解析法相对定向：通过计算相对定向元素建立地面立体模型。 解析法相对定向的共面条件：B·（S1a1 X S2a2）=0 模型连接：将单个模型连接成为航带模型,要将各模型的不同的比例尺归化为统一的比例尺,通常以相邻像对重叠范围内三个连接点的高程应相等为条件,建立统一的以第一个模型的比例尺为基准的航带模型.模型连接的实质是求出相邻模型间的比例归化系数.航带网整体平差的实质是以一条航带模型为平差单元,解求航带的非线性改正系数,即多项式系数.第二篇：摄影测量实习报告-单片空间后方交会
第三篇：摄影测量
第四篇：摄影测量
第五篇：摄影测量
点击咨询 