摄影测量空间后方交会.docx

摄影测量空间后方交会以单张影像空间后方交会方法，求解该像的外方位元素一 、实验数据与理论基础：1、实验数据：航摄仪内方位元素f153.24mm，x0y00，以及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、理论基础(1) 空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，。(2) 每一对像方和物方点可列出2个方程，若有3个已知地面坐标的控制点，可列出6个方程，求取外方位元素改正数Xs，Ys，Zs，。二、数学模型和算法公式1、数学模型：后方交会利用的理论模型为共线方程。共线方程的表达公式为： 其中参数分别为： 旋转矩阵R为2、 由于外方位元素共有6个未知数，根据上述公式可知，至少需要3个不在一条直线上的已知地面点坐标就可以求出像片的外方位元素。3、由于共线方程是非线性方程，为了便于迭代计算，需要按泰勒级数展开，取小值一次项，使之线性化，得式中，(x),(y)为函数的近似值，为6个外方位元素的改正数。它们的系数为函数的偏导数。矩阵形式为： 为了书写方便，令则有公式： 4、 计算中，通常将控制点的地面坐标视为真值，而把相应的像点坐标视为观测值，加入相应的改正数 ，得 ，可列出每个点的误差方程式：5、 最后由、求得外方法元素，得到外方位元素的近似值：三、 基于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 = a1 a2 a3;b1 b2 b3;c1 c2 c3;2、空间后方交会代码clear all;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=a1 a2 a3;b1 b2 b3;c1 c2 c3;%求取最小二乘中的系数矩阵内各个值以及L矩阵的值for i=1:1:4 j=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=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=a1 a2 a3;b1 b2 b3;c1 c2 c3;for i=1:1:4 j=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=Z0+XG(3,1);phi=phi+XG(4,1);omiga=omiga+XG(5,1);k=k+XG