空间后方交会

实习报告

学号 2010302590232

姓名 李雷

指导老师

实习目的

1、掌握空间后方交会的定义和实现算法

(1) 定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。

(2) 算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。

2、了解摄影测量平差的基本过程

(1) 获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。

(2) 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。

(3) 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。

(4) 计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。

(5) 逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。

(6) 逐点计算误差方程式的系数和常数项，组成误差方程式。

(7) 计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。

(8) 解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。

(9) 检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.1′，当3个改正数均小于0.1′时，迭代结束。否则用新的近似值重复(4)～（8）步骤的计算，直到满足要求为止。

(10)流程图

一、程序源代码

typedef struct data{

int no;

double x;

double y;

double X;

double Y;

double Z;

}data;

#include <iostream.h>

#include<math.h>

#include <iomanip.h>

#include "fstream"

#include "HEADFile.h"

double** zhuan(double** a,int n,int m) //矩阵转置

{

double** p=new double*[m];

for(int i=0;i!=m;i++)

p[i]=new double[n];

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

for(int j=0;j<n;j++)

{

p[i][j]=a[j][i];

}

return p;

}

double **MatrixMul(double **m1,int i1,int j1,double **m2,int i2,int j2) //矩阵相乘（矩阵1，行数，列数，矩阵2，行数，列数），函数内已开僻内存，直接返回给定义指针

{

if(i2!=j1)

{

cout<<"error!两矩阵不满足相乘条件！"<<endl;

return NULL;

}

int i,j,k;

double **p;

p=new double*[i1];

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

p[i]=new double[j2];

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

for(j=0;j<j2;j++)

{

p[i][j]=0.0;

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

p[i][j]+=m1[i][k]*m2[k][j];

}

return p;

}

double **InverseMatrix(double** p,int n) //矩阵求逆函数

{

int i,j,k,row,col;

double temp,temp0,kb,kc,kd;

double **e=new double*[n];

for (i=0;i<n;i++) //构造单位矩阵

{

e[i]=new double[n];

for (j=0;j<n;j++)

{

if(i==j)

e[i][j]=1;

else

e[i][j]=0;

}

}

for(i=0;i<n;i++) //逐个找对角线上的元素，并把该元素所在列上的本行以下元素化为0

{

int i0=i;

while(i0<n-1 && p[i0][i]==0 ) //找出某列上第一个不为零的行数

{

i0++;

}

if (p[i0][i]==0)

{

cout<<"The square have no inverse matrix!"<<endl;

break;

}

for (j=0;j<n;j++) //当前行与该行交换

{

temp=p[i0][j];

p[i0][j]=p[i][j];

p[i][j]=temp;

temp0=e[i0][j];

e[i0][j]=e[i][j];

e[i][j]=temp0;

}

kc=p[i][i];

for (j=0;j<n;j++) //交换后的当前行首元素p[i][i]化为1 {

p[i][j]/=kc;

e[i][j]/=kc;

}

if(i<n-1) //若当前行不是矩阵的最后一行，则使当前行以下各行减去本行首元素与当前行的积

{

for(row=i+1;row<n;row++)

{

kb=p[row][i];

for (col=0;col<n;col++)

{

p[row][col]-=kb*(p[i][col]);

e[row][col]-=kb*(e[i][col]);

}

}

}

}

for (i=n-2;i>=0;i--)

{

for(k=i+1;k<n;k++)

{

kd=p[i][k];

for (j=0;j<n;j++)

{

p[i][j]-=kd*p[k][j];

e[i][j]-=kd*e[k][j];

}

}

}

return e;

}

void main()

{

int M,num,i,j;

FILE *fp;

double f;

fp=fopen("E:\\后方交会数据.txt","r");

fscanf(fp,"%d",&num);

data *a=new data[num];

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

{

fscanf(fp,"%d%lf%lf%lf%lf%lf",&a[i].no,&a[i].x,&a[i].y,&a[i].X,&a[i].Y,&a[i].Z);

}

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

{

a[i].x/=1000;

a[i].y/=1000;

}

fscanf(fp,"%lf%d",&f,&M);

f/=1000;

fclose(fp); //读数据

double

phi=0,omega=0,kappa=0,Xs=38437,Ys=27963.155,Zs=7662,Xb=0,Yb=0,Zb=0;

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

{

Xs+=a[i].X; Ys+=a[i].Y;

Zs+=a[i].Z;

}

Xs/=num;

Ys/=num; //计算估计值

Zs=Zs/num+f*M;

double cor1=1;

int Num=2*num;

double** A=new double*[Num]; //定义系数矩阵

for(i=0;i!=Num;i++)

A[i]=new double[6];

double** L=new double*[Num]; //定义误差方程常数量矩阵

for(i=0;i!=Num;i++)

L[i]=new double[1];

double a1,a2,a3,b1,b2,b3,c1,c2,c3;

double** p,**q,**r,**s,**t; //定义下面矩阵运算的数组名

int ab=0;

while(cor1>1E-8) //迭代条件

{

a1=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);

a2=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);

