摄影测量实习报告
——单张影像空间后方交会程序设计
实习时间 2013.5.20-2013.5.24
学生班级 测绘10-2班
学生姓名 刘航
学生学号 1072143212
所在院系 矿业工程学院
指导教师 张会战 邵亚琴
一、实习目的
1．深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
2．掌握空间后方交会的定义和实现算法
(1) 定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。
(2) 算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。
3. 利用Visual C++或者Matlab（或其他熟悉的计算机语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并进行评定精度。
4．通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的动手能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。
二、实习环境
1．硬件环境：Window操作系统
2．软件环境：VC++或Matlab或其他计算机语言
三、实习内容
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素并进行精度评定。
四、实习原理
以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素 Xs，Ys，Zs，ф，ω，κ。共线条件方程如下：
其中：x,y为像点的像平面坐标；x0,y0,f为影像的外方位元素；
Xs,Ys,Zs为摄站点的物方空间坐标；X,Y,Z为物方点的物方空间坐标；
2．精度评定
其中
程序框图
五、实习数据
1．点位略图
2．模拟像片一对：左片号23
3．像片比例尺： 1/30000
4．航摄机主距：f=150mm
5．每张像片有4个控制点
6．各片像点坐标及其地面坐标
六、验证数据
1．已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。
2．已知4对点的影像坐标和地面坐标
3．解算结果
4.数据分析
选取第六张像片进行计算，迭代次数为2次。经过比较发现，计算出的6个外方位元素与所给参考值相比，相差很小，计算结果符合要求：线元素误差小于0.5米；角元素误差30秒。计算其精度，可以通过法方程式中未知数的系数矩阵的逆阵A-1来求解，此时，视像点坐标为等精度不相关观测值。线元素精度mx等，高于0.05米；角元素精度高于0.00003弧度。计算结果都达到标准。在此次计算中，我运用了所给的全部控制点，而空间后方交会所运用的控制点，应该避免位于一个圆柱面上，否则会出现解不唯一的情况。选点时，还需要避免选择的点过于聚集在一起，或位于一条直线上，所选控制点最好分布在像片的四角和中央。并且数量充足，这样有利于提高解算精度。迭代时，所选择控制条件不同，迭代次数略有不同，所以最后结果也会略有不同。一般设置为线元素改正数小于0.01m，角元素改正数小于0.1´。
所提供X Y Z为地面测量坐标，带入共线方程时，需要转换为地面摄影测量坐标，最简单的方法为互换XY的数值，即可达到转换坐标目的。并且其单位为米，而像点坐标的单位为厘米，需要统一坐标单位。
七、程序
输入文件形式如下：
C++源程序如下：
#include <iostream>
#include <fstream>
#include <cmath>
#include <string>
#include <iomanip>
using namespace std;
const int n=6;
void inverse (double c[n][n]);
template<typename T1,typename T2>void transpose (T1*mat1,T2*mat2,int a,int b);
template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);
template<typename T>void input (T*mat,int a,int b);
template<typename T>void output(T*mat,char*s,int a,int b);
int main()
{
ofstream outFile;
cout.precision(5);
double x0=0.0, y0=0.0; double fk=0.15324; //内方位元素
double m=39689; //估算比例尺
double B[4][5]={0.0},R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];
input (B,4,5); //从文件中读取控制点的影像坐标和地面坐标，存入数组B
double Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0;
double X,Y,Z,L[8][1],A[8][6];
//确定未知数的出始值
for(int i=0;i<4;i++)
{Xs=Xs+B[i][2];
Ys=Ys+B[i][3];
Zs=Zs+B[i][4];
}
Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk;
int f=0;
do//迭代计算
{f++;
//组成旋转矩阵
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);
R[2][2]=cos(Q)*cos(W);
//计算系数阵和常数项
for(int i=0,k=0,j=0;i<=3;i++,k++,j++)
{
X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs);
Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs);
Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs);
L[j][0]=B[i][0]-(x0-fk*X/Z);
L[j+1][0]=B[i][1]-(y0-fk*Y/Z);
j++;
A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z;
A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z;
A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;
A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk+fk*cos(K))*cos(W);
A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;
A[k][5]=B[i][1]-y0;
A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z;
A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z;
A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;
A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][1]-y0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin(K))*cos(W);
A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;
A[k+1][5]=-(B[i][0]-x0);
k++;
}
transpose(A,AT,6,8);
multi(AT,A,ATA,6,8,6);
inverse(ATA);
multi(AT,L,ATL,6,8,1);
multi(ATA,ATL,XG,6,6,1);
Xs=Xs+XG[0][0]; Ys=Ys+XG[1][0]; Zs=Zs+XG[2][0];
Q=Q+XG[3][0]; W=W+XG[4][0]; K=K+XG[5][0];
}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265.0);
cout<<"迭代次数为："<<f<<endl;
//精度评定
double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];
multi(A,XG,AXG,8,6,1);
for( i=0;i<8;i++) //计算改正数
V[i][0]=AXG[i][0]-L[i][0];
transpose (V,VT,1,8);
multi(VT,V,VTV,1,8,1);
m0=VTV[0][0]/2;
for(i=0;i<6;i++)
for(int j=0;j<6;j++)
D[i][j]=m0*ATA[i][j];
//屏幕输出误差方程系数阵、常数项、改正数
output(A,"误差方程系数阵A为：",8,6);
output(L,"常数项L为：",8,1);
output(XG,"改正数为：",6,1);
outFile.open("aim.txt",ios::app); //打开并添加aim.txt文件
outFile.precision(10);
//以文件的形式输出像片外方位元素、旋转矩阵、方差阵
outFile<<"一、像片的外方位元素为："<<endl<<endl;
outFile<<setw(10)<<"Xs="<<Xs<<setw(10)<<"Ys="<<Ys<<setw(10)<<"Zs="
<<Zs<<endl;
outFile<<setw(20)<<"航向倾角为："<<Q<<setw(10)<<
"旁向倾角为："<<W<<setw(10)<<"像片旋角为："<<K<<endl;
outFile<<endl<<"二、旋转矩阵R为:"<<endl<<endl;
for( i=0;i<3;i++)
{for(int j=0;j<3;j++)
outFile<<setw(25)<<R[i][j]<<setw(25);
outFile<<endl;
}
outFile<<endl;
outFile<<setw(0)<<"三、精度评定结果为："<<endl;
outFile.precision(5);
for(i=0;i<6;i++)
{for(int j=0;j<6;j++)
outFile<<setw(14)<<D[i][j]<<setw(14);
outFile<<endl;
}
outFile.close();
return 0;
}
template<typename T1,typename T2>void transpose(T1*mat1,T2*mat2,int a,int b)
{ int i,j;
for(i=0;i<b;i++)
for(j=0;j<a;j++)
mat2[j][i]=mat1[i][j];
return;
}
template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c)
{ int i,j,k;
for(i=0;i<a;i++)
{for(j=0;j<c;j++)
{result[i][j]=0;
for(k=0;k<b;k++)
result[i][j]+=mat1[i][k]*mat2[k][j];
}
}
return;
}
template <typename T>void input (T*mat,int a,int b)
{ ifstream inFile;
inFile.open("控制点坐标.txt");
while(!inFile.eof())
{for (int i=0;i<a;i++)
for(int j=0;j<b;j++)
inFile>>mat[i][j];
}
inFile.close();
return;
}
template<typename T>void output(T*mat,char*s,int a,int b)
{ cout<<setw(15)<<s<<endl;
for(int i=0;i<a;i++)
{for(int j=0;j<b;j++)
cout<<setw(13)<<mat[i][j];
cout<<endl;
}
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];
}
程序的结果输出如下：（包括文本输出结果和荧屏输出中间数据）
八、实习心得
通过此次的上机实习，我更进一步的了解了单片空间后方交会的原理与计算步骤，不仅加深了对课堂上所听到的理论知识的理解与认识，并将理论给予实践的考验，理论联系实践，加深了对其的认识。同时还增强了自己动手编程和程序调试的能力，使自己在程序设计和编辑上更加熟悉。这次实习持续时间很长，经历了几次数据的更改，所以程序也几经修改，由最初的直接输入数据到后来可以自行读入数据，并且可以选择计算的像片，功能有所完善，我也在实践的过程中，对空间后方交会有了更深的理解。深刻理解了共线条件方程的运用，各个量的意义，受益匪浅。
虽然在实习的初期对单片空间后方交会的理论知识理解的不是很好，但在实习过程中一步步摸索前进，最后对后方交会的计算过程整体的来龙去脉有了一个大致的认识。在此次实习中也碰到过不少的关于部分计算程序实现的难题，但在老师的指导和同学的帮助下还是完成在整个实习，在此表示真心的感谢。