摄影测量学课程设计报告

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

连续相对定向方法

计算两张像片的相对定向元素

1 / 18

目录：

作业任务……………………………………………………3

计算原理……………………………………………………3

算法流程……………………………………………………13

源程序………………………………………………………14

计算结果……………………………………………………17

结果分析……………………………………………………18

心得体会……………………………………………………18

作业任务：

2 / 18

已知条件：

在一个航带内相邻两张像片上分别量测了6对同名点的像平面坐标，数据如下。

左片 点号 x(mm) 1 2 3 4 5 6 f=24mm 请计算：

采用连续像对相对定向方法，两张像片的相对定向元素。

-6.091 7.098 4.538 6.858 -10.050 -8.0 y(mm) 1.983 0.924 1.068 1.208 -0.514 1.293 x(mm) -5.5 7.694 5.098 7.429 -9.152 -7.441 y(mm) -3.202 -2.830 -2.878 -2.578 -5.2 -3.981 右片 计算原理：

一． 立体像对的相对定向元素

用于描述立体像对中两张像片相对关系的参数，称为相对定向元素。根据所取像空间辅助坐标系的不同，相对定向元素常有两种不同的描述方式。

1.连续像对的相对定向元素

连续像对的相对定向系统，是以左片为基准，求出右片相对于

3 / 18

左片的相对方位元素。以左摄站为原点，建立与左片像空间坐标系一致的像空间辅助坐标系S1X1Y1Z1，右片的像空间辅助坐标系S2X2Y2Z2与S1X1Y1Z1平行，如图1-1所示，在S1X1Y1Z1坐标系中。两张像片是12个方位元素为 左像片： 右像片： XS1YS1ZS1011102,2,2

XS2BX,YS2BY,ZS2BZ（或像空间辅助坐标系）的角方位元2,2,2为右片相对于左片

素；BX为摄影基线的X方向分量，由于X轴接近于摄影基线，

BX远大于BY和BZ，因而可以认为BX只决定模型的比例尺，而

与两张像片的相对关系无关。这样，除BX之外的五个非零元素

BY，BZ，2,2,2可确定两张像片的相对关系，作为连续像对

的相对定向元素。

Z2 S2 Y2 X2

Z1(z1)

S1 Y1(y1) BBZ BX BY2 2 X1(x1) 图1-1连续像对的相对定向元素

连续像对的相对定向系统的特点是以左片为参照，通过解算右片相对

于左片的五个方位元素来确定两张像片之间的相对关系，建立立体模型。

4 / 18

2.单独像对的相对定向元素

单独像对的相对系统以左摄站S1为原点，摄影基线B为X轴，

在左主核面内过S1且垂直於X轴的直线为Z轴，建立像空间辅助坐标系。这时两像片内的12个方位元素可表示为

左像片： 右像片：

XS1YS1ZS101,1,10

XS2B,YS2ZS202,2,2

同样，除B外的五个非零元素1,1,2,2,2确定两像片的相对关系，称为单独像对的相对定向元素。 二、 解析法相对定向原理 1. 相对定向的共面条件

图1-2

在图1-2中，S1a1和S2a2为一条同名光线，与摄影基线B位于同一核面内，即S1a1、S2a2和B三条直线共面。若此，则它们

5 / 18

对应矢量的混合积为零，即 B(S1a1S2a2)0 三矢量在像空间辅助坐标系中的坐标分别为(BX,BY,BZ)，

(X1,Y1,Z1)和(X2,Y2,Z2)，则共面条件方程表示为： BXFX1X2BYY1Y2BZZ10 Z2共面条件方程是否成立是完成相对定向的标准。解析相对定向就是根据共面条件方程解求相对定向元素。 2. 连续像对的相对定向

连续像对的相对定向是以左像片为基准，求出右像片相对于左像片的五个定向元素BY,BZ,2,2,2。在相对定向解析计算时，通常把摄影基线B改写为b，b成为投影基线。这里

Bmb

其中：m为摄影比例尺分母；bX，bY和bZ为投影基线对应的分量。为了统一单位，把bY和bZ两个基线元素改为角度形式表示，如图1-2所示。由图可知，有：

bYbXtanbX （1-3） bXbtanbZZcos式中：和为基线的偏角和倾角。将式（1-3）代入共面条件方程得

