前方后方空間交會(huì)實(shí)驗(yàn)報(bào)告

1、中南大學(xué)本科生課程設(shè)計(jì)(實(shí)踐)任務(wù)書、設(shè)計(jì)報(bào)告 (攝影測(cè)量與遙感概論)題 目空間后方-前方交會(huì) 學(xué)生姓名指導(dǎo)教師鄒崢嶸學(xué) 院地球科學(xué)與信息物理學(xué)院專業(yè)班級(jí)測(cè)繪0902班學(xué)生學(xué)號(hào)一、 實(shí)驗(yàn)?zāi)康耐ㄟ^對(duì)數(shù)字影像空間后交前交的程序設(shè)計(jì)實(shí)驗(yàn)，要求我們進(jìn)一步理解和掌握影像外方位元素的有關(guān)理論、原理和方法。利用計(jì)算機(jī)程序設(shè)計(jì)語言編寫攝影測(cè)量空間交會(huì)軟件進(jìn)行快速確定影像的外方位元素及其精度，然后通過求得的外方位元素求解未知點(diǎn)的地面攝影測(cè)量坐標(biāo)，達(dá)到通過攝影測(cè)量量測(cè)地面地理數(shù)據(jù)的目的。二、 實(shí)驗(yàn)要求 用C、VB或者M(jìn)atlab編寫空間后方交會(huì)-前方交會(huì)計(jì)算機(jī)程序。 提交實(shí)驗(yàn)報(bào)告：程序框圖，程序源代碼、計(jì)算結(jié)果及

2、體會(huì)。 計(jì)算結(jié)果：地面點(diǎn)坐標(biāo)、外方位元素及精度。 完成時(shí)間：2011年11月17日。三、 實(shí)驗(yàn)數(shù)據(jù)點(diǎn)號(hào) 左片 右片 地面攝影測(cè)量坐標(biāo) x y x y X Y Z GCP1 16.012 79.963 -73.93 78.706 5083.205 5852.099 527.925 GCP2 88.56 81.134 -5.252 78.184 5780.02 5906.365 571.549 GCP3 13.362 -79.37 -79.122 -78.879 5210.879 4258.446 461.81 GCP4 82.24 -80.027 -9.887 -80.089 5909.264

3、4314.283 455.484 1 51.758 80.555 -39.953 78.463 2 14.618 -0.231 -76.006 0.036 3 49.88 -0.782 -42.201 -1.022 4 86.14 -1.346 -7.706 -2.112 5 48.035 -79.962 -44.438 -79.736 f=150.000mm，x0=0，y0=0四、 實(shí)驗(yàn)思路 利用后方交會(huì)得出兩張像片各自的外方位元素1) 獲取已知數(shù)據(jù)：從攝影資料中插曲像片比例尺、平均航高、內(nèi)方位元素以及控制點(diǎn)的地面攝影測(cè)量坐標(biāo)及對(duì)應(yīng)的像點(diǎn)坐標(biāo)。2) 確定未知數(shù)的初始值：在豎直攝影的情況下，膠

4、原素的初始值為0，線元素其中Zs=m*f+，Xs=，Ys=。3) 計(jì)算旋轉(zhuǎn)矩陣R。4) 逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值：利用共線方程。5) 組成誤差方程并法化。6) 解求外方位元素。7) 檢查計(jì)算是否收斂。 利用解求出的外方位元素進(jìn)行前方交會(huì)1) 用各自像片的角元素計(jì)算出左右像片的旋轉(zhuǎn)矩陣R1和R2。2) 根據(jù)左右像片的外方位元素計(jì)算攝影基線分量Bx,By，Bz。3) 逐點(diǎn)計(jì)算像點(diǎn)的空間輔助坐標(biāo)。4) 計(jì)算投影系數(shù)。5) 計(jì)算未知點(diǎn)的地面攝影測(cè)量坐標(biāo)。6) 重復(fù)以上步驟完成所有點(diǎn)的地面坐標(biāo)的計(jì)算。五、 實(shí)驗(yàn)過程 程序流程框圖后方交會(huì)函數(shù)確定已知數(shù)據(jù)比例尺m確定各外方位元素初始值計(jì)算旋轉(zhuǎn)矩陣逐點(diǎn)計(jì)算像

