摄影测量（tip1）：空间前方交会与空间后方交会

deja_vu97

已于 2022-05-26 01:36:44 修改

阅读量1.3w

点赞数 11

分类专栏：摄影测量文章标签： matlab 算法

于 2022-05-17 18:30:45 首次发布

本文链接：https://blog.csdn.net/deja_vu97/article/details/124819665

版权

摄影测量专栏收录该内容

4 篇文章 11 订阅

订阅专栏

1. 前方交会

1.1 定义

已知 A, B 平面坐标，在已知点 A, B 分别瞄准未知点 P，从而求出点 P 平面坐标

1.2 分类

-角度前方交会-

通过观测角 A 和 B 求出点 P 坐标的定位方法

-方位角前方交会-

通过观测方位角 $\alpha_{AP}$ 和 $\alpha_{BP}$ 求出点 P 坐标的定位方法

1.3 定位精度

点P的精度为

$m_{P}=m_{\alpha}\frac{\sqrt{2}}{sinP}\sqrt{AP^2+BP^2}$

可见，为了减少测量误差，有效途径为

1、提高方位角观测精度，即减小 mα ；

2、P 点离 A、B 两点不能太远，即减小 AP 和 BP；

3、角 P 最好接近90度，即 sinP 的绝对值尽量大。

2. 后方交会

2.1 定义

已知 A、B 的平面坐标，在未知点 P分别瞄准已知点 A, B，从而求出点 P 平面坐标

2.2 方法

已知 A、B、C三点的坐标，通过测量两个角度 α、β（γ是优角α+β所对的劣角），即可求出这三个角度顶点 P 的坐标。

危险圆：

点 P 在三角形 ABC 的外接圆上时，α,β,γ 将保持不变。如此一来，点 P 的坐标将有无穷个——外接圆上的任意一点均可以是点 P。此时，使用计算公式计算点 P 坐标时，可能会因为除以零而得到无效解。

综上，点 P 靠近外接圆时，很小的观测误差都会引起点 P 位置的较大偏差。因此，称三角形 ABC 的外接圆为危险圆。后方交会时，应避免点 P 离危险圆很近。

3. 侧方交会

已知 A、B 的平面坐标，在某个已知点 A 分别瞄准 B, P，从而求出点 P 平面坐标

4. 空间后方交会

我们可以利用地面控制点的大地坐标及其在像片上的像点，通过空间后方交会的方法来获得摄影测量像片的6个外方位元素。

4.1 原理

共线条件方程

$x=-f\frac{a_1(X-X_S)+b_1(Y-Y_S)+c_1(Z-Z_S)}{a_3(X-X_S)+b_3(Y-Y_S)+c_3(Z-Z_S)}$

$y=-f\frac{a_2(X-X_S)+b_2(Y-Y_S)+c_2(Z-Z_S)}{a_3(X-X_S)+b_3(Y-Y_S)+c_3(Z-Z_S)}$

其中， $(x,y,-f)$ 为某控制点对应像点在像空间坐标系中的坐标， $(X,Y,Z)$ 为某控制点对应像点在像空间辅助坐标系中的坐标； $(X_S,Y_S,Z_S)$ 为待求的外方位线元素， $\begin{bmatrix} a_1 & a_2 & a_3\\ b_1 & b_2 & b_3\\ c_1 & c_2 & c_3 \end{bmatrix}$ 为待求的外方位角元素对应的旋转矩阵（方向余弦矩阵）。

4.2 控制点

由 4.1 可知，一个地面控制点可以列出2个共线方程，单张像片具有6个外方位元素 $X_S,Y_S,Z_S,\varphi,\omega,\kappa$ ，因此对于一张航片，我们需要至少3个控制点来解算外方位元素。

4.3 matlab解算过程

代码如下

% 利用间接平差进行空间后方交会，计算外方位元素
 
clear;clc;
%读取数据
xx=xlsread('空间后方交会已知点坐标.xlsx',1,'B3:B6'); %mm
yy=xlsread('空间后方交会已知点坐标.xlsx',1,'C3:C6'); %mm
XX=xlsread('空间后方交会已知点坐标.xlsx',1,'D3:D6'); %m
YY=xlsread('空间后方交会已知点坐标.xlsx',1,'E3:E6'); %m
ZZ=xlsread('空间后方交会已知点坐标.xlsx',1,'F3:F6'); %m
 