bXFX1X2（1-4）

bXbXY1Y21Z10 Z2Z1bXX1Y1X2Y2Z2式（1-4）中含有五个相对定向因素，其中2,2,2隐含在

6 / 18

该式是一个非线性函数。为了平差计算，将式（1-4）(X2,Y2,Z2)中，

按多元函数泰勒级数展开，取小值一次项，得共面方程的线形公式为：

FF0FFFFFddddd0 (1-5) 式中：F0为函数F的近似值，同时为了书写方便去除了角元素的下标。式中的偏导数计算为

01FbXX1Y1X2Y2FbXX1Y1X2Y2000Z1bX(Z1X2Z2X1) Z21Z1bX(X1Y2X2Y1) Z2FbXX1X20Y1Y2XZ12bXY1Z2Z11Z21Y2 bXbXX1Z1X1Y1 对求偏导得：

F(Z2)bX(Z1Y1)X2bX(Y1X1) bXY1X2bXX1X2bXZ1Z2bXZ2Y1

同理可得：

X2Y2sin Y2X2sinZ2cos 7 / 18

Z2Y2cos  以及

Y21Z21FX2bXbXbX

YZXZXY111111bXYY12bXX1Y2bXZ1Z2bXX1Z2

类似得：

Y21Z21FX2bXbXbX

Y1Z1X1Z1X1Y1bXX2Z1bXZ1Y2bXX1X2bXYY12

将各偏导数代入式（5-4-5），舍去含有和的二次小项，只保留一次小项，同

时等式两边同除以bX得

(Z1X2X1Z2)d(X1Y2Y1X2)dY1X2d(YY12Z1Z2)d

X2Z1dF00 （1-6） bX 顾及点投影系数得：

Z1X2X1Z2bXZ1bZX1bbbXZ1ZX1XZ1 N2N2bXN2X1Y2Y1X2bXY1bYX1bXbYbXYXY1 11N2N2bXN2N2，并近似地取Z2 代入式（1-6），等式两边同乘以Y1Y2,Z1Z2，则式（1-6）可化简为：

F0N2Y2X2Y2Y22bXdbXdN2dZ2NdXNd0 222Z2Z2ZbZ2X2 令QF0N2 bXZ28 / 18

最后得：

Y2X2Y2Y22QbXdbXdN2dZ2N2dX2N2d（1-7） Z2Z2Z2

式（1-7）即为连续法相对定向的解析计算公式。式中

bXX1QbYY1bZZ1XY2Z2F0N2F02

X1Z1bXZ2X1Z2Z1X2X2Z2X1X2X1X2Z1bYbXbZbbZY2XZ2X2Z2X1Z1N1Y1N2Y2bY （1-8） Z1X1Z1X1Z1Y1Z2X2Z2X2Z2 式中：N1Y1为左片投影点在以左摄站为原点的像空间辅助

坐标系中的坐标；N2Y2为右摄站的像空间辅助坐标系中的坐标；bY为两摄站的Y坐标之差。所以Q的几何意义是模型上同名点的Y坐标之差，称为上下视差。

3、单独像对的相对定向

单独像对相对定向时，基线b作为像空间辅助坐标的X轴，bx=b，by=bz=0，相对定向元素为1,1,2,2,2。此时共面

b00条件方程为FX1Y1Z1bX2Y2Z2Y1Z1Y2Z20。

三、相对定向元素的解算

由于存在多余观测，根据最小二乘平差原理，将上下视差Q作为

9 / 18

观测值，可以写出误差方程式，即

Y2X2Y2Y22（1-9） QbXdbXdN2dZ2N2dX2N2dQ Z2Z2Z2用一般符号表示误差方程式为：

adbdcdddedl（1-10）

用矩阵表示误差方程式为

ddedl dd

abcd如在一个像对中量测了n对像点，则可以列出n个误差方程式，即

da1b1c1d1e11l1a2b2c2d2e2dld2 dnanbncndnendln写成一般形式为V式中：

AXL（1-13）

10 / 18

V12na1a2AanddXdddl1lL2lnb1c1Te1b2c2d2e2bncndnend1

相对应的法方程为

APAXAPL0

TT一般情况下像点坐标为等权观测，权阵P是单位矩阵，法方程可

