摄影测量空间后方交会

以单张影像空间后方交会方法，求解该像的外方位元素。

一、实验数据与理论基础：

1、实验数据：

航摄仪内方位元素f＝153.24mm，x0＝y0＝0，以及4对点的影像坐标和相应的地面坐标：

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 = 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;

计算旋转矩阵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=[a1 a2 a3;b1 b2 b3;c1 c2 c3];

求取最小二乘中的系数矩阵内各个值以及l矩阵的值。

for i=1:1:4

j=2*i-1;

z_**a=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_**a;

a(j,2)=(b1*f+b3*x(1,i))/z_**a;

a(j,3)=(c1*f+c3*x(1,i))/z_**a;

a(j+1,1)=(a2*f+a3*y(1,i))/z_**a;

a(j+1,2)=(b2*f+b3*y(1,i))/z_**a;

a(j+1,3)=(c2*f+c3*y(1,i))/z_**a;

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_**a);

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_**a);

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)