5、點(diǎn)坐標(biāo)近似值不滿足限差則重復(fù)計(jì)算逐點(diǎn)計(jì)算誤差方程系數(shù)項(xiàng)，組成誤差系數(shù)矩陣A利用矩陣運(yùn)算求解外方位元素檢查是否滿足限差若滿足則輸出外方位元素將整個(gè)過程作為一個(gè)函數(shù)繼續(xù)進(jìn)行右片的外方位元素求解求解各外方位元素精度空間前方交會(huì)利用已求得的角元素計(jì)算2張像片各自的旋轉(zhuǎn)矩陣?yán)靡亚蟮玫木€元素Xs1，Ys1，Zs1，p0，,w01，k01；Xs2，Ys2，Zs2，p02，w02，k02，計(jì)算基線分量：Bx=Xs2-Xs1；By=Ys2-Ys1；Bz=Zs2-Zs1；輸入像片坐標(biāo)，利用旋轉(zhuǎn)矩陣求解想空間輔助坐標(biāo)計(jì)算點(diǎn)投影系數(shù)：N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1)；N2=(Bx*Z1-

6、Bz*X1)/(X1*Z2-X2*Z1)；計(jì)算地面攝影測(cè)量坐標(biāo)Xt=(N1*X1+Xs1)+(N2*X2+Xs2)/2；Yt=(N1*Y1+Ys1)+(N2*Y2+Ys2)/2；Zt=(N1*Z1+Zs1)+(N2*Z2+Zs2) /2 ；結(jié)束程序 程序中的主要函數(shù)設(shè)計(jì)子函數(shù)（矩陣求積multiply，計(jì)算函數(shù)Resection，矩陣轉(zhuǎn)置transpose，矩陣求逆inMerse1，輸出函數(shù)shuchu，左片的外方位元素求解函數(shù)zuobian。右片的外方位元素求解函數(shù)youbian。） 程序源代碼#include stdio.h#include math.hdouble Xs1,Xs2,Ys1

7、,Ys2,Zs1,Zs2,p01,p02,w01,w02,k01,k02;/求矩陣a的轉(zhuǎn)置矩陣b，a為m行、n列void transpose(double *a, double *b, int m, int n);/矩陣a乘以矩陣b，結(jié)果存儲(chǔ)在c中，a為mn大小，b為nl大小void multiply(double *a, double *b, double *c, int m, int n, int l);/求矩陣a的逆int inMerse1(double *a, int n);/輸出m行、n列的矩陣avoid shuchu(double *a, int m, int n);/計(jì)算并輸出左

8、片的外方位元素void zuobian();/計(jì)算并輸出右片的外方位元素void youbian();void zuobian() FILE *fp = NULL;FILE *fp1 = NULL;if(fp=fopen(F:image.txt,r) = NULL)printf(Open file error!);return;if(fp1=fopen(F:ground.txt,r) = NULL)printf(Open file error!);return;/像點(diǎn)坐標(biāo)和地面點(diǎn)坐標(biāo)double imagecontrol42=0.0;double groundcontrol43=0.0;/攝影

9、比例尺分母double m = 9943;double f=0.15;long i,j,k; for(i=0; i4; i+)for(j=0; j2; j+)fscanf(fp, %lf, &imagecontrolij);imagecontrolij /= 1000.0;for(k=0; k3; k+)fscanf(fp1, %lf, &groundcontrolik);fclose(fp);fclose(fp1); /計(jì)算外方位元素初始值for( i=0;i4;i+)Xs1+=groundcontroli0;Ys1+=groundcontroli1;Zs1+=groundcontroli2

10、;Xs1/=4.0;Ys1/=4.0; Zs1/=4.0;Zs1+=m*f;double R33=0.0;double L3=0.0,L1=0.0,L2=0.0;double L81=0.0,x=0.0,y=0.0;double A86=0.0,AT68=0.0,ATA66=0.0,B68=0.0;double V61=0.0;int n=0;do/計(jì)算旋轉(zhuǎn)矩陣R00=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R01=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R02=(-1)*sin(p0

11、1)*cos(w01);R10=cos(w01)*sin(k01);R11=cos(w01)*cos(k01);R12=(-1)*sin(w01);R20=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R21=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R22=cos(p01)*cos(w01);for(i=0,j=0;j=0.00001|fabs(V40)=0.00001|fabs(V50)=0.00001);R00=cos(p01)*cos(k01)-sin(p01)*sin(w01)*si

