实验一1实验任务理解摄影测量中核心模型-共线方程作用,掌握航空影像中重要的点线面的透视关系以及物方与像方之间的解析关系,单幅影像上像点坐标与相应地面点坐标之间的关系。通过编程实现外方位元素的求解,提升编程能力。2理论模型与方法单张像片的空间后方交会的基本思想:以单幅影像为基础,从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应的像坐标量测值处发,根据共线条件方程,解求该影像在航空摄影时刻的元素SX,SY,SZ,,,。(1)共线方程)()()()()()()()()()()()(33322203331110SASASASASASASASASASASASAZZcYYbXXaZZcYYbXXafyyZZcYYbXXaZZcYYbXXafxx(2)旋转矩阵R123123123coscossinsinsincossinsinsincossincoscossincoscossinsincoscossinsinsinsincossincoscoscosaaabbbccc(3)经过线性化,得到x,y的误差方程式yxaaaZaYaXaxxaaaZaYaXasssysssx)(v)(v262524232221161514131211矩阵形式如下:LAXV系数方程262524232221161514131211aaaaaaaaaaaaA其中:xayxfyfafyxfyxayayxfxfafyxfxyaZycfcaZybfbaZyafaaZxcfcaZxbfbaZxafaa262524161514322332223221311331123111)cossin(/coscos]sin)sincos(/[sin]cossin[/sincos}cos]sincos[/{sin/][/][/][/][/][/][TyyxxL)(,)()()()()()()()()()(333222111SSSSSSSSSZZcYYbXXaZZZcYYbXXaYZZcYYbXXaX近似值计算公式如下:ZYfyyZXfxx//00(4)由最小二乘间接平差原理可得:SSSZYXXLAAAXTT1)(210210210210210210SSSSSSSSSSSSZZZZYYYYXXXX3程序设计本地方仅列出核心代码:%确定初值x0=0;y0=0;f=153.24;m=sqrt(((x(1)-x(2))^2+(y(1)-y(2))^2))/(sqrt(((X(1)-X(2))^2+(Y(1)-Y(2))^2)));Zs=f/m;%3个线性元素Xs=mean(X);Ys=mean(Y);aa=0;%3个外方位角元素初值ww=0;kk=0;cx=zeros(6,1);p=0.1/206264.806247096363;%将0.1秒限差化为弧度aa1=1;ww1=1;kk1=1;k=0;whileabs(aa-aa1)p||abs(ww-ww1)p||abs(kk-kk1)paa1=aa;%赋值ww1=ww;kk1=kk;%计算旋转矩阵a1=cos(aa)*cos(kk)-sin(aa)*sin(ww)*sin(kk);a2=-cos(aa)*sin(kk)-sin(aa)*sin(ww)*cos(kk);a3=-sin(aa)*cos(ww);b1=cos(ww)*sin(kk);b2=cos(ww)*cos(kk);b3=-sin(ww);c1=sin(aa)*cos(kk)+cos(aa)*sin(ww)*sin(kk);c2=-sin(aa)*sin(kk)+cos(aa)*sin(ww)*cos(kk);c3=cos(aa)*cos(ww);R=[a1a2a3;b1b2b3;c1c2c3];%计算误差方程系数fori=1:d1%计算近似值XX=a1*(X(i)-Xs)+b1*(Y(i)-Ys)+c1*(Z(i)-Zs);YY=a2*(X(i)-Xs)+b2*(Y(i)-Ys)+c2*(Z(i)-Zs);ZZ=a3*(X(i)-Xs)+b3*(Y(i)-Ys)+c3*(Z(i)-Zs);a11=(1/ZZ)*(a1*f+a3*(x(i)-x0));a12=(1/ZZ)*(b1*f+b3*(x(i)-x0));a13=(1/ZZ)*(c1*f+c3*(x(i)-x0));a14=(y(i)-y0)*sin(ww)-(((x(i)-x0)/f)*((x(i)-x0)*cos(kk)-(y(i)-y0)*sin(kk))+f*cos(kk))*cos(ww);a15=-f*sin(kk)-((x(i)-x0)/f)*((x(i)-x0)*sin(kk)+(y(i)-y0)*cos(kk));a16=(y(i)-y0);a21=(1/ZZ)*(a2*f+a3*(y(i)-y0));a22=(1/ZZ)*(b2*f+b3*(y(i)-y0));a23=(1/ZZ)*(c2*f+c3*(y(i)-y0));a24=-(x(i)-x0)*sin(ww)-(((y(i)-y0)/f)*((x(i)-x0)*cos(kk)-(y(i)-y0)*sin(kk))-f*sin(kk))*cos(ww);a25=-f*cos(kk)-((y(i)-y0)/f)*((x(i)-x0)*sin(kk)+(y(i)-y0)*cos(kk));a26=-(x(i)-x0);A(2*i-1,:)=[a11a12a13a14a15a16];A(2*i,:)=[a21a22a23a24a25a26];l(2*i-1,:)=x(i)-(x0-f*XX/ZZ);l(2*i,:)=y(i)-(y0-f*YY/ZZ);endcx1=inv(A'*A)*A'*l;%最小二乘平差Xs=Xs+cx1(1);Ys=Ys+cx1(2);Zs=Zs+cx1(3);aa=aa+cx1(4);ww=ww+cx1(5);kk=kk+cx1(6);k=k+1;%统计循环次数endV==A*cx1-l;sigma=sqrt(V'*V/11);4结论与体会本次实验使用了一种相对比较简单的MATLAB语言来编写外方位元素求解,总体来说,过程比较简单,但是在数据检验方面经常不不符,追其根源是在代码书写过程中错误太多,导致结果不对。相对于C语言而言,使用MATLAB编写代码较为简单。实验二1实验任务理解根据立体相对左右倆影像的内、外方位元素和同名像点影像坐标量测值来确定该点的物方空间坐标的过程。学会利用计算机语言实验立体相对的前方交会过程,求解地面点坐标。2理论模型与方法利用点投影系数法实现一对相片的空间前方交会。212121ZNBNZYNBNYXNBNXZYX其中:112111'112122XZZXXBZBNXZZXXBZBNZXZXN和N分别为左像点和右像点投影到地面的系数。在相对定向中的上式中元素都是根据像点坐和相对元素来计算,用左右影像的外方位元素来计算,则121212SSZSSYSSXZZBYYBXXB由左右影像的外方位元素111,,和222,,计算相应的正交矩阵1R和2R,则1111111zyxRZYX,2222222zyxRZYX任一地面点坐标为:21111111211121ZNBXNZZBYNNYBYNYYYXNBXNXXXXSZSYYSSXSS3程序设计%内方位参数X10=0;Y10=0;f=120.000000/1000;%左相片外方位参数Xs1=547350.2538;Ys1=4625202.1458;Zs1=1889.1872;aa1=-0.002580;ww1=-0.003753;kk1=0.217876;%右相片外方位参数Xs2=547811.4390;Ys2=4625350.1678;Zs2=1886.4757;aa2=-0.008686;ww2=0.002983;kk2=0.212013;%左相片旋转矩阵a11=cos(aa1)*cos(kk1)-sin(aa1)*sin(ww1)*sin(kk1);a12=-cos(aa1)*sin(kk1)-sin(aa1)*sin(ww1)*cos(kk1);a13=-sin(aa1)*cos(ww1);b11=cos(ww1)*sin(kk1);b12=cos(ww1)*cos(kk1);b13=-sin(ww1);c11=sin(aa1)*cos(kk1)+cos(aa1)*sin(ww1)*sin(kk1);c12=-sin(aa1)*sin(kk1)+cos(aa1)*sin(ww1)*cos(kk1);c13=cos(aa1)*cos(ww1);R1=[a11a12a13;b11b12b13;c11c12c13];%右相片旋转矩阵a21=cos(aa2)*cos(kk2)-sin(aa2)*sin(ww2)*sin(kk2);a22=-cos(aa2)*sin(kk2)-sin(aa2)*sin(ww2)*cos(kk2);a23=-sin(aa2)*cos(ww2);b21=cos(ww2)*sin(kk2);b22=cos(ww2)*cos(kk2);b23=-sin(ww2);c21=sin(aa2)*cos(kk2)+cos(aa2)*sin(ww2)*sin(kk2);c22=-sin(aa2)*sin(kk2)+cos(aa2)*sin(ww2)*cos(kk2);c23=cos(aa2)*cos(ww2);R2=[a21a22a23;b21b22b23;c21c22c23];%以左相片投影系数为原点%利用点投影系数法进行前方交会Bx=(Xs2-Xs1);By=(Ys2-Ys1);Bz=(Zs2-Zs1);%根据像坐标计算像空间辅助坐标系坐标fori=1:d3%左相片A1=R1*[x1(i);y1(i);-f];X1=A1(1,1);Y1=A1(2,1);Z1=A1(3,1);%右相片A2=R2*[x2(i);y2(i);-f];X2=A2(1,1);Y2=A2(2,1);Z2=A2(3,1);%点投影系数N=(Bx*Z2-Bz*X2)/(X1*Z2-Z1*X2);N1=(Bx*Z1-Bz*X1)/(X1*Z2-Z1*X2);X(i)=Xs1+N*X1;Y(i)=Ys1+0.5*(N*Y1+N1*Y2+By);Z(i)=Zs1+Bz+N1*Z2;end4结论与体会