后方交会求外方位元素

西南交通大学
摄影测量学课程设计报告
单像空间后方交会
求解外方位元素
目录：
作业任务 (3)
计算原理 (3)
算法流程 (6)
源程序 (7)
计算结果 (12)
结果分析 (13)
心得体会 (13)
作业任务：
已知条件
摄影机主距f=153.24mm，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。

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

计算原理：
如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐
（X S,Y S,Z S）；像点a在像空间坐标系中的坐标是（x,y,-f）。

标依次是（X,Y,Z）、
那么由共线条件方程知，
其中a i ，b i ，c i 是只含三个独立参数ψ，ω，κ的九个方向余弦。

在方程中共有六个未知参数X S ,Y S ,Z S,ψ,ω,κ，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。

由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有
用新的符号表示各偏导数后为
s
s
s s s s s
s s s s s dZ Z y dY Y y dX X y d y d y d y y y dZ Z x dY Y x dX X x d x d x d x x x ∂∂+∂∂+∂∂+∂∂+∂∂+∂∂+=∂∂+∂∂+∂∂+∂∂+∂∂+∂∂+=κκωωϕϕκκωωϕϕ)()(
κ
ωϕκωϕd a d a d a dZ a dY a dX a y y d a d a d a dZ a dY a dX a x x s s s s s s 262524232221161514131211)()(++++++=++++++=
其中（x ）、（y ）是函数近似值，ϕd ,s s s dZ dY dX d d ,,,,κω是外方位元素近似
值的改正数，它们的系数为函数的偏导数。

为了便于推导，令
X =)()()(111s s s Z Z c Y Y b X X a -+-+- Y =)()()(222s s s Z Z c Y Y b X X a -+-+- Z
=)()()(333s s s Z Z c Y Y b X X a -+-+-
那么有
⎥⎥
⎥⎦
⎤⎢⎢⎢⎣⎡---=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⎥⎥⎥⎦⎤⎢
⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡s s s T s s s Z Z Y Y X X Z Z Y Y X X c b a c b a c b a Z Y X R 33
322211
1
对于系数，其严密算法（以
s
X x
∂∂,ϕ
∂∂x 为例）如下：
)()()(2
1322132Z X
b Z
Y
b Z X f
Z Y
f b f b X b Y b Z
X
Y b Z b Z f Z Z X X Z f x -++-=⎭⎬⎫
⎩
⎨⎧
----
=∂∂-∂∂-=∂∂ϕ
ϕϕ
[]ωκκκωκωκωωκωcos cos sin )(cos )()
(sin )(cos cos )(sin cos )()(sin )(cos cos 00000000⎭
⎬⎫⎩⎨⎧+------=⎥
⎦
⎤
⎢⎣⎡-+-----+-=f y y x x f x x y y f x x f y y x x y y f
对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为
[])(1
)(1)
()(03131312
2x x a f a Z
a Z X f f a Z X a Z a Z
f X X Z Z X X Z
f X x s
s s
-+=
-=+--=∂∂-∂∂-=∂
∂
⎪⎪
⎪⎪⎪⎭⎪
⎪
⎪
⎪
⎪
⎬⎫
-=+-=-==-=+-=-=-==-==-=x a f y f a f xy a y a f
xy
a f x f a H y a H f a a H x a a H
f
a 2622
2524161522
14232221131211),1(,,),1(,,0,0,
为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。

此时像点坐标（x,y ）作为观测值，加入相应的偶然误差改正数x v ,y v ，可列出每个点的误差方程式
y s s s y x
s s s x l d a d a d a dZ a dY a dX a v l d a d a d a dZ a dY a dX a v -+++++=-+++++=κωϕκωϕ262524232221161514131211
用矩阵表示为
⎥⎦
⎤⎢
⎣⎡=⎥⎦⎤⎢⎣⎡--=⎥⎥⎥
⎥⎥⎥
⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆∆∆∆∆=⎥⎦⎤⎢⎣⎡=2625
24
23
22
21
161514131211
00,,a a a a a a a a a a a a y y x x Z Y X v v s s s y x A l x V κωϕ 那么由l Ax V -=
)()(T
1T l A A A x -=
算法流程：
➢获取已知数据m, x0 , y0 , f, X tp, Y tp, Z tp ➢量测控制点像点坐标x,y
➢确定未知数初值X s0, Y s0, Z s0, ϕ0, ω0, κ0➢组成误差方程式并法化
➢解求外方位元素改正数
➢检查迭代是否收敛
具体的流程程序框图如下：
源程序：
程序代码：
#include<iostream>
#include<cmath>
using namespace std;
const int N=4;
const int n=6;
/*---------------矩阵相乘---------------------*/
void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2) {
int i,j,k;
for(i=0;i<i_1;i++)
for(j=0;j<j_2;j++)
{
result[i*j_2+j]=0.0;
for(k=0;k<j_12;k++)
result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];
}
return;
}
/*---------------矩阵求逆---------------------*/
void inverse(double c[n][n])
{ int i,j,h,k;
double p;
double q[n][12];
for(i=0;i<n;i++)//构造高斯矩阵
for(j=0;j<n;j++)
q[i][j]=c[i][j];
for(i=0;i<n;i++)
for(j=n;j<12;j++)
{if(i+6==j)
q[i][j]=1;
else
q[i][j]=0;}
for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据
for(i=k+1;i<n;i++)
{if(q[i][h]==0)
continue;
p=q[k][h]/q[i][h];
for(j=0;j<12;j++)
{ q[i][j]*=p;
q[i][j]-=q[k][j];
}
}
for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据
for(i=k-1;i>=0;i--)
{if(q[i][h]==0)
continue;
p=q[k][h]/q[i][h];
for(j=0;j<12;j++)
{q[i][j]*=p;
q[i][j]-=q[k][j];}}
for(i=0;i<n;i++)//将对角线上数据化为1
{ p=1.0/q[i][i];
for(j=0;j<12;j++)
q[i][j]*=p;}
for(i=0;i<n;i++) //提取逆矩阵
for(j=0;j<n;j++)
c[i][j]=q[i][j+6];
}
/*---------------矩阵转置---------------------*/
void transpose(double *m1,double *m2,int m,int n) //矩阵转置{
int i,j;
for(i=0;i<m;i++)
for(j=0;j<n;j++)
m2[j*m+i]=m1[i*n+j];
return;
}
void main()
{
double Xs,Ys,Zs,q,w,k;
double a[3],b[3],c[3];
double x0,y0,f;
double x[N],y[N];
double X[N],Y[N],Z[N];
double x1[N],y1[N];
double m;
double L[2*N];
double XX[6];
double A[2*N][6];
double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];
int i,n=0;
double sum=0,m0;
/*---------------输入点地面坐标---------------------*/
X[0]=36589.41;X[1]=37631.08;X[2]=39100.97;X[3]=40426.54;
Y[0]=25273.32;Y[1]=31324.51;Y[2]=24934.98;Y[3]=30319.81;
Z[0]=2195.17;Z[1]=728.69;Z[2]=2386.50;Z[3]=757.31;
/*---------------输入点像片坐标---------------------*/
x[0]=-86.15;x[1]=-53.40;x[2]=-14.78;x[3]=10.46;
y[0]=-68.99;y[1]=82.21;y[2]=-76.63;y[3]=64.43;
cout<<endl;
/*-----------------设定外方位元素初始值--------------*/
x0=0;y0=0;f=153.24;m=40000;
Xs=0;Ys=0;Zs=f*m/1000;
q=0;w=0;k=0;
XX[3]=1;
/*------------------迭代计算--------------------------*/
while((XX[3]>6/206265 || XX[4]>6/206265 || XX[5]>6/206265)&&n<100)
{
/*----------------旋转矩阵R-----------------------*/
a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);
a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);
a[2]=-sin(q)*cos(w);
b[0]=cos(w)*sin(k);
b[1]=cos(w)*cos(k);
b[2]=-sin(w);
c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);
c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);
c[2]=cos(q)*cos(w);
/*-----------------像点坐标计算值------------------*/
for(i=0;i<N;i++)
{
X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs);
Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs);
Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);
x1[i]=x0-f*X0[i]/Z0[i];
y1[i]=y0-f*Y0[i]/Z0[i];
}
/*-------------误差方程中各偏导数的值--------------*/
for(i=0;i<N;i++)
{
A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];
A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i];
A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i];
A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin(k))/f+f*cos(k))
*cos(w);
A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;
A[2*i][5]=y[i]-y0;
L[2*i]=x[i]-x1[i];
A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][1]=((b[1]*f+b[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][2]=((c[1]*f+c[2]*(y[i]-y0)))/Z0[i];
A[1+2*i][3]=-(x[i]-x0)*sin(w)-((y[i]-y0)*((x[i]-x0)*cos(k)-(y[i]-y0)*sin(k))/f-f*sin(k))
*cos(w);
A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;
A[1+2*i][5]=-x[i]+x0;
L[1+2*i]=y[i]-y1[i];
}
/*-------------------解法方程--------------------*/
transpose(&A[0][0],&At[0][0],2*N,6);
mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);
inverse(result1);
mult(&At[0][0],L,&result2[0][0],6,2*N,1);
mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);
Xs+=XX[0];
Ys+=XX[1];
Zs+=XX[2];
q+=XX[3];
w+=XX[4];
k+=XX[5];
n++;
}
/*----------------旋转矩阵R-----------------------*/
cout<<"迭代次数为:"<<n<<endl;
printf("\n像片外方位元素的解\n");
cout<<"航向顷角q:"<<q<<endl;
cout<<"旁向倾角w:"<<w<<endl;
cout<<"像片旋角k:"<<k<<endl;
cout<<"Xs "<<Xs<<endl;
cout<<"Ys "<<Ys<<endl;
cout<<"Zs "<<Zs<<endl;
cout<<endl;
cout<<"所有坐标已经内置,所以不显示了"<<endl; 计算结果:
X S=39795.5
Y S=27476.5
Z S=7572.69
ψ=-0.00398694
ω=0.00211388
κ=-0.067578
结果分析：
由计算结果可知，最后的计算结果是符合 6秒的限差要求的。

由于其中误差未计算，所以精度指标还不明确。

心得体会：
这次的单张航片后方交会计算外方位元素的实习让我对后方交会的原理和解算过程有了更深入的了解，同时对与之相关的一系列知识进行了一个回顾复习。