作业4相对定向元素计算

1、相对定向元素计算实验报告时间：2013年5月24日 # 目录、已知与TOC o 1-5 h z HYPERLINK l bookmark10二、计算原理2 HYPERLINK l bookmark40三、相对定向元素的解算过程5四、图8/、.9七实j10丿I、弋io111丿611 一、已知数据与作业要求在一个航带内相邻两张像片上分别量测了6对同名点的像平面坐标，数据如表1,其中摄影机主距/=24（mm）,要求计算两张像片的相对定向元素。表1,已知点像平面坐标左片右片占号八Vx(nun)y(niin)x(nun)y(niin)1-6.0911.9835.564-3.2027.0980.9247.

2、694-2.8334.5381.0685.098-2.8786.8581.2087.429-2.5785-10.050.514-9.152-5.6426-8.0891.293-7.441-3.981二、计算原理从两个摄站对同一地面摄取一个立体像对时，立体像对中任意物点的两条同名光线都相交于该物点，即存在同名光线对对相交的现彖。若保持两张像片之间的相对位置和姿态关系不变，将两张相片整体移动、旋转和改变基线的长度，同名光线对对相交的特性并不发生变化。解析法相对定向就是根据同名光线对对相交这一立体像对内在的几何关系,通过量测的像点坐标，用解析计算的方法求相对定向元素，建立与地面相似的立体模型，确定模

3、型的三维坐标。相对定向与相片的绝对位置无关，不需要地面控制点。相对定向的共面条件在图1中，和色也和爲色为一条同名光线，与摄影基线B位于同一核面内，即SR、Sm,和2三条直线共面。由空间解析几何知识可知，如果三条直线共面，则它们一一对应矢量的混合积为零，即少佝色爲2）=（2-1-1）三矢量在像空间辅助坐标系中的坐标分别为（2小By,Bz）,（X，X，Zj和（KZZ）,则共面条件方程表示为BxF7X、BzZ=0Z2(2-1-2) # # # #连续像对的相对定向连续像对的相对定向是以左像片为基准，求出右像片相对于左像片的五个定向元素心。在相对定向解析计算时，通常把摄影基线E改写为乩b称为B=mb投

4、影基线，这里(2-2-1)其中：加为摄影比例尺分母；bx,bY,bz为投影基线对应的分量。为了统一单位,把两个基线元素改为角度形式表示，“z图2.bx,bybz空间关系 # # # #由图2可知bY=bxtanbx“bz=cos“tanvbzu # #式中：从u为基线的偏角和倾角。将式（2二二）代入共面条件方程得bxbxAbxv1“uX、XZ=bxX】XZX*儿乙wXrx乙(2-2-3) # # # 式中含有五个相对定向因素，其中隐含在（兀上，7）中，该式是一个非线性函数。为了平差计算，将式（223）按多元函数泰勒级数展开，取小值一次项，得共面方程的线性公式为:Ff煖如羞如養如羞如知0dFMX

5、込-bxZlY2l+bxXXp+bxY2v将各偏导数代入式(2-2-4),舍去含有“和U的二次小项，只保留一次小项，同时等式两边同除以得bx(2-2-5)顾及点投影系数得:牡厂华2迢严二禺人-=切一W=如n2n2 # # # #代入式(2-2-5),等式两边同乘以一垒，并近似地取Y1=Y2,Z1=Z2，则式(2-2-5)可化简为: # # # #bxd/.L-bxdv-z、z、Y2N2d(p-Z2+-kZ2FNN=d3+XNdK-，=0bxzi #最后得:Q=bxdA导bxdu警N2d(PZ2+YNedco+XNdK(2-2-6) # #式(2-2-6)即为连续法相对定向的解析计算公式。同时有

6、下式成立Q=N=_N齐七(2-2-7) # # # 式中：MX为左片投影点在以左摄站为原点的像空间辅助坐标系中的坐标；为右摄站的像空间辅助坐标系中的坐标：如为两摄站的丫坐标之差。所以0的几何意义是模型上同名点的Y坐标之差，称为上下视差，如图3所示。由前方交会公式可知，若同名光线相交于模型点，则bY=NN2g=0o相对定向元素的解算过程在相对定向公式(2-2-6)中，有五个未知数ddu,d(p,dco,dK,而在本算例中有六个标准点位的像平面坐标，因此可以进行迭代求解。具体步骤如下由于存在多余观测，根据最小二乘平差原理，将上下视差0作为观测值，可以写出误差方程式，即ve=M/z-N=d(p-7)

