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

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为单位阵Ix_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) %参数平差值方差阵

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

5. 空间前方交会

已知值：内方位元素 $(x_0,y_0,f)$ ，以及通过空间后方交会求出的外方位元素 $X_S,Y_S,Z_S,\varphi,\omega,\kappa$ ；某个立体像对左右片的像点坐标 $a_1(x_1,y_1), a_2(x_2,y_2)$

立体像对空间前方交会待求量：立体像对 $a_1,a_2$ 的对应点在地面摄影测量坐标系中的坐标A

5.1 原理

5.2 计算步骤

1. 利用外方位角元素计算左右像片的旋转矩阵 $R_1,R_2$

2. 利用外方位线元素计算左右像片构成的摄影基线分量 $B_u,B_v,B_w$

$B_u=X_{S2}-X_{S1}$

$B_v=Y_{S2}-Y_{S1}$

$B_w=Z_{S2}-Z_{S1}$

3. 逐点计算各像点的像空间辅助坐标 $(u_1,v_1,w_1),(u_2,v_2,w_2)$

$\begin{bmatrix} u_1\\ v_1\\ w_1 \end{bmatrix} = R_1 \begin{bmatrix} x_1\\ y_1\\ -f \end{bmatrix}$

$\begin{bmatrix} u_2\\ v_2\\ w_2 \end{bmatrix} = R_2 \begin{bmatrix} x_2\\ y_2\\ -f \end{bmatrix}$

4. 逐点计算各像点的投影系数

$N_1=\frac{B_uw_2-B_wu_2}{u_1w_2-u_2w_1}$

$N_2=\frac{B_uw_1-B_wu_1}{u_1w_2-u_2w_1}$

5. 逐点计算各像点的地面摄影测量坐标系的坐标

$X=X_{S1}+N_1u_1=X_{S2}+N_2u_2$

$Y=Y_{S1}+N_1v_1=Y_{S2}+N_2v_2$

$Z=Z_{S1}+N_1w_1=Z_{S2}+N_2w_2$

另，同名光线不可能在空间完全相交，为确保空间光线的坐标之差仅在Y方向，最后只需对Y1和Y2取平均（上下视差），即，

$Y=\frac{1}{2}[(Y_{S1}+N_1v_1)+(Y_{S2}+N_2v_2)]$