前方后方空间交会试验报告

1、中南大学本科生课程设计（实践）任务书、设计报告（摄影测量与遥感概论）题 目 空间后方-前方交会学生姓名指导教师邹峥蝶学 院地球科学与信息物理学院专业班级测绘0902班学生学号实验目的通过对数字影像空间后交前交的程序设计实验，要求我们进理解和掌握影像外方位元素的有关理论、原理和方法。利用计算机程 序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位 元素及其精度，然后通过求得的外方位元素求解未知点的地面摄影测量坐标,达到通过摄影测量量测地面地理数据的目的。二实验要求A用C、VB或者Matlab编写空间后方交会前方交会计算机程序。提交实验报告:程序框图,程序源代码、计算结果及体会。计算结果:

2、地面点坐标、外方位元素及精度。完成时间：2011年月17日。三.实验数据占号左片右片地面摄影测量坐标XyXyXYZGCP116. 01279. 963-73.9378. 7065083. 2055852. 099527. 925GCP288. 5681.134-5.25278.1845780. 025906. 365571. 549GCP313. 362-79. 37-79.122-7& 8795210. 8794258.446461.81GCP482. 24-80. 027-9. 887-80. 0895909. 2644314. 283455. 484151. 75880. 555-39.

3、 95378. 463214. 618-0. 231-7& 0060.036349. 88-0. 782-42. 201-1. 022486.14-1. 346-7. 706-2.112548. 035-79. 962-44.438-79. 736f=150000mm z xo=0 z yo=0!1!实验思路禾I用后方交会得出两张像片各自的夕卜方位元素1）获取已知数据:从摄影资料中插曲像片比例尺 平均航高、内方位元素以及控制点的地面摄影测量坐标及对应的像点坐标。2)确定未知数的初始值:在竖直摄影的情况下，胶原素的初始值lyz lyx lyY3)计算旋转矩阵Ro为0,线元素其中Zs=m*f+4A

4、 , Xs=4A ，丫s厶4）逐点计算像点坐标的近似值:利用共线方程。5）组成误差方程并法化。6）解求外方位元素。刀检查计算是否收敛。利用解求出的外方位元素进行前方交会1）用各自像片的角元素计算出左右像片的旋转矩阵R1和R2。2）根据左右像片的外方位元素计算摄影基线分量Bx,By , Bze3）逐点计算像点的空间辅助坐标。4)计算投影系数。5)计算未知点的地面摄影测量坐标。6)重复以上步骤完成所有点的地面坐标的计算。五.实验过程程摩流程框图不 满 足 限 差 则 重 复 计 算将整个过程作为一 个函数继续进行右 片的外方位元素求 解求解各外方位元素 精度Z利用己求得的角元素 计算2张像片各自的

5、 旋转矩阵利用已求得的线元素畑Yd纣、pO,、w01, kOl： Xs2, Ys2, Zs2, p02, w02, k02,计算基线分量:Bx=Xs2-Xsi：Bv=Ys2-Ysl；Bz=Zs2-Zsl；J 丿输入像片坐标，利用旋 转矩阵求解想空间辅 助坐标计算点投影系数：Nl=(Bx*Z2-Bz*X2)/(Xl*Z2-X2*Zl)；N2=(Bx*Z l-Bz*Xl )/(X 1 *Z2-X2 *Z1)；计算地面摄影测量坐标Xt=(Nl*Xl+XsI)+(N2*X2+Xs2)/2 :Yt=(N 1 *Y1 + Ys 1 )+(N2 *Y2+Ys2)/2 :结束程序程摩中的主要函数设计子函数（矩

6、阵求积multiply ,计算函数Resection ,矩阵转置 transpose ,矩阵求逆inMersel ,输出函数shuchu z左片的外方位 元素求解函数zuobiane右片的夕卜方位元素求解函数youbiane ） 程摩源代码#include stdio-h1#include ”math.hdouble Xsl,Xs2,Ysl,Ys2,Zsl,Zs2,p01,p02,w01,w023d)l,k02;/求矩阵a的转置矩阵b, a为m行、n列void transpose(double *a，double int m, mt n);矩阵a乘以矩阵b,结果存储在cP, a为mXn大小，b

