空间后方交会程序————————————————————————————————作者：————————————————————————————————日期：ﻩ一． 实验目的:掌握摄影测量空间后方交会的原理,利用计算机编程语言实现空间后方交会外方位元素的解算。

二． 仪器用具及已知数据文件：计算机wind ｏws xp 系统,编程软件(VI ＳＵA Ｌ C ＋+6．0)，地面控制点在摄影测量坐标系中的坐标及其像点坐标文件shu ｊu.txt 。

三． 实验内容:单张影像的空间后方交会:利用已知地面控制点数据及相应像点坐标根据共线方程反求影像的外方位元素。

数学模型：共线条件方程式:)(3)(3)(3)(1)(1)(1Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f x -+-+--+-+--=)(3)(3)(3)(2)(2)(2Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f y -+-+--+-+--= 求解过程:（1)获取已知数据。

从航摄资料中查取平均航高与摄影机主距;获取控制点的地面测量坐标并转换为地面摄影测量坐标。

（2)量测控制点的像点坐标并做系统改正。

(３)确定未知数的初始值。

在竖直摄影且地面控制点大致分布均匀的情况下,按如下方法确定初始值,即:n X X S ∑=0,n Y Y S ∑=0,nZ mf Z S ∑=0φ =ω=κ=0式中;ｍ为摄影比例尺分母;ｎ为控制点个数。

（４)用三个角元素的初始值,计算个方向余弦，组成旋转矩阵R 。

(5)逐点计算像点坐标的近似值。

利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x ）、（y ）。

（6)逐点计算误差方程式的系数和常数项,组成误差方程式。

(7）计算法方程的系数矩阵A A T 和常数项l A T ,组成法方程式。

（8）解法方程，求得外方位元素的改正数dXs ,S dY ，s dZ ,d φ,ｄω，d κ。

（９）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。

(10)将求得的外方位元素改正数与规定的限差比较，若小于限差则迭代结束。

否则用新的近似值重复（4)~(９)，直到满足要求为止。

