空間后方交會程序

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

2、豎直攝影且地面控制點大致分布均勻的情況下，按如下方法確定初始值，即：, =0 式中；m為攝影比例尺分母；n為控制點個數(shù)。 （4）用三個角元素的初始值，計算個方向余弦，組成旋轉(zhuǎn)矩陣R。 （5）逐點計算像點坐標的近似值。利用未知數(shù)的近似值和控制點的地面坐標代入共線方程式，逐點計算像點坐標的近似值（x）、（y）。 （6）逐點計算誤差方程式的系數(shù)和常數(shù)項，組成誤差方程式。（7）計算法方程的系數(shù)矩陣和常數(shù)項，組成法方程式。 （8）解法方程，求得外方位元素的改正數(shù),，d,d,d。（9）用前次迭代取得的近似值，加本次迭代的改正數(shù)，計算外方位元素的新值。（10）將求得的外方位元素改正數(shù)與規(guī)定的限差比較，若小于

3、限差則迭代結(jié)束。否則用新的近似值重復（4）（9），直到滿足要求為止。四實驗結(jié)果：程序的源代碼如下所示：#include#include#include#include#include#define N 4void turn(double *A,double A2,int m,int n) /計算矩陣的轉(zhuǎn)置int i,j; for(i=0;im;i+)for(j=0;jn;j+)A2j*m+i=Ai*n+j;void mulAB(double *A,double *B,double *C,int am,int an,int bm,int bn) /計算兩矩陣相乘int i,j,l,u;if(an

4、!=bm)printf(error!cannot do the multiplication.n);return;for(i=0;iam;i+)for(j=0;jbn;j+)u=i*bn+j;Cu=0.0;for(l=0;lan;l+)Cu+=Ai*an+l*Bl*bn+j;return;double *inv(double *a,int n) /計算矩陣的逆，本程序的難點 /采用高斯-約旦-全選主元法int *is,*js,i,j,k,l,u,v; double d,p; is=(int*)malloc(n*sizeof(int); js=(int*)malloc(n*sizeof(int)

5、; for (k=0; k=n-1; k+) d=0.0; for (i=k;in;i+) for (j=k;jd) 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 (j=0;jn;j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0;in;i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/

6、al; for (j=0;jn;j+) if (j!=k) u=k*n+j; au=au*al; for (i=0;in;i+) if (i!=k) for (j=0;jn;j+) if (j!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0;i=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;in;i+) u=i*n+k; v=i*n+isk; p=au; au=av; av=p; free(is);free(j

7、s); return a;void printmatrix(double M,int m,int n) /矩陣的輸出int i,j;for(i=0;im;i+)for(j=0;jn;j+)printf(%.5f,Mi*n+j);printf(n);printf(n);main() /主函數(shù)，空間后方交會的計算FILE *fp; /定義一個文件指針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,*AT

8、A_=NULL,l2*N,ATl6,V6; if(fp=fopen(e:shuju.txt,r)=NULL) /使fp指向被打開的shuju.txt文件 /并判斷文件是否打開成功printf(nerror on open shuju.txtn);getch();exit(1); for(i=0;iN;i+)fscanf(fp,%lf%lf%lf%lf%lf,&xi,&yi,&Xi,&Yi,&Zi); /將文件中的數(shù)據(jù)賦給主函數(shù)定義的變量printf(原始數(shù)據(jù):n);printf(xttyttXttYttZttn); /輸出文件中的原始數(shù)據(jù)for(i=0;iN;i+)printf(%lf %lf

9、 %lf %lf %lfn,xi,yi,Xi,Yi,Zi);printf(n請輸入攝影機主距和攝影比例尺分母；f,m:);scanf(%lf,%lf,&f,&m); /輸入f,mf=f/1000.0;for(i=0;iN;i+)xi/=1000.0;yi/=1000.0;S1+=Xi;S2+=Yi;S3+=Zi; Xs0=S1/N;Ys0=S2/N; /計算外方位元素的初始值Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; while(j3) printf(-第%d次計算-n,j+1);a0=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);a1=-cos(t

10、)*sin(k)-sin(t)*sin(w)*cos(k);a2=-sin(t)*cos(w);b0=cos(w)*sin(k);b1=cos(w)*cos(k); /計算旋轉(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)*cos(k);c2=cos(t)*cos(w);for(i=0;iN;i+)x0i=-f*(a0*(Xi-Xs0)+b0*(Yi-Ys0)+c0*(Zi-Zs0)/(a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0); /計算像點坐標近似

11、值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;iN;i+)Ai*12+0=(a0*f+a2*xi)/Gi;Ai*12+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); /計算誤差方程的系數(shù)陣以及l(fā)Ai*12+4=-f*sin(k)

12、-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)*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,

13、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);mulAB(AT,A,ATA,6,2*N,2*N,6); /間接平差法計算外方位元素,中間計算的矩陣可以根據(jù)需要/printf(output matrix: ATAn); /選擇性的輸出，這里將其屏蔽，為了在打印輸出結(jié)果的時候/printmatrix(ATA,6,6); /節(jié)約資源ATA_=inv(ATA,6);/printf(output mat

14、rix: ATA_n);/printmatrix(ATA_,6,6);mulAB(AT,l,ATl,6,2*N,2*N,1);/printf(outpit matrinx: ATln);/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次計算的外方位元素為：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=%.5lf,Ys=%.5lf