7、为nXl大小 void multiply(double *a, double *b, double *c, int hl int n. int 1);求矩阵a的逆mt uiMersel (double *a, int n);输出m行、n列的矩阵avoid shuchu(double *a, mt m, int n);计算并输出左片的外方位元素void zuobian();计算并输出右片的外方位元素void youbian();void zuobianQFILE *fp = NULL;FILE *fpl =NULL;if(fp=fbpen(,T:iniage.txt,Hf,) = NULL) p

8、rmtf(nOpen file error!H); return;if(Q)l=fbpen(,T:grouiid.txt,；,r,t) = NULL)prmtf(nOpen file error!n); return;像点坐标和地面点坐标double unagecontrol4 2=0.0;double groundcontrol43=0.0;摄影比例尺分母double m = 9943;double f=0.15;long ij,k;fbi(i=O; i4; i+)for(j=0;j2;j+)fscanf(fp. M%lf &imagecontrolij); unagecontrolij /

9、= 1000.0;for(k=0; k3; k+)fscanf(fpl, %lf &groundcontrolik);fclose(fp);fclose(fpl);计算外方位元素初始值fbr( i=0;i4;i+)Xs l+=groundcontioli 0;Y s l+=groundcontioli 1;Zsl +=groundcontroli 2;Xsl/=4.0;Ysl/=4.0;Zsl/=4.0;Zsl+=m*f;double R33=0.0;double L3=0QLl=0OL2=00;double L8l=0.0,x=0.0,y=0.0;double A86=0.0,AT 68=0

10、.0IA 6 6=0.0.B68=0.0;double V6l=0.0;int n=0;计算旋转矩阵R0 0=cos(p0 l)*cos(kO 1 )-siii(pO l)*siii(wO l)*sm(kO 1);R0l=(-l)*cos(p01)*sin(k01)-sm(p01)*sin(w01)*cos(k01);R0 2=(-l) *sin(pO 1) *cos(wO 1);R1 0=cos(w0 l)*sin(k01);Rl 1 =cos(wO l)*cos(kO 1);Rl2=(-l)*sin(w01);R2 0=siii(p0 l)*cos(kO l)+cos(pO l)*siii

11、(wO 1 )*sin(kO 1);R2 1=(-1) *sin(pO 1 )*sin(kO l)+cos(pO 1 )*sin(wO 1 )*cos(kO 1);R2 2=cos(p01 )*cos(wO 1);fbr(i=Oj=O:j =0.00001|fabs(V40)=0.00001|fabs(V50)=0.00001);R00=cos(p01)*cos(k01)-sm(p01)*sin(w01)*siii(k01);R0l=(-l)*cos(p01)*sin(k01)-sm(p01)*sm(w01)*cos(k01); R02=(-l)*sm(p01)*cos(w01);R1 0=c

12、os(w01 )*siii(kO 1);Rll =cos(wO l)*cos(kO 1);Rl2=(-l)*sm(w01);R20=sin(p01)*cos(k01)+cos(p01)*sm(w01)*sm(k01);R2l=(-l)*sm(p01)*sm(k01)+cos(p01)*sin(w01)*cos(k01);R22=cos(p01)*cos(w0i);进行未知数的精度评定double AV8l,X8l,XTl8,XTXlljiio.D664iii6; multiply(&A00,&V00,&AV00,&6,l);for(i=0;i8;i+)Xi0=AVi0-Li0;transpos

13、e(&X0 0 & XT0 0 1);multiply(&XT00h&X00h&XTX00丄&1); mo=XTX00/2;for(i=0;i6;i-H-)for(j=0j6j-H-)Di|j=mo*ATAi|j;fbr(i=0;i6;i-H-)mii=sqrt(Dii);pmirff左片结果为：nn”)； pmitf(”旋转矩阵R为： shuchu(&R00,3,3);printf(“外方为元素为:n”); pmitf(”Xs 1 =%lfiT；Xs 1); pmitfC Ys 1 =%lfiT； Ys 1); pnntfC，Zsl=%lfn,Zsl); pnntfC，p01=%lfn,p0

14、1);wO 1 =%lfjiw01); pnntfC，k01=%lfn,k01);各外方位元素精度为：iT);4,mi5)；pmuff迭代次数为：dnnnn”,ii); fclose(fp);void youbian()FILE *fp = NULL;FILE *fpl =NULL;i(fp=fbpen(,T:Vniage2.txt，7，rH) NULL)pimtf(nOpen file error!H); return;i(fpl=fbpenCT:ground.txtHrH) = NULL)pimtf(nOpen file error!H); return;像点坐标和地面点坐标double

