作業(yè)4空間后方交會

1、攝影測量作業(yè)報告之空間后方交會陳聞亞 20080729作業(yè)報告空間后方交會專 業(yè)： 測繪工程 班 級： 2008級(1)班 姓 名： 陳聞亞 指導(dǎo)教師： 陳強(qiáng) 2010 年 4 月 16 日1 作業(yè)任務(wù) -32 作業(yè)思想 -33 作業(yè)條件及數(shù)據(jù)-34 作業(yè)過程 -35 源程序 - 46 計算結(jié)果 - 177心得體會與建議- 171 作業(yè)任務(wù)計算近似垂直攝影情況下后方交會解。即利用攝影測量空間后方交會的方法，獲取相片的6個外方位元素。限差為0.1。2作業(yè)思想利用攝影測量空間后方交會的方法求解。該方法的基本思想是利用至少三個一直地面控制點的坐標(biāo)a（xa，ya，za）、b（xb，yb，zb）c（xc

2、，yc，zc），與其影像上對應(yīng)的三個像點的影像坐標(biāo)a（xa，ya）、b（xb，yb）、c（xc，yc），根據(jù)共線方程，反求該相片的外方位元素xs、ys、zs、。3作業(yè)條件及數(shù)據(jù)已知攝影機(jī)主距f=153.24mm，四對點的像點坐標(biāo)與相應(yīng)的地面坐標(biāo)列入下表：表1點號像點坐標(biāo)地面坐標(biāo)x（mm）y(mm)x(m)y(m)z(m)1-86.15-68.9936589.4125273.322195.172-53.40-82.2137631.0831324.51.728.69314.7876.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31

3、4作業(yè)過程 41 獲取已知數(shù)據(jù) 相片比例尺1/m=1：10000，內(nèi)方位元素f=153.24mm，x0，y0；獲取控制點的地面測量坐標(biāo)xt、yt、zt。42 量測控制點的像點坐標(biāo)： 本次作業(yè)中為已知。見表1。43 確定未知數(shù)的初始值： 在近似垂直攝影情況下，膠原素的初始值為0，即0 = 0 = 0=0；線元素中，zs0=h=mf=1532.4m，xs0、ys0的取值可用四個控制點坐標(biāo)的平均值，即：xs0= =38437.00ys0= =89106.6244 計算旋轉(zhuǎn)矩陣r： 利用膠原素的近似值計算方向余弦值，組成r陣。45 逐點計算像點坐標(biāo)的近似值： 利用未知數(shù)的近似值按共線方程式計算控制點像

4、點坐標(biāo)的近似值（x）（y）。46 組成誤差方程： 逐點計算誤差方程式的系數(shù)和常數(shù)項。47 組成法方程式： 計算法方程的系數(shù)矩陣ata與常數(shù)項atl。48 求解外方位元素： 根據(jù)法方程，由式x=（ata）-1 atl解求外方位元素改正數(shù)，并與相應(yīng)的近似值求和，得到外方位元素新的近似值。49 求解外方位元素： 將求得的外方位元素的改正數(shù)與規(guī)定的限差（0.1）比較，小于限差則計算終止，否則用新的近似值重復(fù)第4.4至4.8步驟的計算，知道滿足要求為止。5 源程序#include #include #include const double precision=1e-5;typedef double d

