1、一 實(shí)驗(yàn)?zāi)康模?掌握攝影測量空間前方交會(huì)的原理，利用計(jì)算機(jī)編程語言實(shí)現(xiàn)空間前方交會(huì)外方位元素的解算。二 儀器用具及數(shù)據(jù)文件： 計(jì)算機(jī)windows xp 系統(tǒng)，編程軟件VISUAL C+6.0，地面控制點(diǎn)在攝影測量坐標(biāo)系中的坐標(biāo)及其像點(diǎn)坐標(biāo)文件shuju.txt。 三 實(shí)驗(yàn)內(nèi)容： 單張影像的空間前方交會(huì)：利用地面控制點(diǎn)數(shù)據(jù)及相應(yīng)像點(diǎn)坐標(biāo)根據(jù)共線方程反求影像的外方位元素。 數(shù)學(xué)模型：共線條件方程式： 求解過程：1獲取數(shù)據(jù)。從航攝資料中查取平均航高與攝影機(jī)主距；獲取控制點(diǎn)的地面測量坐標(biāo)并轉(zhuǎn)換為地面攝影測量坐標(biāo)。2量測控制點(diǎn)的像點(diǎn)坐標(biāo)并做系統(tǒng)改正。3確定未知數(shù)的初始值。在豎直攝影且地面控制點(diǎn)大致分布

2、均勻的情況下，按如下方法確定初始值，即：, =0 式中；m為攝影比例尺分母；n為控制點(diǎn)個(gè)數(shù)。 4用三個(gè)角元素的初始值，計(jì)算個(gè)方向余弦，組成旋轉(zhuǎn)矩陣R。 5逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值。利用未知數(shù)的近似值和控制點(diǎn)的地面坐標(biāo)代入共線方程式，逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值x、y。 6逐點(diǎn)計(jì)算誤差方程式的系數(shù)和常數(shù)項(xiàng)，組成誤差方程式。7計(jì)算法方程的系數(shù)矩陣和常數(shù)項(xiàng)，組成法方程式。 8解法方程，求得外方位元素的改正數(shù),，d,d,d。9用前次迭代取得的近似值，加本次迭代的改正數(shù)，計(jì)算外方位元素的新值。10將求得的外方位元素改正數(shù)與規(guī)定的限差比擬，假設(shè)小于限差那么迭代結(jié)束。否那么用新的近似值重復(fù)49，直到滿足要求為止

3、。四實(shí)驗(yàn)結(jié)果：程序的源代碼如下所示：#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<math.h>#include<conio.h>#define N 4void turn(double *A,double A2,int m,int n) /計(jì)算矩陣的轉(zhuǎn)置int i,j; for(i=0;i<m;i+)for(j=0;j<n;j+)A2j*m+i=Ai*n+j;void mulAB(double *A,double *B,double *C,int

4、 am,int an,int bm,int bn) /計(jì)算兩矩陣相乘int i,j,l,u;if(an!=bm)printf("error!cannot do the multiplication.n");return;for(i=0;i<am;i+)for(j=0;j<bn;j+)u=i*bn+j;Cu=0.0;for(l=0;l<an;l+)Cu+=Ai*an+l*Bl*bn+j;return;double *inv(double *a,int n) /計(jì)算矩陣的逆，本程序的難點(diǎn) /采用高斯-約旦-全選主元法int *is,*js,i,j,k,l,u,