15、unagecontrol42=0.0;double groundcontiol43=0.0;摄影比例尺分母double m = 10177;double f=0.15;long ij,k;fbi(i=O; i4; i+)for(j=0:j2;j+)fscanf(fp, %lf ； &imagecontrolij); unagecontrolij /= 1000.0;for(k=0; k3; k+)fscanf(fpl, %lf &groundcontrolik); fclose(fp); fclose(fpl);计算外方位元素初始值fbr( i=0;i4;i+)Xs2+=groundcontr

16、oli 0;Y s2+=groundcontioli 1;Zs2+=groundcontroli 2;Xs2/=4.0;Ys2/=4.0;Zs2/=4.0;Zs2+=m*f;double R33=0.0;double L3=0QLl=0OL2=00;double L8l=0.0,x=0.0,y=0.0;double A86=0.0,AT 68=0.0,AIA 6 6=0.0.B68=0.0;double V6l=0.0;int n=0;doR00=cos(p02)*cos(k02)-sin(p02)*siii(w02)*sin(k02);R0l=(-l)*cos(p02)*sin(k02)-s

17、in(p02)*sin(w02)*cos(k02);R0 2=(-l )*sin(p02) *cos(w02);R1 0=cos(w02)*sin(k02);Rl I=cos(w02)*cos(k02);Rl2=(-l)*sin(w02);R20=siii(p02)*cos(k02)+cos(p02)*siii(w02)*sin(k02);R2l=(-l)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R2 2=cos(p02)*cos(w02);fbr(i=Oj=O:j =0.00001|fabs(V40)=0.00001|fabs(V50)=0.

18、00001);R00=cos(p02)*cos(k02)-sm(p02)*sin(w02)*siii(k02);R0l=(-l)*cos(p02)*sin(k02)-sm(p02)*sm(w02)*cos(k02);R02=(-l)*sm(p02)*cos(w02);R1 0=cos(w02)*sin(k02);Rll =cos(w02)*cos(k02);Rl2=(-l)*sm(w02);R20=sin(p02)*cos(k02)+cos(p02)*sm(w02)*siii(k02);R2l=(-l)*sm(p02)*sm(k02)+cos(p02)*sm(w02)*cos(k02);R22

19、=cos(p02)*cos(w02);进行未知数的精度评定double AV8l,X8l,XTl8,XTXlljno.D664iii6;multiply( &A00.&V00,&AV00,&6,l);for(i=0;i8;i+)Xi0=AVi0-Li0;transpose(&X0 0 & XT0 0 81);multiply(&XT00h&X00h&XTX00丄&1);mo=XTX00/2;for(i=0;i6;i-H-)for(j=0j6j-H-)Di|j=mo*ATAi|j;fbr(i=0;i6;i-H-)mii=sqrt(Dii);pnntf(-K转矩阵 R 为：iiH)； shuchu

20、(&R0033)； pmitf(”外方为元素为：n”)； pnntf(MXs2=%lfn,Xs2); pnntf(MYs2=%lfn,Ys2); pnntf(MZs2=%lfnH,Zs2); pnntfC，pO2=%lfn,pO2); pnntf(MwO2=%lfji,w02); pnntfC，kO2=%lfn,kO2);各外方位元素精度为：iT);4,mi5)；pmitf(”迭代次数为：dm”,n); fclose(fp);/求矩阵a的转置矩阵b, a为m行、n列void transpose(double *a, double *b. iiit m, mt n)intij;for(i=0;i

21、m;i+)for(j=0jnj+)bj*m+i = ai*n+j；矩阵a乘以矩阵b,结果存储在cP, a为mXn大小，b为nXl大小 void multiply(double *a, double *b, double *c, int m. int n, int 1) int i,j,k;double t;for(i=0;im;i+)for(j=Ojlj+)t=0;fbr(k=O;kn;k-H-)t += ai*n+k*bk*14j；/求矩阵a的逆mt uiMersel (double *a. uit n)mt *is, *js, i,j,k, 1, u, M; double d、p;is =

22、 new mtn;js = new mtn;fbr(k=O; k=n-l; k+)d=0.0;fbr(i=k; i=n-l; i+)for(j=k;jd) d=p: isk=i;jsk=j; if(d+lO=lO)delete is;delete js;pimtf(nEnoi; not iiiMerse!nn); return 0;for(j=O;j=n-l;j+)u=k*n-+j; M=isk*n+j;p=au; au=aM; aM=p;1 扎jsk!=k)fbr(i=O; i=n-l; i+)u=i*n+k;M=i*n+jsk;p=au;au=aM;aM=p；l=k*n+k;al=1.0/

