摄影测量空间后方交会

摄影测量空间后方交会

以单张影像空间后方交会方法，求解该像的外方位元素一、实验数据与理论基础：1、实验数据：航摄仪内方位元素f=153.24mm,x0=y0=0,以及4对点的影像坐标和相应的地面坐标：影像坐标地面坐标x(mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.312、理论基础空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,e,3,K。每一对像方和物方点可列出2个方程，若有3个已知地面坐标的控制点，可列出6个方程，求取外方位元素改正数△Xs,AYs,AZs,A^,A«,Ako二、数学模型和算法公式1、数学模型：后方交会利用的理论模型为共线方程。共线方程的表达公式为：qa(X-X)+b(Y-Y)+c(Z-Z)TOC\o"1-5"\h\zX=—f74S1_74S1——74S—a(X—X)+b(Y—Y)+c(Z—Z)3AS34S34Sa(X—X)+b(Y—Y)+c(Z—Z)y=—fAS2_AS2——AS—a(X—X)+b(Y—Y)+c(Z—Z)3AS3AS3AS其中参数分别为：a]=cos甲cosk一sin甲sin①sinka?=—cos9sink—sin中sinwsinka^=—sin9coswb=cos①sinkb=coswcoskb=—sinwi23C]=sin中cosk+cos中sin①sinkc?=—sin甲sink+cos甲sinwcoskc=cos9cosw旋转矩阵R为