12、n(k01);R01=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R02=(-1)*sin(p01)*cos(w01);R10=cos(w01)*sin(k01);R11=cos(w01)*cos(k01);R12=(-1)*sin(w01);R20=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R21=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R22=cos(p01)*cos(w01);/進(jìn)行未知數(shù)的精度評(píng)定double AV81,X8

13、1,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i8;i+)Xi0=AVi0-Li0;transpose(&X00,&XT00,8,1);multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2; for(i=0;i6;i+)for(j=0;j6;j+)Dij=mo*ATAij; for(i=0;i6;i+)mii=sqrt(Dii);printf(左片結(jié)果為：nn);printf(旋轉(zhuǎn)矩陣R為：n);shuchu(&R00,3,3);printf(外方為元素為：n);printf(

14、Xs1=%lfn,Xs1);printf(Ys1=%lfn,Ys1);printf(Zs1=%lfn,Zs1);printf(p01=%lfn,p01);printf(w01=%lfn,w01);printf(k01=%lfn,k01);printf(各外方位元素精度為：n); printf(m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn,mi0,mi1,mi2,mi3,mi4,mi5);printf(迭代次數(shù)為：%dnnnn,n);fclose(fp);void youbian()FILE *fp = NULL;FILE *fp1 = NULL;if(f

15、p=fopen(F:image2.txt,r) = NULL)printf(Open file error!);return;if(fp1=fopen(F:ground.txt,r) = NULL)printf(Open file error!);return;/像點(diǎn)坐標(biāo)和地面點(diǎn)坐標(biāo)double imagecontrol42=0.0;double groundcontrol43=0.0;/攝影比例尺分母double m = 10177;double f=0.15;long i,j,k; for(i=0; i4; i+)for(j=0; j2; j+)fscanf(fp, %lf, &image

16、controlij);imagecontrolij /= 1000.0;for(k=0; k3; k+)fscanf(fp1, %lf, &groundcontrolik);fclose(fp);fclose(fp1); /計(jì)算外方位元素初始值for( i=0;i4;i+)Xs2+=groundcontroli0;Ys2+=groundcontroli1;Zs2+=groundcontroli2;Xs2/=4.0;Ys2/=4.0; Zs2/=4.0;Zs2+=m*f;double R33=0.0;double L3=0.0,L1=0.0,L2=0.0;double L81=0.0,x=0.0

17、,y=0.0;double A86=0.0,AT68=0.0,ATA66=0.0,B68=0.0;double V61=0.0;int n=0;doR00=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R01=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R02=(-1)*sin(p02)*cos(w02);R10=cos(w02)*sin(k02);R11=cos(w02)*cos(k02);R12=(-1)*sin(w02);R20=sin(p02)*cos(k02)+cos(p02)*sin

18、(w02)*sin(k02);R21=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R22=cos(p02)*cos(w02);for(i=0,j=0;j=0.00001|fabs(V40)=0.00001|fabs(V50)=0.00001);R00=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R01=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R02=(-1)*sin(p02)*cos(w02);R10=cos(w02)*sin(k02

19、);R11=cos(w02)*cos(k02);R12=(-1)*sin(w02);R20=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R21=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R22=cos(p02)*cos(w02);/進(jìn)行未知數(shù)的精度評(píng)定double AV81,X81,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i8;i+)Xi0=AVi0-Li0;transpose(&X00,&XT00,8,1)

20、;multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2; for(i=0;i6;i+)for(j=0;j6;j+)Dij=mo*ATAij; for(i=0;i6;i+)mii=sqrt(Dii);printf(右片結(jié)果為：nn);printf(旋轉(zhuǎn)矩陣R為：n);shuchu(&R00,3,3);printf(外方為元素為：n);printf(Xs2=%lfn,Xs2);printf(Ys2=%lfn,Ys2);printf(Zs2=%lfn,Zs2);printf(p02=%lfn,p02);printf(w02=%lfn,w02);printf(k02

21、=%lfn,k02);printf(各外方位元素精度為：n); printf(m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn,mi0,mi1,mi2,mi3,mi4,mi5);printf(迭代次數(shù)為：%dn,n);fclose(fp);/求矩陣a的轉(zhuǎn)置矩陣b，a為m行、n列void transpose(double *a, double *b, int m, int n)int i,j;for(i=0;im;i+)for(j=0;jn;j+)bj*m+i = ai*n+j;/矩陣a乘以矩陣b，結(jié)果存儲(chǔ)在c中，a為mn大小，b為nl大小void multi