四．实验结果：程序的源代码如下所示:#iｎcｌude<stdiｏ.h>＃include<sｔdliｂ.h>#iｎclude＜ｍaｌloc．h>#inｃlude<ｍaｔｈ.h＞#ｉncｌude<coｎiｏ.h＞#ｄefｉｎe N 4voｉd turｎ(doｕble *Ａ,dｏｕble Ａ2[],int ｍ,int n）//计算矩阵的转置｛int i，ｊ;ﻩfｏｒ（i=0;i<m;ｉ+＋）for(j=0;j<n;j+＋)ﻩA2[ｊ*m＋i］=A[i＊ｎ＋j];｝vｏid mulAB(doｕble *Ａ,ｄoublｅ*Ｂ,ｄouｂle ＊C,inｔａm,int ａｎ,ｉnｔｂm，int bn) ／/计算两矩阵相乘{int i，j,l,u;iｆ(an！=bm）{ﻩpｒintf("erroｒ!ｃanｎot do thｅｍultiplicａｔiｏn．＼n");ﻩﻩreturn;ﻩ}ｆor(i=0;i<am;i＋＋)ﻩfor(j＝0;j<bn;j＋+）ﻩ{ﻩﻩu=ｉ＊bn+ｊ;C[ｕ］=0．０;ﻩfor(l=0；l<ａn；l+＋)ﻩﻩﻩC[u]+=A[i*ａn+l］*B［l＊bn+j];ﻩ}ﻩrｅturｎ;}douｂle *ｉnv(doｕblｅ*ａ,ｉｎt n) //计算矩阵的逆,本程序的难点｛//采用高斯-约旦－全选主元法iｎt *ｉs,*jｓ,ｉ,ｊ,ｋ,l，u,ｖ;douｂｌeｄ，p;ｉs＝(int*)mａlloc(n*sizｅof（ｉnt));js=(int＊)malloc（n*sｉzｅｏf(int））;ｆoｒ(k=０; ｋ<=ｎ-1; k++){ﻩd=0.０；fｏr(i＝ｋ;ｉ<ｎ;ｉ+＋)fｏｒ(j=ｋ;j<n;j++）{ l＝i*n+j;ﻩp=fabｓ(a［l]);ｉf (p>d)ﻩ{ﻩﻩﻩﻩd=p;ﻩﻩis［ｋ]=i;ﻩﻩjs[k]＝j；ﻩ}}ｉf (d+１.0=＝1.0）{ｆrｅe(is）;fｒｅe(js);printf("error not ｉnｖ\n");return NULL;｝if (ｉs[k]!＝k）for （ｊ＝0；ｊ<ｎ;j+＋){u=k*n+j;ﻩv=is［k]*n+j;p=ａ[u］;ﻩａ[u］＝a[ｖ];ﻩﻩa[v］=p;｝if(js[k]！=k)for （i=0;i<ｎ；i++)｛ｕ=i*n+ｋ;ﻩv=ｉ*n+js[k］;p＝a[u];ﻩa[ｕ]=a[v];ﻩa[v]=p；}l=k*n+ｋ;a[l]=1.0/a［ｌ］;foｒ(j=0；j<n;j++)iｆ(j!=k）｛ﻩu=ｋ＊n+j;ﻩﻩａ［ｕ］=a[u]*a[l]；ﻩ｝foｒ(i=0;i<n;i++)iｆ(i!=k)ｆor (ｊ=0;j＜n;ｊ+＋)ｉf (ｊ!=ｋ){ｕ＝i*n＋j;ａ［u]=a［u]-a[i*n+k］*a[k*n+ｊ];｝for(i＝0;ｉ<n;i＋＋)ｉf (i！=k){ u=ｉ＊n+k；ﻩa［u]=-ａ[u]＊a[l］;ﻩﻩ}｝ﻩfor (k=n－1;k>=0;k--){ if（jｓ[ｋ］!=k）ｆｏr (j=0;j<＝ｎ－1；j+＋){ u=k＊ｎ+j；ﻩﻩv=ｊｓ[ｋ]*n＋j；p=a［u];ﻩａ[u]=ａ［v］;ﻩﻩﻩa［ｖ]=p;}if(is[ｋ]!=ｋ)for (i=0;i＜n;ｉ++){u=i*n+k;ﻩv＝i*n+is[ｋ];p＝ａ［u];ﻩa［u］＝a［ｖ];ﻩﻩa[v］＝ｐ；｝｝ｆrｅe(ｉs）;ﻩfree(js);rｅtｕrn ａ;}voｉd ｐrintmａtrix(dｏuble M[］,ｉnt ｍ,int n) ／/矩阵的输出｛int i,j;for(i=０;i＜ｍ；ｉ+＋){for(j＝0;j＜n;j＋+)printf("%.５f＂，M[i*ｎ+j]）；ﻩpｒinｔf(＂\ｎ"）;}pｒiｎtf("\n"）;｝mａin() /／主函数,空间后方交会的计算{ﻩFＩLE*ｆp；//定义一个文件指针ｆpﻩint m，i，j=０；ﻩdｏubｌe f,t,w,k,Ｓ１=0．0,S2=0.0，S3=０.0,x[Ｎ],y[Ｎ］,x0[N］，y0［Ｎ],X［N］,Y[Ｎ],Z ［N］，Ｘｓ０，Ys0,Ｚs0;ﻩｄｏuｂle a[3],ｂ[3]，c[3],A［2*N＊6],AT[6*2＊Ｎ],ＡTA[6*6]，*ATＡ＿=NULL,l[２*N],ＡTl[6],V［6];if(（fp=fopen（"e：\\ｓｈuju.txt","r"))=＝NUＬＬ）/／使fp指向被打开的shujｕ.tｘt文件ﻩ{ //并判断文件是否打开成功ﻩｐriｎtｆ（＂\nerroｒon ｏpｅn sｈｕju．ｔｘt\n"）;geｔcｈ（）；eｘｉt(１);｝fｏｒ（i=0;i<Ｎ;i＋+）ﻩﻩfscanf(fp,"%lｆ％lf％lf%lf%lf",＆ｘ[i］,＆ｙ［i],&X［ｉ],＆Y[ｉ],&Ｚ[i])；//将文件中的数据赋给主函数定义的变量ﻩpｒｉnｔf("原始数据:＼n");pｒｉnｔf("ｘ\ｔ\ty\t\t\X\ｔ＼tY＼t＼tZ\ｔ\ｔ\n"); /／输出文件中的原始数据ｆor（i＝0；i<N;i++)ﻩﻩｐriｎtｆ("%lf %ｌf %lf %lf ％lf\n"，x[ｉ]，y[i],X[i],Y[i]，Z［i]);pｒintｆ(＂＼n请输入摄影机主距和摄影比例尺分母；f，m：");ｓcanf("％lf,%lf",&f,&ｍ); //输入f，mf=f/10０0.0;fｏr(i=0;i＜N;i++）ﻩ{ﻩx［ｉ]／=100０．0;ﻩﻩy［i]/=１00０.0;S1＋=X[i]；S２＋＝Y［i];ﻩS３+＝Z[i];}Xs０＝S１/N；ﻩYs0=S2/Ｎ; //计算外方位元素的初始值ﻩZs０=m*f+Ｓ３／N;t=0.0;w=０.0;ｋ=０.0；whiｌｅ(j<3)｛pｒｉnｔf(＂--－----－-－-----－－-－-－－-－--------－第%d次计算-----－-－---－－-－－-----－-－－---－-＼n"，j+1);ﻩa[0]=ｃｏs(t）＊cｏs(k）－siｎ（ｔ)＊sin(w)＊sｉｎ（k);a[1]=-cｏs（t)*ｓin(k)-ｓin（t)*sin(w)*cos(ｋ）;ﻩa［2]=-sin（t)*cｏs(ｗ);ｂ［０］=cos(ｗ)＊ｓin(k）；ﻩb[1]=cｏｓ(ｗ)＊ｃos(k); /／计算旋转矩阵b[2］=-sin(w);c[0]＝ｓin(t)*ｃos(k）+cos(t)*ｓin(w）*ｓin（k)；c［1]=-sin(ｔ)*sin(k)+cos（t)*ｓｉｎ(w）*coｓ（ｋ）;c[２]=ｃos(t)*cos(w);fｏr(ｉ=０;i<N；i++)ﻩ{ﻩx0[i]＝－f*（a［０]*（X[i]-Xｓ0)＋b[０］*(Ｙ[i]-Yｓ0)＋ｃ[0]*(Z[ｉ]-Zs０）)／（a[2]*（X[i]-Xs0)+ｂ[2]*(Ｙ［i]-Ｙs0）+c［2]*(Z[i]-Zs0)）; /／计算像点坐标近似值y0[i]＝-f*(a[1］*(X[i]-Ｘs0）+b[１]*(Ｙ[i］-Ys０)+c［1]*(Z［i]-Zs0))/(ａ[２]*(X［i］-Xｓ0）＋ｂ[2]*(Ｙ[i]-Ｙs0）＋c［2］*(Z[i]-Ｚs0));ﻩG[i］=ａ［２]*（X[i]-Xs０)＋b[2］*(Y［ｉ]-Ys０)+ｃ［2］*(Z［ｉ]－Zｓ0)；ﻩ｝for（i=0;i<N;i＋+)｛ﻩA[i*1２+0]=（a[0]*f+a[2］＊x［ｉ])/G［i]；A［ｉ*12＋1]=（b[0]＊ｆ+b[2]＊x［ｉ］)/G［i］;A［i*12+2］＝(c［0]*f＋ｃ[2］*x[i］)/G[i]；ﻩﻩA［i*1２+３]=y[ｉ]*sin(w)-(x[i]＊(x[i]*cｏs(k)-ｙ[i]*sｉn(ｋ）)/f＋f*cos(k）)＊cos(w）; ／/计算误差方程的系数阵以及lﻩA[i*1２+4]=-f＊siｎ(k)-ｘ[i］*（x[i］*sｉn(k)+ｙ[i］*cos(k））/f；ﻩＡ［i*1２+5]=ｙ[i];ﻩA[i*12+6］=(a[1]*ｆ+a[2]*y[i])/G［ｉ]；ﻩA［ｉ＊12+7]=（ｂ[1]*ｆ+b［2]*y[i]）／G［ｉ];ﻩﻩA［i*1２+8］=(c[1］*ｆ＋ｃ[2]＊y[ｉ］）/G[i]；ﻩＡ[i＊１２+９]＝-x[ｉ］*sｉn(ｗ)-(y[i］＊(x[i］＊cos（k）-y[i］*ｓｉｎ(k))/f-ｆ*sin(k））*cos(w);ﻩA[i＊12+10]＝-ｆ＊coｓ（k)-y[ｉ]*(x[i]*siｎ(k)＋ｙ[i］＊cos(k))/f;ﻩA[i*12+11]=-ｘ[i];ﻩﻩｌ[i＊2+0]=x［ｉ］－ｘ0[i];ﻩﻩl［ｉ*2+１]=ｙ[ｉ]－ｙ0[i］;ﻩ}// printf("ｏutput ｍaｔrix: Ａ\n"）；／/ﻩｐrintmatrｉx(Ａ,2*N,6）;//ﻩｐｒintf（"output matrix: ｌ\ｎ")；／/prｉｎtmatrix(l，2*N，1);turn(Ａ,AT,２*Ｎ,6);／/ prinｔf（"ｏutｐut matrix: AＴ\n")；／/ prｉnｔｍａtrｉx（A T,6,2*Ｎ)；mulAB(ＡＴ，A,ATA,６,2＊N,2*N,6); ／/间接平差法计算外方位元素,中间计算的矩阵可以根据需要//ﻩpｒiｎｔf(＂output matrｉx：ATA\n")；//选择性的输出,这里将其屏蔽,为了在打印输出结果的时候／/ﻩprintmatrix(AＴA,6，6);//节约资源ﻩATA_=ｉnv(ＡＴA,6);//ﻩprｉntf（"ｏuｔput mａtｒiｘ:A TA＿\n＂);//ﻩpｒintmatｒｉx（AＴＡ_,6,６)；mｕlAB(AT,l,ATl,６,２*N，2*Ｎ，1)；//ﻩprｉntf("ｏutpit matrinｘ: ＡTl\n");// printmａtriｘ(ATl，６,１)；mulAB（ATＡ_,A Tl,Ｖ,6，6,６,１)；// pｒinｔf（"ouｔpuｔmatriｘ: V\n");／/ﻩpriｎtｍatrix(Ｖ,６,1）；ﻩXs0+=V［0];Yｓ0+=Ｖ[1］;ﻩZs０+=V[2];ﻩt＋=Ｖ[3];ﻩw+=V[4］;k＋=V[5］;prｉｎtf（"第%ｄ次计算的外方位元素为：\n",++j）；ｐｒintf("Ｘｓ＝%.5lｆ,Ys=%.5lｆ，Zs=%.５lf，ｔ＝%.５lf,w=％．5ｌf,k=%.５ｌf\ｎ",Xs０,Ys0,Zs0,t,ｗ,k);}prｉntｆ(＂\n外方位元素为:\n");ｐrintｆ（"Xｓ＝％.5lf,Yｓ=％.5ｌf,Ｚs＝％．5ｌf,t=%.5lf,ｗ=％.5lf,k＝%．5lｆ\n",Xｓ0,Yｓ0,Zs０，ｔ,ｗ,ｋ)；ｆclｏｓe(ｆｐ);｝程序运行的结果为：五.实验总结：通过这次的实验我学到了很多的东西,通过编程加深了对摄影测量空间后方交会相关知识的理解。