5、ouble5;int inputdata(int &num, double *&data,double &m,double &f);int resection(const int &num,const double *&data,const double &m,const double &f);int inversematrix(double *matrix,const int &row);int main(int argc, char* argv)double *data=null;int num;double f(0),m(0);if(inputdata(num,data,m,f) if

6、(data!=null) delete data; return 1;if(resection(num,data,m,f) if (data!=null) delete data; return 1;if (data!=null) delete data;printf(解算完畢.n);do printf(計算結(jié)果保存于結(jié)果.txt文件中n 請選擇操作（輸入p打開結(jié)果數(shù)據(jù)，r打開原始數(shù)據(jù)，其它退出程序）：); fflush(stdin); /刷新輸入流 char order=getchar(); if (p=order | p=order) system(結(jié)果.txt); else if (r=

7、order | r=order) system(data.txt); else break; system(cls);while(1);system(pause);return 0;/*函數(shù)名：inputdata *函數(shù)介紹：從文件（data.txt）中讀取數(shù)據(jù),*文件格式如下： *點數(shù) m(未知寫作0)* 內(nèi)方位元素(f x0 y0)*編號 x y x y z*實例：4 0153.24 0 01 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100

8、.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31*參數(shù)：(in/out)num(點數(shù))，*(in/out)data(存放數(shù)據(jù)),m,f,x0,y0*返回值：int ，0成功，1文件打開失敗，2控制點個*數(shù)不足，3文件格式錯誤*/int inputdata(int &num, double *&data,double &m,double &f)double x0,y0;file *fp_input;if (!(fp_input=fopen(data.txt,r) return 1;fscanf(fp_input,%d%lf,&n

9、um,&m);if (num4) return 2;fscanf(fp_input,%lf%lf%lf,&f,&x0,&y0);f/=1000;if (m0 | f0) return 3;data=new doublenum;double *temp= new doublenum-1;double scale=0;int i;for (i=0;inum;i+) /讀取數(shù)據(jù)，忽略編號 if(fscanf(fp_input,%*d%lf%lf%lf%lf%lf, &datai0,&datai1,&datai2, &datai3,&datai4)!=5) return 3; /單位換算成m data

10、i0/=1000.0; datai1/=1000.0;/如果m未知則歸算其值if (0=m) for (i=0;inum-1;i+) tempi=(datai2-datai+12)/(datai0-datai+10)+ (datai3-datai+13)/(datai1-datai+11); scale+=tempi/2.0; m=scale/(num-1);fclose(fp_input);delete temp;return 0;/*函數(shù)名：matrixmul *函數(shù)介紹：求兩個矩陣的積,*參數(shù)：jz1(第一個矩陣),row(第一個矩陣行數(shù)),*jz2(第二個矩陣),row(第二個矩陣列數(shù)

11、),com(第一個*矩陣列數(shù)),(out)jgjz(存放結(jié)果矩陣)*返回值：void*/void matrixmul(double *jz1,const int &row,double *jz2, const int &line,const int &com,double *jgjz)for (int i=0;irow;i+) for (int j=0;jline;j+) double temp=0; for (int k=0;kcom;k+) temp+=*(jz1+i*com+k)*(*(jz2+k*line+j); *(jgjz+i*line+j)=temp; /*函數(shù)名：output

12、*函數(shù)介紹：向結(jié)果.txt文件輸出數(shù)據(jù)*參數(shù)：q協(xié)因數(shù)陣，m精度，m0單位權(quán)中誤差，6個外*方位元素，旋轉(zhuǎn)矩陣*返回值：int，0成功，1失敗*/int output(const double *&q,const double *&m,const double &m0, const double &xs,const double &ys,const double &zs, const double &phi,const double &omega, const double &kappa,const double *r)file *fp_out;if (!(fp_out=fopen(結(jié)果.tx

13、t,w) return 1;file *fp_input;if (!(fp_input=fopen(data.txt,r) return 1;fprintf(fp_out,* * * *n);fprintf(fp_out,n空間后方交會程序(cc+)n測繪一班n 學(xué)號：20080729n姓名：陳聞亞nn);fprintf(fp_out,* * * *n);fprintf(fp_out,已知數(shù)據(jù)：nn已知點數(shù)：);int num;double temp,x,y;fscanf(fp_input,%d%lf,&num,&temp);fprintf(fp_out,%dn,num);fprintf(fp

14、_out,攝影比例尺（0表示其值位置）：);fprintf(fp_out,%10.0lfn,temp);fprintf(fp_out,內(nèi)方位元素(f x0 y0)：);fscanf(fp_input,%lf%lf%lf,&temp,&x,&y);fprintf(fp_out,%10lft%10lft%10lfn,temp,x,y);for (int i=0;inum;i+) double temp5; fscanf(fp_input,%*d%lf%lf%lf%lf%lf, &temp0,&temp1,&temp2,&temp3,&temp4); fprintf(fp_out,%3dt%10lf

15、t%10lft%10lft%10lft%10lfn, i+1,temp0,temp1,temp2,temp3,temp4);fclose(fp_input);fprintf(fp_out,* * * *n);fprintf(fp_out,計算結(jié)果如下：nn外方位元素：n);fprintf(fp_out,txs=%10lfn,xs);fprintf(fp_out,tys=%10lfn,ys);fprintf(fp_out,tzs=%10lfn,zs);fprintf(fp_out,tphi=%10lfn,phi);fprintf(fp_out,tomega=%10lfn,omega);fprin

16、tf(fp_out,tkappa=%10lfnn,kappa);fprintf(fp_out,旋轉(zhuǎn)矩陣：n);for (i=0;i3;i+) fprintf(fp_out,t); for (int j=0;j3;j+) fprintf(fp_out,%10lft,*(r+i*3+j); fprintf(fp_out,n);fprintf(fp_out,n單位權(quán)中誤差：%10lfnn,m0);fprintf(fp_out,協(xié)因數(shù)陣：n);for (i=0;i6;i+) fprintf(fp_out,t); for (int j=0;j6;j+) fprintf(fp_out,%20lft,*(q

17、+i*6+j); fprintf(fp_out,n);fprintf(fp_out,n外方位元素精度：);for (i=0;i6;i+) fprintf(fp_out,%10lft,mi);fprintf(fp_out,n);fprintf(fp_out,* * * *n);fclose(fp_out);return 0;/*函數(shù)名：resection *函數(shù)介紹：計算*參數(shù)：num(點數(shù))，data(數(shù)據(jù))，m,f(焦距),x0,y0*返回值：int，0成功，其它失敗*/int resection(const int &num,const double *&data,const double

18、 &m, const double &f)double xs=0,ys=0,zs=0;int i,j;/設(shè)置初始值for (i=0;inum;i+) xs+=datai2; ys+=datai3;xs/=num;ys/=num;zs=m*f;double phi(0),omega(0),kappa(0);double r33=0.0;double *l=new double2*num;typedef double double66;double6 *a=new double62*num;double *at=new double2*num*6;double *ata=new double6*6

19、;double *atl=new double6;double *xg=new double6;/迭代計算do /旋轉(zhuǎn)矩陣 r00=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa); r01=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa); r02=-sin(phi)*cos(omega); r10=cos(omega)*sin(kappa); r11=cos(omega)*cos(kappa); r12=-sin(omega); r20=sin(phi)*cos(kappa)+cos(

20、phi)*sin(omega)*sin(kappa); r21=-sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa); r22=cos(phi)*cos(omega); for (i=0;inum;i+) double x=r00*(datai2-xs)+r10*(datai3-ys)+ r20*(datai4-zs); double y=r01*(datai2-xs)+r11*(datai3-ys)+ r21*(datai4-zs); double z=r02*(datai2-xs)+r12*(datai3-ys)+ r22*(datai4-

21、zs); double xxx,yyy; xxx=-f*x/z; yyy=-f*y/z; /常數(shù)項 l2*i=datai0-(-f*x/z); l2*i+1=datai1-(-f*y/z); a2*i0=(r00*f+r02*(xxx)/z; a2*i1=(r10*f+r12*(xxx)/z; a2*i2=(r20*f+r22*(xxx)/z; a2*i3=(yyy)*sin(omega)-(xxx)/f)* (xxx)*cos(kappa)-(yyy)*sin(kappa)+ f*cos(kappa)*cos(omega); a2*i4=-f*sin(kappa)-(xxx)/f)*(xxx

22、)* sin(kappa)+(yyy)*cos(kappa); a2*i5=(yyy); a2*i+10=(r01*f+r02*(yyy)/z; a2*i+11=(r11*f+r12*(yyy)/z; a2*i+12=(r21*f+r22*(yyy)/z; a2*i+13=-(xxx)*sin(omega)-(yyy)/f)* (xxx)*cos(kappa)-(yyy)*sin(kappa)- f*sin(kappa)*cos(omega); a2*i+14=-f*cos(kappa)-(yyy)/f)*(xxx)* sin(kappa)+(yyy)*cos(kappa); a2*i+15=

23、-(xxx); /求矩陣a的轉(zhuǎn)置矩陣at for (i=0;i2*num;i+) for (j=0;j=precision |fabs(xg1)=precision | fabs(xg2)=precision |fabs(xg3)=precision | fabs(xg4)=precision | (xg5)=precision);/注：協(xié)因數(shù)陣，旋轉(zhuǎn)矩陣等計算本應(yīng)該使用最后外方位元素值，/由于變換很小忽略double *q=ata;double *v=new double2*num;matrixmul(&a00,2*num,xg,1,6,v);double vtv=0;for(i=0;i2*

24、num;i+) vi-=li; vtv+=vi*vi;double m0=sqrt(vtv/(2*num-6);double *mm=new double6;for (i=0;i6;i+) mmi=sqrt(*(q+i*6+i)*m0;output(q,mm,m0,xs,ys,zs,phi,omega,kappa,&r00);delete l;delete a;delete at;delete ata;delete atl;delete xg;delete mm;delete v;return 0;void swap(double &a,double &b)double temp=a;a=b;

25、b=temp;/*函數(shù)名：inversematrix *函數(shù)介紹：求矩陣的逆（高斯-約當(dāng)法） *輸入?yún)?shù)：(in/out)matrix（矩陣首地址），*（in）row（矩陣階數(shù)）*輸出參數(shù)：matrix（原矩陣的逆矩陣）*返回值：int ，0成功，1失敗*調(diào)用函數(shù)：swap(double&,double&)*/int inversematrix(double *matrix,const int &row)double *m=new doublerow*row;double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;irow;i+) for (j=0;

26、jrow;j+) *pt=*ptemp; ptemp+; pt+; int k;int *is=new introw,*js=new introw;for (k=0;krow;k+) double max=0; /全選主元 /尋找最大元素 for (i=k;irow;i+) for (j=k;jmax) max=*(m+i*row+j); isk=i; jsk=j; if (0 = max) return 1; /行交換 if (isk!=k) for (i=0;irow;i+) swap(*(m+k*row+i),*(m+isk*row+i); /列交換 if (jsk!=k) for (i=0;irow;i+) swap(*(m+i*row+k),*(m+i*row+jsk); *(m+k*row+k)=1/(*(m+k*row+k); for (j=0;jrow;j+) if (j!=k) *(m+k*row+j)*=*(m+k*row+k); for (i=0;irow;i+) if (i!=k) for (j=0;jrow;j+) if(j!=k) *(m+i*row+j)-=*(m+i*row+k)*(m+