22、ply(double *a, double *b, double *c, int m, int n, int l)int i,j,k;double t;for(i=0;im;i+)for(j=0;jl;j+)t=0;for(k=0;kn;k+)t += ai*n+k*bk*l+j;ci*l+j=t;/求矩陣a的逆int inMerse1(double *a, int n)int *is, *js, i, j, k, l, u, M;double d,p;is = new intn;js = new intn;for(k=0; k=n-1; k+)d=0.0;for(i=k; i=n-1; i+

23、)for(j=k; jd) d=p; isk=i; jsk=j;if(d+1.0=1.0)delete is;delete js;printf(Error, not inMerse!n);return 0;if(isk != k)for(j=0; j=n-1; j+)u=k*n+j; M=isk*n+j;p=au; au=aM; aM=p; if(jsk!=k)for(i=0; i=n-1; i+)u=i*n+k;M=i*n+jsk;p=au; au=aM; aM=p;l=k*n+k;al=1.0/al;for(j=0; j=n-1; j+)if(j!=k)u=k*n+j; au=au*al;

24、for(i=0; i=n-1; i+)if(i!=k)for(j=0; j=n-1; 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; M=jsk*n+j;p=au; au=aM; aM=p; if(isk!=k)for(i=0; i=n-1; i+) u=i*n+k; M=i*n+isk;p=au; au=aM; aM=p;delete is;delete js;return 1;/輸出m行、n列的矩陣avoid shuchu(double *a,

25、int m, int n)int i,j;for(i=0;im;i+)for(j=0;jn;j+)printf(%9lf ,ai*n+j);printf(n);void main()zuobian();youbian();double a11,a12,a13,a21,a22,a23,b11,b12,b13,b21,b22,b23,c11,c12,c13,c21,c22,c23;double x1,y1,x2,y2,X1,Y1,Z1,X2,Y2,Z2,N1,N2,Xt,Yt,Zt,Bx,By,Bz;double f=0.15;scanf(x1=%lf,y1=%lf,x2=%lf,y2=%lf,&

26、x1,&y1,&x2,&y2); /scanf(y1=%lf,&y1);/printf(x2=); /scanf(x2=%lf,&x2);/旋轉(zhuǎn)矩陣a11=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01); a12=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);a13=(-1)*sin(p01)*cos(w01);a21=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02); a22=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*

27、cos(k02);a23=(-1)*sin(p02)*cos(w02);b11=cos(w01)*sin(k01);b12=cos(w01)*cos(k01);b13=(-1)*sin(w01);b21=cos(w02)*sin(k02);b22=cos(w02)*cos(k02);b23=(-1)*sin(w02);c11=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);c12=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);c13=cos(p01)*cos(w01);c21=sin(p02)*co

28、s(k02)+cos(p02)*sin(w02)*sin(k02);c22=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);c23=cos(p02)*cos(w02);/計(jì)算各待定點(diǎn)的像空間輔助坐標(biāo)X1=a11*x1/1000+a12*y1/1000-a13*f;Y1=b11*x1/1000+b12*y1/1000-b13*f;Z1=c11*x1/1000+c12*y1/1000-c13*f;X2=a21*x2/1000+a22*y2/1000-a23*f;Y2=b21*x2/1000+b22*y2/1000-b23*f;Z2=c21*x2/1000+c22*y2/1000-c23*f;/計(jì)算攝影基線分量Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;/計(jì)算點(diǎn)投影系數(shù)N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1);N2=(Bx*Z1-Bz*X1)/(X1*Z2-X2*Z1);/計(jì)算待定點(diǎn)的地面攝影測(cè)量坐標(biāo)Xt=(N1*X1+Xs1)+(N2*X2+Xs2)/2;Yt