TTAAXAL0 化简为

T1TX(AA)AL 法方程的解为

解X即为相对定向元素近似值的改正数。由于误差方程式是根据泰勒级数展开的一次项近似公式，因此定向元素要用迭代法解求，具体过程如下：

（1） 原始数据的输入及像点坐标的预处理； （2） 确定相对定位元素的初始值；

（3） 计算右片的方向余弦值，组成旋转矩阵R1，计算左片各像点

的像空间辅助坐标(X,Y1,Z1)；

1（4） 计算右片的方向余弦值，组成旋转矩阵R2，计算基线分量bY和

bZ；

（5） 计算右片各像点的像空间辅助坐标(X2,Y2,Z2)，计算各像点的

11 / 18

点投影系数N1,N2和上下视差Q；

（6） 逐点组成误差方程式并法化，完成法方程系数矩阵和常数项

矩阵的计算；

（7） 解法方程，求出相对定位元素的改正数； （8） 计算相对定位元素的新值：

0d,0d,0d,0d,0d

4（9） 检查所有改正数是否小于限值0.310rad，如满足条件，则结

束相对定向计算。重复（4）~（9）。以上步骤的程序框图如图1-4所示。

一、 模型点坐标的计算

当完成相对定向，正确求解出相对定向元素后，就可用空间前方交会公式计算出模型点的坐标，建立与地面相似的数字立体模型，以左摄站为原点，其大小和方位均是任意的。对于任一模型点，有

X1x1YRy111Z1fN1X2x2YRy22 2Z2fbXZ2bZX2bZbXN2X1Z1

X1Z2X2Z1X1Z2X2Z1对于单独像对的相对定向，有bXbZ0，相应的点投影系数为

N1bZ2bZ1 N2X1Z2X2Z1X1Z2X2Z1模型内左、右摄站任一模型点在像空间辅助坐标系中的坐标： 左摄站坐标为：

12 / 18

XS10YS10（1-12） ZS10右摄站坐标为：

XS2XS1bXbXYS2YS1bYbY（1-13） ZS2ZS1bZbZ任一模型点坐标：

XmXS1N1X1N1X1Y1YNYYNY11S222m2S11 YN1Y1N2Y2bYS1212N1Y1N2Y2bYZZS1N1Z1N1Z1m

算法流程：

具体的流程图如下：

输入原始数据及预处理 确定初始值bX(x1x2)1，000000 R1和计算左片(X,Y,Z)计算bY、bZ、右片的R2 计算右片各像点的(X2,Y2,Z2) 计算各点N1,N2和Q 逐点组成误差方程式并法化 否 所有法化点完否 13 / 18

解法方程，求改正数d,d,d,d,d 计算定向元素新值 否 改正数是否小于限差 计算结束

源程序：

#include #include usingnamespace std; constint N=8;

void Mult(double R[8][8], double b[8][8],double c[8][8], int m, int t, int p); void inverse(double c[N][N],int t); void main()

{double R[N][N],lx[N],ly[N],rx[N],ry[N],A[N][N],l[N][N],ATA[N][N]; double f=0.024,Q=0,W=0,K=0,x=0,y=0,X[8][8],bu,bv,bw; int i,j,n; int q=0,k;

cout<<\"请输入同名像点的对数\"<<\" \"; cin>>n;

lx[1]=0.001983,lx[2]=0.000924,lx[3]=0.001068,lx[4]=0.001208,lx[5]=-0.000514,lx[6]=0.001293, ly[1]=-0.006091,ly[2]=0.007098,ly[3]=0.004538,ly[4]=0.006858,ly[5]=-0.010050,ly[6]=-0.0080, rx[1]=-0.003202,rx[2]=-0.002830,rx[3]=-0.002878,rx[4]=-0.002578,rx[5]=-0.0052,rx[6]=-0.003981, ry[1]=-0.0055,ry[2]=0.007694,ry[3]=0.005098,ry[4]=0.007429,ry[5]=-0.009152,ry[6]=-0.007441, bu=rx[1]-lx[1];

do{ q++; bv=bu*x;bw=bu*y;

double R[8][8],b[8][8],l[8][8],c[8][8],A[8][8],AT[8][8],Z[8][8],ATA[8][8],d[8][8];

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R[0][2]=-sin(Q)*cos(W); R[1][0]=cos(W)*sin(K); R[1][1]=cos(W)*cos(K); R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K);

R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);

