摄影测量空间后方交会.docx
- 文档编号:25644192
- 上传时间:2023-06-11
- 格式:DOCX
- 页数:10
- 大小:201.80KB
摄影测量空间后方交会.docx
《摄影测量空间后方交会.docx》由会员分享,可在线阅读,更多相关《摄影测量空间后方交会.docx(10页珍藏版)》请在冰豆网上搜索。
摄影测量空间后方交会
摄影测量空间后方交会
以单张影像空间后方交会方法,求解该像的外方位元素
一、实验数据与理论基础:
1、实验数据:
航摄仪内方位元素f=153.24mm,x0=y0=0,以及4对点的影像坐标和相应的地面坐标:
影像坐标
地面坐标
x(mm)
y(mm)
X(m)
Y(m)
Z(m)
1
-86.15
-68.99
36589.41
25273.32
2195.17
2
-53.40
82.21
37631.08
31324.51
728.69
3
-14.78
-76.63
39100.97
24934.98
2386.50
4
10.46
64.43
40426.54
30319.81
757.31
2、理论基础
(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=[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;
%计算旋转矩阵R
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];
%求取最小二乘中的系数矩阵内各个值以及L矩阵的值
fori=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=[a1a2a3;b1b2b3;c1c2c3];
fori=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);
end;
formatshortg
R
formatlongg
X0,Y0,Z0
四、实验结果
五、实验体会
本次作业通过实验让我更加清楚地的认识到什么是空间后方交会,并在完成作业的同时巩固了关于空间后方交会基本公式和误差方程式以及共线方程式的知识。
同时由于编程基础薄弱,因此基于MATLAB的程序编写,基本参照老师给的文档,不过在自己动手去编的过程中,还是收获很多的,增强了个人的动手能力,希望以后可以通过自己编出程序。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 摄影 测量 空间 后方 交会