后方交会求外方位元素

1 / 13 西南交通大学 摄影测量学课程设计报告

单像空间后方交会 求解外方位元素 2 / 13

目录：

作业任务…………………………………………………………3 计算原理…………………………………………………………3 算法流程…………………………………………………………6 源程序……………………………………………………………7 计算结果…………………………………………………………12 结果分析…………………………………………………………13 心得体会…………………………………………………………13 3 / 13

作业任务： 已知条件 摄影机主距f=153.24mm，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。

以单像空间后方交会方法，求解该像片的外方位元素。 计算原理： 如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。那么由共线条件方程知，

点号 像点坐标 地面坐标 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

)()()()()()()()()()()()(33322203331110ssssssssssssZZcYYbXXaZZcYYbXXafyyZZcYYbXXaZZcYYbXXafxx

 4 / 13

其中ai，bi，ci是只含三个独立参数ψ，ω，κ的九个方向余弦。

在方程中共有六个未知参数XS,YS,ZS,ψ,ω,κ，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有

用新的符号表示各偏导数后为 ssssss

ssssss

dZZydYYydXXydydydyyydZZxdYYxdXXxdxdxdxxx)(

)( 5 / 13

dadadadZadYadXayydadadadZadYadXaxxssssss262524232221161514131211)()(

其中（x）、（y）是函数近似值，

d

,

sssdZdYdXdd,,,,是外方位元素近似

值的改正数，它们的系数为函数的偏导数。为了便于推导，令 X=)()()(111sssZZcYYbXXa

Y=)()()(222sssZZcYYbXXa Z=)()()(333sssZZcYYbXXa

那么有





sssTsssZZYYXXZZYYXXcbacbacbaZYXR

333222

111

对于系数，其严密算法（以sXx,x为例）如下：

)()()(21322132ZXbZYbZXfZYfbfbXbYbZXYbZbZfZZXXZfx





coscossin)(cos)()(sin)(coscos)(sincos)()(sin)(coscos00000000fyyxxfxxyyfxxf

yyxxyyf

对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为

)(1)(1)()(031313122xxafaZaZXffaZXaZaZfXXZZXXZfXxsss



 6 / 13







xafyfafxyayafxyafxfaHyaHfaaHxaaHfa2622252416152214232221131211),1(,,),1(,,0,0,

为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数xv,yv，可列出每个点的误差方程式

ysssyxsssxldadadadZadYadXavldadadadZadYadXav262524232221

161514131211

用矩阵表示为



26252423222116151413121100,,aaaaaaaaaaaayyxxZYXvvsssyxA

lxV

那么由lAxV )()(T1TlAAAx

算法流程： 7 / 13

➢ 获取已知数据 m, x0 , y0 , f , Xtp, Ytp, Ztp

➢ 量测控制点像点坐标 x,y

➢ 确定未知数初值 Xs0, Ys0, Zs0, 0, 0, 0

➢ 组成误差方程式并法化

➢ 解求外方位元素改正数 ➢ 检查迭代是否收敛 具体的流程程序框图如下：

是 否 否 是

否 是

源程序：

输入原始数据 像点坐标计算，系统误差改正 确定外方位元素初始值

组成旋转矩阵R 逐点组成误差方程式并法化 所以像点完否 解法方程，求外方位元素改正数 计算改正后的外方位元素 外方位元素改正数是否小于限差 输出计算成果，计算结束

迭代次数小于

n

结束并显示错误信息 8 / 13

程序代码：

#include #include 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;ifor(j=0;j{ result[i*j_2+j]=0.0; for(k=0;kresult[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;ifor(j=0;jq[i][j]=c[i][j]; for(i=0;ifor(j=n;j<12;j++) {if(i+6==j) q[i][j]=1; else q[i][j]=0;} for(h=k=0;kfor(i=k+1;i{if(q[i][h]==0) continue; p=q[k][h]/q[i][h]; for(j=0;j<12;j++) { q[i][j]*=p; 9 / 13

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{ p=1.0/q[i][i]; for(j=0;j<12;j++) q[i][j]*=p;} for(i=0;ifor(j=0;jc[i][j]=q[i][j+6]; }

/*---------------矩阵转置---------------------*/ void transpose(double *m1,double *m2,int m,int n) //矩阵转置 { int i,j; for(i=0;ifor(j=0;jm2[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];