摄影测量学后前交会实验报告

江西理工大学应用科学学院

摄 影 测 量 实 验 报 告

班级：测绘101班 姓名：邬志刚

学号：08090110125

一、实验目的：

掌握空间后方交会和前方交会的原理，根据所给控制点的地面物方坐标以及相应的像点在像平面坐标系中的坐标和同名像点在左右像片上的坐标，利用计算机编程语言实现空间后方交会、前方交会的过程，完成所给立体像对中两张像片各自的6个外方位元素的解求和所给立体像对上5对同名点对应的地面物方点的坐标计算。

二、实验原理：

如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标

依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。

那么由共线条件方程可知：

a(X Xs)) b(Y Ys)) cc(Z Zs))1(X X1(Y Y1(Z Za b 1s1s1sx x0 f

x x0 f(X X) b(Y Y) c(Z Z)

3(X Xs) b3(Y Ys) c3(Z Zs)a

3s3s3s

a(X Xs)) b(Y Ys)) cc(Z Zs))2(X X2(Y Y2(Z Za b 2s2s2syy 0 f yy0 fa(X Xs)) b(Y Ys)) cc(Z Zs))3(X X3(Y Y3(Z Z a b 3s3s3s

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

在方程中共有6个未知参数Xs,Ys,Zs, , , ,所以有三个不在一条直

线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级

数展开，取一次项得到线型表达式

，

如下：

x (x)

d d d dXs dYs dZs Xs Ys Zs

y y y y y y

y (y) d d d dXs dYs dZs

Xs Ys Zs

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

x x x x x x

x (x) a11dXs a12dYs a13dZs a14d a15d a16d

y (y) a21dXs a22dYs a23dZs a24d a25d a26d

其中（x）、（y）是函数近似值，d ,d ,d ,dXs,dYs,dZs是外方位元素

近似值的改正数，它们的系数为函数的偏导数。为了便于推导，令：

X=a1(XY=a2(XZ=a3(X

Xs) b1(Y Ys) c1(Z Zs)

Xs) b2(Y Ys) c2(Z Zs) Xs) b3(Y Ys) c3(Z Zs)

那么有

X Y Z a1 a2 a3

b1b2b3

c1 X Xs X Xs T c2Y Ys RY Ys

c3 Z Zs Z Zs

x Xs

对于系数，其严密算法（以

,

x

为例）如下：

x X

s

1

fZ

2

(

X X

s

Z X

Z X

s

X) 1

fZ

2

( a1Z a3X)

a11

a

21于直影言，

a14a24

a1f a3(x x0) (a1f fa3)

ZZZ

fx ,a12 0,a13

HH

fy

0,a22 ,a23

HH 2

xxy f(1 2),a15 ,a16 y

ff

2 xyy

,a25 f(1 2),a26 x

ff

对竖摄而像

片的角方位元素都是小值，因而各系数的近似值为：

x

fZ

(

X

X ZZ

)

f X

(b1Y b2X) b2Z b3Y

Z Z b2f b3f

YZ

f

XZ

(b1

YZ

b2

XZ

)

fcos cos (y y0)sin (x x0) (y y0) (x x0)

cos sin cos cos

ff (x x0)

(y y0)sin

f

(x x0)cos

fcos cos (y y0)sin

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

vx a11dXvy a21dX

ss

a12dYs a13dZs a14d a15d a16d lx a22dYs a23dZs a24d a25d a26d ly

用矩阵表示为：

具体的

a11

A

a21 vx V ,

vy

Xs Ys Zs x ,

a12a22

a13a23

a14a24

x0 x

l 0

y y

流程图框：

a15a25

a16 a26

四、实验源程序：

#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;

int 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]=5083.205;X[1]=5780.02;X[2]=5210.879;X[3]=5909.264;

Y[0]=5852.099;Y[1]=5906.365;Y[2]=4258.446;Y[3]=4314.283;

Z[0]=527.925;Z[1]=571.549;Z[2]=461.81;Z[3]=455.484;

/*---------------输入点像片坐标---------------------*/ x[0]=16.012;x[1]=88.56;x[2]=13.362;x[3]=82,24; y[0]=79.963;y[1]=81.134;y[2]=-79.37;y[3]=-80.027;

cout<<endl;

/*-------------误差方程中各偏导数的值--------------*/ 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];

/*------------------迭代计算--------------------------*/

while((XX[3]>6/206265 || XX[4]>6/206265 || XX[5]>6/206265)&&n<100) {

/*-----------------设定外方位元素初始值--------------*/ x0=0;y0=0;f=152;m=10000; Xs=0;Ys=0;Zs=f*m/10000;

q=0;w=0;k=0; XX[3]=1;

/*----------------旋转矩阵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];

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;

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]; L[2*i]=x[i]-x1[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; }

五、实验计算结果：

六、实验结果分析：

由计算结果可知，最后的计算结果是符合 6秒的限差要求的。由

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

。

七、实验感想：

写内容：（字数要求600字左右）

上课的时候，老师一讲好像什么都明白，但是每当合上课本的时候，好像什么都忘记了，为了能有个比较完整的思路，我一边看书上求解的方法，一边把求解的过程用流程图在草稿纸上粗略的表达了一下，这对我后面的编程有极大的帮助。在对课本中知识有了一个系统完善的理解以后，接下来就是一个艰巨的过程，

对于理论知识我们多看几遍就可以有个粗略的理

解，但对于实践来说，首先我们必须对理论有个透彻的理解，然后要将它转化到实际应用中去。在刚刚编写代码的时候，只是慢慢的不断地摸索，对于一个语句的定义，一个结构的格式我都只能到书本和网络上一步一步的求解。这在刚刚开始的时候对我来说是个极大的考验。经过几天的临时突击，对于编程的语言我已经有了一个大概的了解，在设计整个程序所需要的，我基本已经掌握了。接下来就是具体功能编写了。根据前面所总结的流程图，我对整个程序的框架有了一个比较完善的把握。整个大的程序中，我将它们割成几个小的功能，对于每个功能，一步一步来实现。

。 注意：1、实验要求在“双向解析摄影测量”（1、空间后方交会——

前方交会；2、解析法相对定向——绝对定向；3、双像光束法整体结算）中任意选择一种进行编程计算；

2、所有同学运行完成后保留源代码； 3、所有同学必须按照以上格式进行书写；

4、要求主标题用2号黑体、副标题用4号宋体、正文用小4

号宋体，行间距24磅；

5、要求用A4纸单面打印，装订成册；

6、上交时间：本学期第15周星期三，统一上交学习委员，按

学号排好；