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

1、作業(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，YC，ZC），與其影像上對應(yīng)的三個像點的影像坐標(biāo)a（

2、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.314作業(yè)過程 41 獲取已知數(shù)據(jù) 相片比例尺1/m=1：

3、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ù)的近似值按共線方程式計算控制點像點坐標(biāo)的近似值（x）（y）。46 組成誤差方程： 逐點

4、計算誤差方程式的系數(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 DOUBLE5;int InputData(int &N

5、um, 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 (Data!=NULL) delete Data; r

6、eturn 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=order | r=order) system(dat

7、a.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.97 24934.98 2386.504 10.46

8、 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,&Num,&m);if (Num4) return 2;f

9、scanf(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 Datai0/=1000.0; Datai1/=1000.0;

10、/如果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ù)),com(第一個*矩陣列數(shù)),(out)JgJz(存

11、放結(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 *函數(shù)介紹：向結(jié)果.txt文件輸出數(shù)據(jù)*參數(shù)：Q協(xié)因數(shù)

12、陣，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é)果.txt,w) return 1;FILE *fp_inpu

13、t;if (!(fp_input=fopen(data.txt,r) return 1;fprintf(fp_out,* * * *n);fprintf(fp_out,n空間后方交會程序(CC+)n測繪一班n 學(xué)號：n姓名：陳聞亞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_out,攝影比例尺（0表示其值位置）：);fprintf(fp_ou

14、t,%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%10lft%10lft%10lft%10lft%10lfn, i+1,temp

15、0,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);fprintf(fp_out,tKappa=%10lfnn,Kappa);fpr

16、intf(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+i*6+j); fprintf(fp_out,n);fprintf(

17、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 &m, const double &f)double Xs=0,Ys

18、=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;double *ATL=new double6;double *Xg

19、=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(Phi)*sin(Omega)*sin(Kappa); R21=-si

20、n(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-Zs); double xxx,yyy; xxx=-f*X/Z; yy

21、y=-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)* sin(Kappa)+(yyy)*cos(Kappa); A2*

22、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=-(xxx); /求矩陣A的轉(zhuǎn)置矩陣AT for (i=0;i2*Nu

23、m;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*Num;i+) Vi-=Li; VTV+=Vi*Vi;double m

24、0=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;b=temp;/*函數(shù)名：InverseMatrix *函數(shù)介紹：求矩

25、陣的逆（高斯-約當(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;jrow;j+) *pt=*ptemp; ptemp+; pt+; i

26、nt 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+