7、+N2dco+X2N2dK-Q(3-1-1)用一般符号表示误差方程式为：v=axdt+bxdv+cxd(p+dxdco+exdK-l(3-1-2)式中NrG=XqN?)=Qyv=Vq用矩阵表示误差方程式为v=abcdedjLLdud(pdcodK在本题中有6对像点，因此可以列出6个误差方程式，即VqdelV.brsdrCcV3ci3b3c3d3e3V4Q4bqC4勺V5ci5b5c5d5e5?6.方6Q“6dud(pdcodK # # (3-1-3)写成一般形式为6x1=Ax55xl_L&相对应的法方程为atpax-atpl=o一般情况卞像点坐标为等权观测，权阵P是单位矩阵，法方程可化简为at

8、ax-atl=o法方程的解为X=(ATAy1ATL解X即为相对定向元素近似值的改正数。由于误差方程式是根据泰勒级数展开的一次项近似公式，因此定向元素要用迭代法解求，具体过程如下：原始数据的输入及像点坐标的预处理在本算例中，由于题目给出的是像平面坐标，因此不需要进行从框标坐标到像平面坐标的转换。在由像平面坐标转换为像空间坐标时，将其第三维坐标统一设为摄影物镜像方节点到像平面的垂直距离，即-=-/=-24(mm)；确定相对定位元素的初始值在本算例中，为简化计算，将其视为航线的第一个像对，因此左像片的角元素为零。同时，右像片的三个角元素的初始值取零，基线分量幼,乞同样取零，即“0二q,二他=列二心二

9、0,取标准点位第一点的左右视差，即=(X_七)；计算左片的方向余弦值，组成旋转矩阵尺，计算左片6个像点的像空间辅助坐标(禺兀引;计算右片的方向余弦值，组成旋转矩阵尼，计算基线分量by和其中鸟可通过(P3K的初始值或其迭代后的值算得，by,乞的计算方法如卞by=bxxjL9bz=bxxv；计算右片各像点的像空间辅助坐标(兀,込,),计算各像点的点投影系数绚,弘和上下视差0,其中Q=NH-N0.0269、0.0045、Q.Q124rad。4.最终计算成果为=0.09110.0177/u=-0.03140.0221?0=-0042300269/adft)=-0.03010.00457-K=0.097

10、500124说七、实验心得与体会本次试验让我感到收获很大。1在知识方面，我学会了如何求解相对定向中相关的元素，并在理解的基础上做到了熟练应用。在编程方面，这是本次实验中我收获最大的地方。在本实验的C+程序中，我定义了一个类：POINT,它有3各数据成员，以及若干成员函数，我将6个像点的12个像空间坐标以及其对应的像空间辅助坐标全部定义为该类的对彖，这样使数据结构和程序的结构十分明了。此外我还用到了C+中的模板函数，使定义的函数有更广泛的适用性。为了保证实验的正确性，我在编写了C+程序的基础上，写了一个MATLAB的程序作为对照，从而保证了结果准确性。虽然已经完成了本次实验，但还是存在一些需要改

11、进的地方，比如我对于C+程序中在函数间传递二维数组的方法掌握的还不是很熟练，程序总体的模块化也不是很好，看起来比较零散，易读性较差。我还是会继续努力，提高自己的能力，使自己称为一个更优秀的人，感谢老师给了我这次锻炼自己的机会。八、程序截图及代码1-C+程序截图第1次迭代改正数分别为，0090703-0014503-0.055356-0.026160.098183第2次迭代改正数分别为：-0.0051741-0.0160450.011862-0.00360040.001093第3次迭代（改正数分别为：00047759-000137280.0012993-0.00031752-0.0015611第

12、4次迭代改正数分别为：000081188000043311-43554e-005-2.7233e-005-0.00024275第5次迭代改正数分别为：1.7762e-0058.3368e-005-34398e-0051.5655e-0066.155e-00?第6次迭代（改正数分别为：-1.5051e-0051.8367e-006-3.8592e-0067.7004e-0075.1413e-006u.u.w.K最终结果为:0.0911一00314一00423-0.03010.0975图4,程序截图 2.C+程序代码#includeBiostream#includeiomanip#include#

13、include#includeusingnamespacestd,#defineN5boolGauss(doubleA5,doubleB5,intn),采用部分主元的高斯消去法求方阵A的逆矩templatevoidtranspose(Tlmatl5,Tlmat26,inta,intb)intij;for(i=0,ia,i+)for()=0,jb,j+)mat2ji=matlij,templatevoidmulti655(Tlmatl6,Tlmat25,Tlmat35,inta,intb,intc)intijkfor(i=0,ia,i+)For()=0,jvc,j+)mat3ij=0;for(k

14、=0,kb,k+)mat3ij+=matlik*mat2kj,templatevoidmulti556(Tlmatl5,Tlmat26,Tlmat36,inta,intb,intc)intijkfor(i=0,ia,i+)for(j=0,jvoidshowmat(T*mat,inta,intb)fa(inti=0,ia,i+)for(intj=0,jxl0=point(l.983/1000,-6.091/1000,-24.0/1000),xll=point(0.924/1000,7.098/1000,-240/1000),xl2=point(l.068/1000,4.538/1000,-240

15、/1000),xl3=point(l.208/1000,6.858/1000,-240/1000),xl4=point(-0514/1000,-10.05/1000,-24.0/1000),xl5=point(l.293/1000,-8.089/1000,-24.0/1000),xil0=point(-3.202/1000,-5.564/1000,-24.0/1000),xill=point(-l83/1000,7.694/1000,-240/1000),xr2=point(-2.878/1000,5098/1000,-24.0/1000),xr3=point(?578/1000,7429/1

16、000,-240/1000),xr4=point(-5.642/1000,-9.152/1000,-24.0/1000);xr5=point(-3.981/1000,-7.441/1000,-24.0/1000);b=point(xl0.getx0-xi0getx(),0,0);Rl00=cos(X02)*cos(X04)-sin(X02)*sin(X03)*sin(X04)Rl0l=-cos(X02)*sin(X04)-sin(X02)*sin(X03)*cos(X04),Rl02=-sin(X02)*cos(X0),Rll0=cos(X03)*sin(X04);Rlll=cos(X03)*

17、cos(X04);Rll2-sin(X03D;Rl20=sin(X02)*cos(X04)+cos(X02)*sin(X03)*sin(X04);R121=-sin(X02)*sin(X04)+cos(X02)*sin(X03)*cos(X04);Rl22=cos(X02)*cos(X03);初始化旋转矩阵R23R200=cos(X02)*cos(X04)-sin(X02)*sin(X03)*sin(X04)R20l=-cos(X02)*sin(X0)-sin(X02)*sin(X03)*cos(X04),R202=-sin(X02)*cos(X03),R2l0=cos(X03)*sin(X0

18、4);R2ll=cos(X03)*cos(X04);R2l2-sin(X03D;R220=sin(X02)*cos(X04)+cos(X02)*sin(X03)*sin(X04);R22l=-sin(X02)*sin(X04)+cos(X02)*sin(X03)*cos(X04);R222=cos(X02)*cos(X03);pointXL6,XR6，/XL,XR分别存放相空间辅助坐标系下的个像点坐标for(inti=0,i=0.00003|X10=000003|X20=0.00003|X30=000003|X40=0.00003|X50=000003)time+;R200=cos(X02)*

19、cos(X04)-sin(X02)*sin(X03)*sin(X04);R20l=-cos(X02)*sin(X04)-sin(X02)*sin(X03)*cos(X04);R202=-sin(X02)*cos(X03),R210=cos(X03)*sin(X04),R211=cos(X03)*cos(X04),R2l2=-sin(X03),R220=sin(X02)*cos(X04)+cos(X02)*sin(X03)*sin(X04);R221=-sin(X02)*sin(X04)+cos(X02)*sin(X03)*cos(X04),R222=cos(X02)*cos(X03),for(

20、inti=0,ibz=bx*X0l,b=point(bx,by,bz)doubleA65=0,l6l=0,Nl61,N261,/构造系数矩阵和常数项。fa(inti=0,ii+)Nlil=(bgetx0*XRi.getzO-bgetz0*XRigetx0)/(XLigetx0*XRi-getz0-XR1.getx0*XL1getzO),N2il=(bgetx0*XL1.getzO-bgetz0*XL1.getx0)/(XL1.getx0*XR1getz0-XRi.getx0*XLigetzO)；Ai0M).getx0,Ail=-XRi,getyO/XR1getzO*bgetxO；Ai2=-XR

21、i,getx0*XRi.getyO/XRi.getz0*N2il,Ai3=-QCRigetz0+XR1.getyO*XR1getyO/XR1getzO)*N2il,Ai4=XRi.getxO*N2il；li0=Nlil*XLi.gety0-N2il*XRigetyO-bgetyO.doubleAr56=0,ATA5=0JNVATA55=0,INVATAAr=0,tianspose(A,AT,6,5),multi655(ATAATA,565),Gauss(ATA,INVArA,5)coutprecision(5),multi556(INVATAAT,INVATAAT,5,5,6);foiint1

22、=O,15,1+)XiO=O;for(intk=0,k6,k+)XiO+=INVATA宾Tik*lk0,coutendl,coutB第time次迭代u,v,W,w,K改正数分别为:endl;for(inti=0,i5,i+)XO1+=X1O;coutsetvz(14)Xi0,1ii.coutendl,coutfixed,cout.precision(4),coutBu,v,W,w,K最终结果为：n;for(inti=0,i5,i+)coutsetvz(l4)X01ii.getcharO,retuin0,boolGauss(doubleA5,doubleB5,mtn)int1,j,k,doubl

23、emax,temp,doubletNN,临时矩阵将A矩阵存放在临时矩阵tnn中for(i=0,in;i+)for(j=0,jn;j+)/初始化B矩阵为单位阵fcr(1=0,1vn,计+)for(j=0,jn,j+)(double)l:0,fcr(1=0,1vn,计+)寻找主元max=tii;k=i；for(j=i+l,jfabs(max)max=tji,k=J,如果主元所在行不是第i行，进行行交换if(k!=1)foi(j=0,jn,j+)tenp=tij,tij=tkj,tkj=tenp,/B伴随交换tenp=Bij;Bij=Bkj,Bkj=tenp;判断主元是否为0,若是，则矩阵A不是满秩

24、矩阵，不存在逆矩阵if(tii=O)coutThereisnoinveisematrix!,retiUTifalse,消去A的第1列除去1行以外的各行元素ten)=tii,for(j=0,jn;j+)tij=tij/tenip,主对角线上的元素变为1Bij=Bij/temp,伴随计算for(j=0,j第n行/不是第i行ten=tji,for(k=0,k-sin(X0(3)*cos(X0(4),R1(2,l)=cos(X0(4)*sin(X0(5),R12)=cos(X0(4)*cos(X0(5),Rl(2,3)=-sin(X0(4);Rl(3,l)=sin(X0(3)*cos(X0(5)+co

25、s(X0(3)*sin(X0(4)*sin(X0(5);R1(3,2)=-sin(X0(3)*sin(X0(5)+cos(X0(3)*sin(X0(4)*cos(X0(5);R1(3,3)=cos(X0(3)*cos(X0(4),%各像点的初始值,XL,XR都为3*6矩阵%xl=-6.0917.0984.5386.858-10.050-8.089,1.9830.9241.0681.208-0.5141.293,-24-24-24-24-24-24J/1000,%刃=5.5647.6945.0987.429-9.152-7.441,-3.202-2.830-2.878-2.578-5.642-3

26、.981;-24-24-24-24-24-24J/1000,xl=1.9830.9241.0681.208-0.5141.293,-6.0917.0984.5386.858-10.050-8.089,-24-24-24-24-24-24,)/1000,xr=3.202-2.830-2.878-2.578-5.642-3.981,-5.5647.6945.0987.429-9.152-7.441,-24-24-24-24-24-24J/1000,%组成旋转矩阵R2,X03,wX04,KX05n=0,bx=xl（l）xr（l）”%计算左像片在像空间辅助坐标系卜的坐标Xl,YlZh得到的坐标放在矩阵

27、XL中，矩阵大小为3*6for1=1:6XL(：ti)=Rl*xlC,i)；endR2(l,l)=cos(X0(3)*cos(X0(5)-sin(X0(3)*sin(X0(4)*sin(X0(5);R2(l,2)=.cos(X0(3)*sin(X0(5)-sin(X0(3)*sin(X0(4)*cos(X0(5);R2(l,3)=-sin(X0(3)*cos(X0(4),R2(2tl)=cos(X0(4)*sin(X0(5);R2(22)=cos(X0(4)*cos(X0(5),R2(2,3)=-sin(X0(4),R2(3,l)=sin(XO(3)*cos(XO(5)+cos(X0(3)*s

28、in(X0(4)*sin(X0(5);R2(32)=-sin(X0(3)*sin(X0(5)+cos(X0(3)*sin(X0(4)*cos(X0(5);R2(33)=cos(X0(3)*cos(X0(4);%计算右像片在像空间辅助坐标系卜的坐标X2,Y2,Z2,得到的坐标放在矩阵XR中，矩阵大小为3*6for1=1:6endwhile(n=0|abs(X(l)=0.00003|abs(X(2)=0.00003|abs(X(3)=0.00003|abs(X(4)=000003|abs(X(5)=000003)R2(l,l)=cos(X0(3)*cos(X0(5)-sin(X0(3)*sin(X

29、0(4)*sin(X0(5);R2(l,2)=.cos(X0(3)*sin(X0(5)-sin(X0(3)*sin(X0(4)*cos(X0(5);R2(l,3)=-sin(X0(3)*cos(X0(4),R2(2tl)=cos(X0(4)*sin(X0(5);R2(22)=cos(X0(4)*cos(X0(5),R2(2,3)=-sin(X0(4),R2(3,l)=sin(X0(3)*cos(X0(5)+cos(X0(3)*sin(X0(4)*sin(X0(5);R2(32)=-sin(X0(3)*sin(X0(5)+cos(X0(3)*sin(X0(4)*cos(X0(5);R2(33)=cos(X0(3)*cos(X0(4);%计算右像片在