a3=-sin(phi)*cos(omega);

b1=cos(omega)*sin(kappa);

b2=cos(omega)*cos(kappa);

b3=-sin(omega);

c1=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);

c2=-sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa);

c3=cos(phi)*cos(omega);

for(i=0,j=0;i<Num;i+=2,j++) /*给系数矩阵赋值

*/

{

Xb=a1*(a[j].X-Xs)+b1*(a[j].Y-Ys)+c1*(a[j].Z-Zs);

Yb=a2*(a[j].X-Xs)+b2*(a[j].Y-Ys)+c2*(a[j].Z-Zs);

Zb=a3*(a[j].X-Xs)+b3*(a[j].Y-Ys)+c3*(a[j].Z-Zs);

A[i][0]=(a1*f+a3*a[j].x)/Zb;

A[i][1]=(b1*f+b3*a[j].x)/Zb; A[i][2]=(c1*f+c3*a[j].x)/Zb;

A[i][3]=a[j].y*sin(omega)-((a[j].x*(a[j].x*cos(kappa)-a[j].y*sin(kappa))/f)+f*cos(kappa))*cos(o

mega);

A[i][4]=-f*sin(kappa)-a[j].x*(a[j].x*sin(kappa)+a[j].y*cos(kappa))/f;

A[i][5]=a[j].y;

A[i+1][0]=(a2*f+a3*a[j].y)/Zb;

A[i+1][1]=(b2*f+b3*a[j].y)/Zb;

A[i+1][2]=(c2*f+c3*a[j].y)/Zb;

A[i+1][3]=-a[j].x*sin(omega)-(a[j].y*(a[j].x*cos(kappa)-a[j].y*sin(kappa))/f-f*sin(kappa))*cos(o

mega);

A[i+1][4]=-f*cos(kappa)-a[j].y*(a[j].x*sin(kappa)+a[j].y*cos(kappa))/f;

A[i+1][5]=-a[j].x;

L[i][0]=a[j].x+f*Xb/Zb;

L[i+1][0]=a[j].y+f*Yb/Zb;

}

p=zhuan(A,Num,6); //矩阵转置

q=MatrixMul(p,6,Num,A,Num,6); //矩阵相乘

r=InverseMatrix(q,6); //矩阵求逆，其结果同时是下面要

求的Nbb矩阵的逆，即参数值的协方差矩阵

s=MatrixMul(r,6,6,p,6,Num);

t=MatrixMul(s,6,Num,L,Num,1);

Xs+=t[0][0];

Ys+=t[1][0];

Zs+=t[2][0];

phi+=t[3][0];

omega+=t[4][0];

kappa+=t[5][0];

cor1=fabs(t[3][0]);

ab++;

if (ab>30)

{

cout<<"The data may be incorrect!"<<endl;

exit(1);

}

}

cout<<"迭代次数："<<ab<<endl; //查看迭代次数 double **V,sigama2;

double **m1,**m2; //定义矩阵运算的中间变量

V=MatrixMul(A,Num,6,t,6,1);

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

{

V[i][0]-=L[i][0];

}

m1=zhuan(V,Num,1);

m2=MatrixMul(m1,1,Num,V,Num,1);

sigama2=m2[0][0]/(Num-6); //单位权方差

double d[6]; //定义六个参数的中误差的数组

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

{

d[i]=sqrt(r[i][i]*sigama2);

}

cout<<setprecision(7)<<"平差结果及其中误差 :"<<endl<<"Xs="<<Xs<<" 中误差： "<<d[0]<<endl<<"Ys="<<Ys<<" 中误差：

"<<d[1]<<endl<<"Zs="<<Zs<<" 中误差："<<d[2]<<endl;

cout<<"phi="<<phi<<" 中误差："<<d[3]<<endl<<"omega="<<omega<<" 中误差："<<d[4]<<endl<<"kappa="<<kappa<<" 中误差："<<d[5]<<endl;

}

二、实习数据

三、;计算结果

四、实习心得

编写程序前应该先把实验的原理弄清楚，编写程序时要细致小心，同时，在编写每一小段程序时都要有注释，以便运行出错时方便找出错误并加以修改。程序运行正确后，再把程序进行修饰，把不必要的注释去掉，增加代码的美观性。通过实习，让我对空间后方交会的原理与方法有了更深一步的理解，并且在此过程中锻炼了自己的编程能力。通过实习也让我发现了自己的不足，知道自己有很多东西没有实际的掌握。我会不断的努力学习，使自己的知识掌握得更加牢固。