14 / 18

R[2][2]=cos(Q)*cos(W);

b[0][0]=lx[1];b[1][0]=ly[1];b[2][0]=-f;b[0][2]=lx[3];b[1][2]=ly[3];b[2][2]=-f;b[0][4]=lx[5];b[1][4]=ly[5];b[2][4]=-f; f;

d[0][0]=rx[1];d[1][0]=ry[1];d[2][0]=-f;d[0][2]=rx[3];d[1][2]=ry[3];d[2][2]=-f;d[0][4]=rx[5];d[1][4]=ry[5];d[2][4]=-f; =-f;

Mult(R,d,c,3,3,6);

{l[j][0]=b[1][j]*(bu*c[2][j]-bw*c[0][j])/(c[2][j]*b[0][j]-c[0][j]*b[2][j])-c[1][j]*(bu*b[2][j]-bw*b[0][j])/(c[2][j]*b A[j][0]=-c[0][j]*c[1][j]*(bu*b[2][j]-bw*b[0][j])/((c[2][j]*b[0][j]-c[0][j]*b[2][j])*c[2][j]); A[j][1]=-(c[2][j]+c[1][j]*c[1][j]/c[2][j])*(bu*b[2][j]-bw*b[0][j])/(c[2][j]*b[0][j]-c[0][j]*b[2][j]); A[j][2]=(bu*b[2][j]-bw*b[0][j])/(c[2][j]*b[0][j]-c[0][j]*b[2][j])*c[0][j]; A[j][3]=bu;

A[j][4]=-bu*c[1][j]/c[2][j]; }

for(k=0;k<=4;k++)

{AT[k][i]=A[i][k];}

for(j=0;j<=5;j++) [0][j]-c[0][j]*b[2][j])-bv;

d[0][1]=rx[2];d[1][1]=ry[2];d[2][1]=-f;d[0][3]=rx[4];d[1][3]=ry[4];d[2][3]=-f;d[0][5]=rx[6];d[1][5]=ry[6];d[2][5]b[0][1]=lx[2];b[1][1]=ly[2];b[2][1]=-f;b[0][3]=lx[4];b[1][3]=ly[4];b[2][3]=-f;b[0][5]=lx[6];b[1][5]=ly[6];b[2][5]=-

for(i=0;i<=5;i++)

Mult(AT,A,ATA,5,6,5); inverse(ATA,5); Mult(ATA,AT,Z,5,5,6); Mult(Z,l,X,5,6,1);

Q+=X[0][0];W+=X[1][0];K+=X[2][0];

x+=X[3][0];y+=X[4][0];}

while(fabs(X[0][0])>=0.00003||fabs(X[1][0])>=0.00003||fabs(X[2][0])>=0.00003||fabs(X[3][0])>=0.00003||fabs(X[4][0])>=0.00003);

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R[0][2]=-sin(Q)*cos(W);

R[1][0]=cos(W)*sin(K); R[1][1]=cos(W)*cos(K); R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K); R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);

15 / 18

cout<<\"迭代次数为：\"<cout<<\"φ为：\"<R[2][2]=cos(Q)*cos(W);

}

void Mult(double VT[1][6], double V[6][1],double VTV[1][1], int m, int t, int p) { double sum; for(int i=0;i{

sum = 0; for(int k=0;k{

sum +=VT[i][k]*V[k][j]; }

VTV[i][j] = sum; } }

} //矩阵相乘子函数//

void Mult(double R[8][8], double b[8][8],double c[8][8], int m, int t, int p) { double sum; for(int i=0;i{

sum = 0; for(int k=0;k{

sum +=R[i][k]*b[k][j]; }

c[i][j] = sum; } }

} //矩阵相乘子函数//

void inverse(double c[N][N],int t) { int i,j,h,k; double p; double q[6][12];

16 / 18

for(i=0;ielse

q[i][j]=0;}

for(h=k=0;k{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=t-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;ifor(i=0;i计算结果：

17 / 18

数据分析：

经检验，数据合格。

心得体会：

参考文献：《摄影测量学》——林君建，苍桂华主编

18 / 18