%已知参数
ff=153.24; %焦距(mm)
xxo=0;     %mm
yyo=0;     %mm  
mb=40000; %航摄比例尺
 
%框标坐标系->像空间坐标系
xx=xx-xxo;
yy=yy-yyo;
 
%根据共线方程求间接平差 系数矩阵B 和 常数阵l
syms ph ka om Xs Ys Zs X Y Z f x y
% %旋转矩阵参数
Rph=[cos(ph) 0 -sin(ph);0 1 0;sin(ph) 0 cos(ph)];
Rom=[1 0 0;0 cos(om) -sin(om);0 sin(om) cos(om)];
Rka=[cos(ka) -sin(ka) 0;sin(ka) cos(ka) 0;0 0 1];
R=Rph*Rom*Rka;
a1=R(1,1);a2=R(1,2);a3=R(1,3);
b1=R(2,1);b2=R(2,2);b3=R(2,3);
c1=R(3,1);c2=R(3,2);c3=R(3,3);
% %泰勒展开常数项
x0=-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs))./(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs));
y0=-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs))./(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs));
% %系数阵B
bx_ph=diff(x0,ph);
bx_om=diff(x0,om);
bx_ka=diff(x0,ka);
bx_Xs=diff(x0,Xs);
bx_Ys=diff(x0,Ys);
bx_Zs=diff(x0,Zs);
by_ph=diff(y0,ph);
by_om=diff(y0,om);
by_ka=diff(y0,ka);
by_Xs=diff(y0,Xs);
by_Ys=diff(y0,Ys);
by_Zs=diff(y0,Zs);
b=[ bx_ph,bx_om,bx_ka,bx_Xs,bx_Ys,bx_Zs;
    by_ph,by_om,by_ka,by_Xs,by_Ys,by_Zs];
% %常数阵l
l=[x-x0;y-y0];
 
% %赋值常数：X Y Z f x y
B=[];
L=[];
for i=1:4  %4个已知点数据
    bb=subs(b,[f,x,y,X,Y,Z],[ff,xx(i),yy(i),XX(i),YY(i),ZZ(i)]);
    B=[B;bb];
    ll=subs(l,[f,x,y,X,Y,Z],[ff,xx(i),yy(i),XX(i),YY(i),ZZ(i)]);
    L=[L;ll];
end
 
%% 平差
%初始化 ph ka om Xs Ys Zs 
phi=0;
omega=0;
kappa=0;
XXs=mean(XX);
YYs=mean(YY);
ZZs=mean(ZZ); %单位m 已知数据为大地坐标，因此不用航高作为Z初始值
 
%迭代平差
x_hat=[0,0,0,mean(XX),mean(YY),mean(ZZ)];
times=0; %迭代次数
v_max=10;%改正数最大值
while(v_max>=1e-6)
    B0=double(subs(B,[ph,om,ka,Xs,Ys,Zs],[phi,omega,kappa,XXs,YYs,ZZs]));
    L0=double(subs(L,[ph,om,ka,Xs,Ys,Zs],[phi,omega,kappa,XXs,YYs,ZZs]));
    %平差
    dx=(B0'*B0)\(B0'*L0); %等精度观测，权阵P为单位阵I
    x_hat=x_hat+dx;
    %迭代参数计算
    v_max=max(abs(dx));
    times=times+1
    %下次迭代自变量
    phi=x_hat(1);
    omega=x_hat(2);
    kappa=x_hat(3);    
    XXs=x_hat(4);
    YYs=x_hat(5);
    ZZs=x_hat(6);
end
%验证结果
disp('检验计算:')
L0=double(subs(L,[ph,om,ka,Xs,Ys,Zs],[phi,omega,kappa,XXs,YYs,ZZs]))
%输出参数
disp('外方位元素平差值:')
phi   %弧度
omega %弧度
kappa %弧度
XXs   %米
YYs   %米
ZZs   %米
%精度评定
disp('精度评定:')
disp('参数平差值方差阵:')
V=B0*dx-L0;
r=8-3;%多余观测数
sigma0=(V'*V)/r; %单位权方差
format short
Dxx=sigma0*inv(B0'*B0) %参数平差值方差阵

代码中待解算的原始数据为四个控制点的框标坐标和地面坐标。