aaa123R=bbb123ccc123」2、由于外方位元素共有6个未知数，根据上述公式可知，至少需要3个不在一条直线上的已知地面点坐标就可以求出像片的外方位元素。3、由于共线方程是非线性方程，为了便于迭代计算，需要按泰勒级数展开，取小值一次项，使之线性化，得dxdxx=(x)+~dX~dXS+dYTdys+ssy=(y)+j+令“dYsdxdx.dx.dxdxx=(x)+~dX~dXS+dYTdys+ssy=(y)+j+令“dYsaxdx+s式中，(劝(y)为函数的近似值，dX,dY,dZ,式中，(劝(y)为函数的近似值，dX,dY,dZ,d甲,do,dK为6个外方位元素的改正数。它ss们的系数为函数的偏导数。矩阵形式为:X=\dXSdYsdZsd中dodK〕aiia21a12a22a13a23a14a24a15a25a16a26」为了书写方便，令X=a(X-X)+b(Y-Y)+c(Z-Z)AS1AS1ASY=a(X-X)+b(Y-Y)+c(Z-Z)AS2AS2ASZ=a(X-X)+b(Y-Y)+c(Z-Z)AS3AS3AS则有公式：a11a12a13a21a22a23dxdXsdxda11a12a13a21a22a23dxdXsdxdYsdxdZsdydXsdydYsdydZs=只a1f=字[b1f=Z\cf+cx]=Z\af+ay]-=[bf+by]a14a15a16a24a25a26dx一x,一.一厂,=ysino———(xcosk-ysink)+fcoskcosod©Lfdxx=-fsink-——(xsink+ycosk)dofdx-dy_=-xsino-d©dy/=-fcosk-(xsink+ycosk)dofdy=-xdKf(xcosk-ysink)-fsinkcosoV=adX+adY+adZ+ad©+ado+adK-1x11S12S13S141516xV=adX+adY+adZ+ad©+ado+adK-1y21S22S23S242526y5、最后由V=AX-L、X=(AtA)-1(AtL)求得外方法元素，得到外方位元素的近似值:X=X0+AX1+AX2+..Y=Yo+AY1+AY2+...Z=Zo+AZ1+AZ2+...9=9o+A91+A92+..•①二①0+Am+A①2+...K=K0+AK1+AK2+...三、基于MATLAB程序代码1、旋转矩阵代码function[R]=Rotation(P,W,K)TO_RAD=pi/180;P=P*TO_RAD;W=W*TO_RAD;K=K*TO_RAD;a1=cos(P)*cos(K)-sin(P)*sin(W)*sin(K);a2=-cos(P)*sin(K)-sin(P)*sin(W)*cos(K);a3=-sin(P)*cos(W);b1=cos(W)*sin(K);b2=cos(W)*cos(K);b3=-sin(W);c1=sin(P)*cos(K)+cos(P)*sin(W)*sin(K);c2=-sin(P)*sin(K)+cos(P)*sin(W)*cos(K);c3=cos(P)*cos(W);R=[a1a2a3;b1b2b3;c1c2c3];2、空间后方交会代码clearall;clc;%输入控制点坐标x=[-86.15,-53.40,-14.78,10.46]/1000;y=[-68.99,82.21,-76.63,64.43]/1000;X=[36589.41,37631.08,39100.97,40426.54];Y=[25273.32,31324.51,24934.98,30319.81];Z=[2195.17,728.96,2386.50,757.31];%输入焦距f，外方位元素以及内方位元素初始值，n为迭代次数x0=0.0;y0=0.0;phi=0.0;omiga=0.0;k=0.0;m=44811.00;f=153.24/1000;X0=mean(X);Y0=mean(Y);Z0=mean(Z)+m*f;%定义最小二乘所需变量；XG=zeros(6,1);A=zeros(8,6);L=zeros(8,1);n=0;phi=phi*pi/180;omiga=omiga*pi/180;k=k*pi/180;n=n+1;%计算旋转矩阵Ra1=cos(phi)*cos(k)-sin(phi)*sin(omiga)*sin(k);a2=-cos(phi)*sin(k)-sin(phi)*sin(omiga)*cos(k);a3=-sin(phi)*cos(omiga);b1=cos(omiga)*sin(k);b2=cos(omiga)*cos(k);b3=-sin(omiga);c1=sin(phi)*cos(k)+cos(phi)*sin(omiga)*sin(k);c2=-sin(phi)*sin(k)+cos(phi)*sin(omiga)*cos(k);c3=cos(phi)*cos(omiga);R=[a1a2a3;b1b2b3;c1c2c3];%求取最小二乘中的系数矩阵内各个值以及L矩阵的值fori=1:1:4j=2*i-1;Z_Ava=a3*(X(1,i)-X0)+b3*(Y(1,i)-Y0)+c3*(Z(1,i)-Z0);A(j,1)=(a1*f+a3*x(1,i))/Z_Ava;A(j,2)=(b1*f+b3*x(1,i))/Z_Ava;A(j,3)=(c1*f+c3*x(1,i))/Z_Ava;A(j+1,1)=(a2火f+a3火y(1,i))/Z_Ava;A(j+1,2)=(b2火f+b3火y(1,i))/Z_Ava;A(j+1,3)=(c2火f+c3火y(1,i))/Z_Ava;A(j,4)=y(1,i)*sin(omiga)-(x(1,i)/f*(x(1,i)*cos(k)-y(1,i)*sin(k))+f*cos(k))*cos(omiga);A(j,5)=-f*sin(k)-x(1,i)/f*(x(1,i)*sin(k)+y(1,i)*cos(k));A(j+1,4)=-x(1,i)*sin(omiga)-(y(1,i)/f*(x(1,i)*cos(k)-y(1,i)*sin(k))-f*sin(k))*cos(omiga);A(j+1,5)=-f*cos(k)-y(1,i)/f*(x(1,i)*sin(k)+y(1,i)*cos(k));A(j+1,6)=-x(1,i);L(j,1)=x(1,i)-(x0-f*(a1*(X(1,i)-X0)+b1*(Y(1,i)-Y0)+c1*(Z(1,i)-Z0))/Z_Ava);L(j+1,1)=y(1,i)-(y0-f*(a2*(X(1,i)-X0)+b2*(Y(1,i)-Y0)+c2*(Z(1,i)-Z0))/Z_Ava);end;%根据最小得到的公式求取观测值XG=(inv(A'*A))*(A'*L);%求取地面点坐标X0=X0+XG(1,1);Y0=Y0+XG(2,1);Z0=Z0+XG(3,1);phi=phi+XG(4,1);omiga=omiga+XG(5,1);k=k+XG(6,1);、%对计算误差进行判断，在误差范围内，则继续迭代，不在误差范围内，则推出循环。while(XG(4,1)>=6.0/206265.0||XG(5,1)>=6.0/206265.0||XG(6,1)>=6.0/206265.0)n=n+1;a1=cos(phi)*cos(k)-sin(phi)*sin(omiga)*sin(k);a2=-cos(phi)*sin(k)-sin(phi)*sin(omiga)*cos(k);a3=-sin(phi)*cos(omiga);b1=cos(omiga)*sin(k);b2=cos(omiga)*cos(k);b3=-sin(omiga);c1=sin(phi)*cos(k)+cos(phi)*sin(omiga)*sin(k);c2=-sin(phi)*sin(k)+cos(phi)*sin(omiga)*cos(k);c3=cos(phi)*cos(omiga);R=[a1a2a3;b1b2b3;c1c2c3];fori=1:1:4j=2*i-1;Z_Ava=a3*(X(1,i)-X0)+b3*(Y(1,i)-Y0)+c3*(Z(1,i)-Z0);A(j,1)=(a1*f+a3*x(1,i))/Z_Ava;A(j,2)=(b1*f+b3*x(1,i))/Z_Ava;A(j,3)=(c1*f+c3*x(1,i))/Z_Ava;A(j+1,1)=(a2火f+a3火y(1,i))/Z_Ava;A(j+1,2)=(b2火f+b3火y(1,i))/Z_Ava;A(j+1,3)=(c2火f+c3火y(1,i))/Z_Ava;A(j,4)=y(1,i)*sin(omiga)-(x(1,i)/f*(x(1,i)*cos(k)-y(1,i)*sin(k))+f*cos(k))*cos(omiga);A(j,5)=-f*sin(k)-x(1,i)/f*(x(1,i)*sin(k)+y(1,i)*cos(k));A(j,6)=y(1,i);A(j+1,4)=-x(1,i)*sin(omiga)-(y(1,i)/f*(x(1,i)*cos(k)-y(1,i)*sin(k))-f*sin(k))*cos(omiga);A(j+1,5)=-f*cos(k)-y(1,i)/f*(x(1,i)*sin(k)+y(1,i)*cos(k));A(j+1,6)=-x(1,i);L(j,1)=x(1,i)-(x0-f*(a1*(X(1,i)-X0)+b1*(Y(1,i)-Y0)+c1*(Z(1,i)-Z0))/Z_Ava);L(j+1,1)=y(1,i)-(y0-f*(a2*(X(1,i)-X0)+b2*(Y(1,i)-Y0)+c2*(Z(1,i)-Z0))/Z_Ava);end;%完成最终迭代，输出地面点坐标及旋转矩阵XG=(inv(A'*A))*(A'*L);X0=X0+XG(1,1);Y0=Y0+XG(2,1);Z0=