攝影測(cè)量空間后方交會(huì)

攝影測(cè)量空間后方交會(huì)以單張影像空間后方交會(huì)方法，求解該像的外方位元素一 、實(shí)驗(yàn)數(shù)據(jù)與理論基礎(chǔ)：1、實(shí)驗(yàn)數(shù)據(jù)：航攝儀內(nèi)方位元素f153.24mm，x0y00，以及4對(duì)點(diǎn)的影像坐標(biāo)和相應(yīng)的地面坐標(biāo)：影像坐標(biāo)地面坐標(biāo)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、理論基礎(chǔ)(1) 空間后方交會(huì)是以單幅影像為基礎(chǔ)，從該影像所覆蓋地面范圍內(nèi)若干控制點(diǎn)的已知地面坐標(biāo)和相應(yīng)點(diǎn)的像坐標(biāo)量測(cè)值出發(fā)，根據(jù)共線條件方程，解求該影像在航空攝影時(shí)刻的外方位元素Xs，Ys，Zs，。(2) 每一對(duì)像方和物方點(diǎn)可列出2個(gè)方程，若有3個(gè)已知地面坐標(biāo)的控制點(diǎn)，可列出6個(gè)方程，求取外方位元素改正數(shù)Xs，Ys，Zs，。二、數(shù)學(xué)模型和算法公式1、數(shù)學(xué)模型：后方交會(huì)利用的理論模型為共線方程。共線方程的表達(dá)公式為： 其中參數(shù)分別為： 旋轉(zhuǎn)矩陣R為2、 由于外方位元素共有6個(gè)未知數(shù)，根據(jù)上述公式可知，至少需要3個(gè)不在一條直線上的已知地面點(diǎn)坐標(biāo)就可以求出像片的外方位元素。3、由于共線方程是非線性方程，為了便于迭代計(jì)算，需要按泰勒級(jí)數(shù)展開(kāi)，取小值一次項(xiàng)，使之線性化，得式中，(x),(y)為函數(shù)的近似值，為6個(gè)外方位元素的改正數(shù)。它們的系數(shù)為函數(shù)的偏導(dǎo)數(shù)。矩陣形式為： 為了書(shū)寫(xiě)方便，令則有公式： 4、 計(jì)算中，通常將控制點(diǎn)的地面坐標(biāo)視為真值，而把相應(yīng)的像點(diǎn)坐標(biāo)視為觀測(cè)值，加入相應(yīng)的改正數(shù) ，得 ，可列出每個(gè)點(diǎn)的誤差方程式：5、 最后由、求得外方法元素，得到外方位元素的近似值：三、 基于MATLAB程序代碼1、旋轉(zhuǎn)矩陣代碼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、空間后方交會(huì)代碼clear all;clc; %輸入控制點(diǎn)坐標(biāo)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èi)方位元素初始值，n為迭代次數(shù)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;%計(jì)算旋轉(zhuǎn)矩陣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;%求取最小二乘中的系數(shù)矩陣內(nèi)各個(gè)值以及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;%根據(jù)最小得到的公式求取觀測(cè)值XG=(inv(A*A)*(A*L); %求取地面點(diǎn)坐標(biāo)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);、%對(duì)計(jì)算誤差進(jìn)行判斷，在誤差范圍內(nèi)，則繼續(xù)迭代，不在誤差范圍內(nèi)，則推出循環(huán)。while(XG(4,1)=6.0/.0|XG(5,1)=6.0/.0|XG(6,1)=6.0/.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;%完成最終迭代，輸出地面點(diǎn)坐標(biāo)及旋轉(zhuǎn)矩陣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(