23、al;for(j=0;j=n-l;j-H-)】f(j!=k)u=k*n+j; au=au*al;for(i=0: i=n-l; i+)if(i!=k)他(J=O；j=n-l；j+)if(j!=k)u=i*n4j;au=au-ai*n-rk*ak*n+j； fbr(i=O; i=0; k-)if (jsk!=k)for(j=0;j=n-l;j+)u=k*n+j; M=jsk*n4j； p=au; au=aM; aM=p;if(isk!=k)fbr(i=O; i=n-l; i+)u=i*n+k; M=i*n+isk; p=au; au=aM; aM=p;二+UZ5七6强)七目d(+?vro-LOJ

24、) (+=UIOL)尼) a 二 ur .U二 ur/Rqnoppqonqs P-OA c蛊發启耳U ,归UT壬耀、宀-;amp】 sen utRP s 二utRPA%7X$二 A$二 乂冬邑 fHeATfHexy-HMT%HTXkueosq.oARqnop N耳Aqxqzxx7N TN7Z7X7XG TNT AUX7A,咲M二 x RqnopHeTeeo二 eqcno&Tq 二。neq7eqeq,2q7Tq 二 Tqcnee&ee 二 eque-足二 eEqnop olnnqnoA oUCTqonz().SST P-OAall=cos(pO l)*cos(kO 1 )-sin(pO l)*s

25、in(wO 1 )*sin(kO 1); al2=(-l)*cos(p01)*siii(k01)-sin(p01)*sm(w01)*cos(k01); a 13=(-l)*siii(pO 1 )*cos(vO 1);a21=cos(p02)*cos(k02)-siii(p02)*siii(w02)*sm(k02); a22=(-l)*cos(p02)*siii(k02)-sin(p02)*siii(w02)*cos(k02); a23=(-l)*sin(p02)*cos(w02);b 1 l=cos(wO l)*siii(kO 1);b 12=cos(w0 l)*cos(kO 1);bl3=(

26、-l)*sin(w01);b21 =cos(w02)*siii(k02);b22=cos(w02)*cos(k02);b23=(-1 )*sin(wO2);c ll=sin(pO l)*cos(kO 1 )+cos(pO l)*sin(wO l)*sm(kO 1); C12=(-l)*sm(p01)*siii(k01)+cos(p01)*sm(w01)*cos(k01); c 13=cos(pO 1 )*cos(wO 1);c21=siii(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02); c22=(-l)*sin(p02)*sm(k02)+cos(p02)*

27、siii(w02)*cos(k02); c23=cos(p02)*cos(w02);计算各待定点的像空间辅助坐标Xl=all*xl/1000+al2*yl/1000-al3*f;Yl=bll*xl/1000+bl2*yl/1000-bl3*f; Zl=cll*xl/1000+cl2*yl/1000-cl3*f;X2=a21 *x2/1000+a22*v2/1000-a2 3 *f;Y2=b21 *x2/1000+b22 *v2/1000-b23 *f;Z2=c21 *x2/1000+c22*v2/l 000-c23 *f;计算摄影基线分量Bx=Xs2-Xsl;Bv=Ys2-Ysl;JJBz=Z

28、s2-Zsl;计算点投影系数Nl=(Bx*Z2-Bz*X2)/(Xl*Z2-X2*Zl);N2=(Bx*Z l-Bz*Xl )/(X 1 *Z2-X2 *Z 1);计算待定点的地面摄影测量坐标Xt=(Nl *X 1 +Xs 1 )+(N2 *X2+Xs2)/2;Yt=(Nl * Y1 +Ys 1 )+(N2 * Y2+Ys2)/2;Zt=(N l*Zl+Zsl )+(N2 *Z2+Zs2)/2-1500;pmHfCXt=%lfuYt=%lfnZt=%lfir；XEYtZt)；实验结果左片结果为:旋转矩阵R为:0995467 -0.095109 -0.0002150.0950630995047 -0.0290600.0029780.0289080.999578外方为元素为：Xsl=4999.770153?sl=4999.728862Zsl=2000.002341p01=0.000215w01=0.029064k01=0.095247各外方位元素精度为:ml=0.2444?夕m2=0.318316m3=0.107253n4=0.000155m5=0.000153=0.000058迭代次数为：3n6=0