5、v; double d,p; is=(int*)malloc(n*sizeof(int); js=(int*)malloc(n*sizeof(int); for (k=0; k<=n-1; k+) d=0.0; for (i=k;i<n;i+) for (j=k;j<n;j+) l=i*n+j; p=fabs(al); if (p>d) d=p; isk=i; jsk=j; if (d+1.0=1.0) free(is); free(js); printf("error not invn"); return NULL; if (isk!=k) for

6、 (j=0;j<n;j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0;i<n;i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/al; for (j=0;j<n;j+) if (j!=k) u=k*n+j; au=au*al; for (i=0;i<n;i+) if (i!=k) for (j=0;j<n;j+) if (j!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0;i<

7、n;i+) if (i!=k) u=i*n+k; au=-au*al; for (k=n-1;k>=0;k-) if (jsk!=k) for (j=0;j<=n-1;j+) u=k*n+j; v=jsk*n+j; p=au; au=av; av=p; if (isk!=k) for (i=0;i<n;i+) u=i*n+k; v=i*n+isk; p=au; au=av; av=p; free(is);free(js); return a;void printmatrix(double M,int m,int n) /矩陣的輸出int i,j;for(i=0;i<m;

8、i+)for(j=0;j<n;j+)printf("%.5f",Mi*n+j);printf("n");printf("n");main() /主函數(shù)，空間前方交會(huì)的計(jì)算FILE *fp; /定義一個(gè)文件指針fpint m,i,j=0;double f,t,w,k,S1=0.0,S2=0.0,S3=0.0,xN,yN,x0N,y0N,XN,YN,ZN,Xs0,Ys0,Zs0;double a3,b3,c3,A2*N*6,AT6*2*N,ATA6*6,*ATA_=NULL,l2*N,ATl6,V6; if(fp=fopen(&qu

9、ot;e:shuju.txt","r")=NULL) /使fp指向被翻開的shuju.txt文件 /并判斷文件是否翻開成功printf("nerror on open shuju.txtn");getch();exit(1); for(i=0;i<N;i+)fscanf(fp,"%lf%lf%lf%lf%lf",&xi,&yi,&Xi,&Yi,&Zi); /將文件中的數(shù)據(jù)賦給主函數(shù)定義的變量printf("原始數(shù)據(jù):n");printf("xttyt

10、tXttYttZttn"); /輸出文件中的原始數(shù)據(jù)for(i=0;i<N;i+)printf("%lf %lf %lf %lf %lfn",xi,yi,Xi,Yi,Zi);printf("n請(qǐng)輸入攝影機(jī)主距和攝影比例尺分母；f,m:");scanf("%lf,%lf",&f,&m); /輸入f,mf=f/1000.0;for(i=0;i<N;i+)xi/=1000.0;yi/=1000.0;S1+=Xi;S2+=Yi;S3+=Zi; Xs0=S1/N;Ys0=S2/N; /計(jì)算外方位元素的初始值

11、Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; while(j<3) printf("-第%d次計(jì)算-n",j+1);a0=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);a1=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k);a2=-sin(t)*cos(w);b0=cos(w)*sin(k);b1=cos(w)*cos(k); /計(jì)算旋轉(zhuǎn)矩陣b2=-sin(w);c0=sin(t)*cos(k)+cos(t)*sin(w)*sin(k);c1=-sin(t)*sin(k)+cos(t)*sin(w)*co

12、s(k);c2=cos(t)*cos(w);for(i=0;i<N;i+)x0i=-f*(a0*(Xi-Xs0)+b0*(Yi-Ys0)+c0*(Zi-Zs0)/(a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0); /計(jì)算像點(diǎn)坐標(biāo)近似值y0i=-f*(a1*(Xi-Xs0)+b1*(Yi-Ys0)+c1*(Zi-Zs0)/(a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0);Gi=a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0);for(i=0;i<N;i+)Ai*12+0=(a0*f+a2*xi)/Gi;Ai*12+

13、1=(b0*f+b2*xi)/Gi;Ai*12+2=(c0*f+c2*xi)/Gi;Ai*12+3=yi*sin(w)-(xi*(xi*cos(k)-yi*sin(k)/f+f*cos(k)*cos(w); /計(jì)算誤差方程的系數(shù)陣以及l(fā)Ai*12+4=-f*sin(k)-xi*(xi*sin(k)+yi*cos(k)/f;Ai*12+5=yi;Ai*12+6=(a1*f+a2*yi)/Gi;Ai*12+7=(b1*f+b2*yi)/Gi;Ai*12+8=(c1*f+c2*yi)/Gi;Ai*12+9=-xi*sin(w)-(yi*(xi*cos(k)-yi*sin(k)/f-f*sin(k)*

14、cos(w);Ai*12+10=-f*cos(k)-yi*(xi*sin(k)+yi*cos(k)/f;Ai*12+11=-xi;li*2+0=xi-x0i;li*2+1=yi-y0i;/printf("output matrix: An");/printmatrix(A,2*N,6);/printf("output matrix: ln");/ printmatrix(l,2*N,1);turn(A,AT,2*N,6);/printf("output matrix: ATn");/printmatrix(AT,6,2*N);mulA

15、B(AT,A,ATA,6,2*N,2*N,6); /間接平差法計(jì)算外方位元素,中間計(jì)算的矩陣可以根據(jù)需要/printf("output matrix: ATAn"); /選擇性的輸出，這里將其屏蔽，為了在打印輸出結(jié)果的時(shí)候/printmatrix(ATA,6,6); /節(jié)約資源ATA_=inv(ATA,6);/printf("output matrix: ATA_n");/printmatrix(ATA_,6,6);mulAB(AT,l,ATl,6,2*N,2*N,1);/printf("outpit matrinx: ATln");

16、/printmatrix(ATl,6,1);mulAB(ATA_,ATl,V,6,6,6,1);/printf("output matrix: Vn");/printmatrix(V,6,1);Xs0+=V0;Ys0+=V1;Zs0+=V2;t+=V3;w+=V4;k+=V5;printf("第%d次計(jì)算的外方位元素為：n",+j); printf("Xs=%.5lf,Ys=%.5lf,Zs=%.5lf,t=%.5lf,w=%.5lf,k=%.5lfn",Xs0,Ys0,Zs0,t,w,k); printf("n外方位元素為：n"); printf("Xs=%.5