摄影测量实习报告.docx
《摄影测量实习报告.docx》由会员分享，可在线阅读，更多相关《摄影测量实习报告.docx（15页珍藏版）》请在冰豆网上搜索。
摄影测量实习报告
摄影测量实习报告
——单张影像空间后方交会程序设计
一、实习目的
深入理解单张影像空间后方交会呃原理，体会再有多余观测情况下，用最小二乘平差方法编程实现解求影像外方为元素的过程。
利用VisualC++或者Matlab（或其他熟悉的计算机语言）边学一个完整的单张影像空间后方交会程序，通过对提供的试验数据进行计算，输出相片的外方为元素并进行评定精度。
通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的主动能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。
二、实习环境
1、硬件环境：
Windows操作系统；
2、软件环境：
VC++活Matlab或其他计算机语言。
三、实习内容
利用一定数量的地面控制点，根据共线条件方程求解相片外方为元素并进行精度评定。
四、实习原理
1．共线方程
2．精度评定
其中
五、实习数据
2．模拟像片一对：
左片号23右片号24
2．像片比例尺：
1/30000
3．航摄机主距：
f=150mm
4．每张像片有4个控制点
5．点位略图
3．Δ6Δ12
4．
5．⊙5⊙8
6．
7．Δ4Δ10
8．Δ3Δ9
9．
10．⊙5⊙8
11．
12．Δ1Δ7
6．各片像点坐标及其地面坐标
片
号
点
号
像点坐标（mm）
地面坐标（m）
x
y
X
Y
Z
23
1
-91.596
-74.859
100000.00
137500.00
11.00
3
-94.230
81.446
100000.00
142500.00
36.00
7
95.207
-75.521
106000.00
137500.00
42.00
9
96.797
83.077
106000.00
142500.00
56.00
24
4
-102.695
-79.618
103000.00
137500.00
90.00
6
-99.904
81.754
103000.00
142500.00
31.00
10
86.890
-77.540
109000.00
137500.00
7.00
12
88.904
76.257
109000.00
142500.00
5.00
七、验证数据
1．已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。
2．已知4对点的影像坐标和地面坐标
影像坐标
地面坐标
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
八、程序设计如下
usingSystem;
usingSystem.Collections.Generic;
usingSystem.ComponentModel;
usingSystem.Data;
usingSystem.Drawing;
usingSystem.Text;
usingSystem.Windows.Forms;
usingSystem.IO;
usingSystem.Data.OleDb;
namespace单像空间后方交会
{
publicpartialclass后方交会:
Form
{
public后方交会（）
{
InitializeComponent（）;
}
//定义已知数据，x[]、y[]存储像点坐标，X[]、Y[]、Z[]存储地面点坐标
double[]x=newdouble[4],y=newdouble[4],X=newdouble[4],Y=newdouble[4],Z=newdouble[4];
doublef,d;//f为主距,d为迭代精度
privatevoidopenDataToolStripMenuItem_Click（objectsender,EventArgse）
{
////打开*.txt文件数据
if（ofdOpenFile.ShowDialog（）==DialogResult.OK）
{
StreamReadermystreamreader=newStreamReader（ofdOpenFile.FileName）;//将文件载入到mystreamreader
stringLine;//读取的每行信息将记录在该字符串中
double[]Temp=newdouble[22];//将每一行的信息存放在该一维数组中
inttemp=0;
//将*.text文件中的信息读如到Temp[]一维数组中
while（（Line=mystreamreader.ReadLine（））!
=null）
{
Temp[temp++]=Convert.ToDouble（Line）;
}
mystreamreader.Close（）;
for（inti=0;i<4;i++）
{
//根据所读的数据赋值已知数据
x[i]=Temp[5*i];
y[i]=Temp[5*i+1];
X[i]=Temp[5*i+2];
Y[i]=Temp[5*i+3];
Z[i]=Temp[5*i+4];
//将已知数据添加到列表视图中
ListViewItema;
a=lst已知数据.Items.Add（x[i].ToString（））;
a.SubItems.Add（y[i].ToString（））;
a.SubItems.Add（X[i].ToString（））;
a.SubItems.Add（Y[i].ToString（））;
a.SubItems.Add（Z[i].ToString（））;
}
//将所读数据赋值给f和d，并显示到文本框中
f=Temp[20];
d=Temp[21];
txt主距f.Text=f.ToString（）+"m";
txt限差d.Text=d.ToString（）+"m";
}
}
privatevoidbtnCalculate_Click（objectsender,EventArgse）
{
/////////////////初始计算///////////////////
////已知数据
//double[]x={-0.08615,-0.0534,-0.01478,0.01046};
//double[]y={-0.06899,0.08221,-0.07663,0.06443};
//double[]X={36859.41,37631.08,39100.97,40426.54};
//double[]Y={25273.32,31324.51,24934.98,30319.81};
//double[]Z={2195.17,728.69,2386.50,757.31};
//doublef=0.15324,d=0.000001;
constintc_N=100;//迭代次数
intn=0;//迭代次数统计变量
doubleH,m=0;//H为航高，m为比例尺
////求m的值
intt=0;
doublesum;
double[]mid=newdouble[6];//存储比例尺m的值
for（inti=0;i<4;i++）//求六对线段的长度比
{
for（intj=i+1;j<4;j++）
{
mid[t]=（Math.Sqrt（（X[i]-X[j]）*（X[i]-X[j]）+（Y[i]-Y[j]）*（Y[i]-Y[j]）））/
（Math.Sqrt（（x[i]-x[j]）*（x[i]-x[j]）+（y[i]-y[j]）*（y[i]-y[j]）））;
t++;
}
}
sum=0.0;
for（inti=0;i{
sum=sum+mid[i];
}
m=sum/t;//比例尺为各段长度比的均值
H=f*m;//计算航高
lbl比例尺m.Text="比例尺m="+m.ToString（）;//显示比例尺
lblH.Text="航高H="+H.ToString（）+"m";//显示航高
////定义外方位元素，并附初值
doubleXs,Ys,Zs,φ=0,ω=0,κ=0;
Xs=（X[0]+X[1]+X[2]+X[3]）/4.0;
Ys=（Y[0]+Y[1]+Y[2]+Y[3]）/4.0;
Zs=（Z[0]+Z[1]+Z[2]+Z[3]）/4.0+m*f;
////定义x,y近似值
double[]_x=newdouble[4];
double[]_y=newdouble[4];
////定义共线方程中的分子分母项，便于计算
double[]_X=newdouble[4];
double[]_Y=newdouble[4];
double[]_Z=newdouble[4];
////定义旋转矩阵R
double[,]R=newdouble[3,3];
////定义解方程所用矩阵
double[,]A=newdouble[8,6];//A系数阵
double[,]AT=newdouble[6,8];//A系数阵转置
double[,]ATA=newdouble[6,6];//A的转置与A的乘积
double[,]temp=newdouble[6,12];//临时矩阵
double[,]ATAR=newdouble[6,6];//A的转置与A的乘积的逆
double[,]ATARAT=newdouble[6,8];//A的转置与A的乘积的逆的转置
double[,]L=newdouble[8,1];//误差方程中的常数项
double[,]XX=newdouble[6,1];//X向量
//////////////////开始迭代///////////////////
do
{
////计算旋转矩阵
R[0,0]=Math.Cos（φ）*Math.Cos（κ）-Math.Sin（φ）*Math.Sin（ω）*Math.Sin（κ）;//a1
R[0,1]=-Math.Cos（φ）*Math.Sin（κ）-Math.Sin（φ）*Math.Sin（ω）*Math.Cos（κ）;//a2
R[0,2]=-Math.Sin（φ）*Math.Cos（ω）;//a3
R[1,0]=Math.Cos（ω）*Math.Sin（κ）;//b1
R[1,1]=Math.Cos（ω）*Math.Cos（κ）;//b2
R[1,2]=-Math.Sin（ω）;//b3
R[2,0]=Math.Sin（φ）*Math.Cos（κ）+Math.Cos（φ）*Math.Sin（ω）*Math.Sin（κ）;//c1
R[2,1]=-Math.Sin（φ）*Math.Sin（κ）+Math.Cos（φ）*Math.Sin（ω）*Math.Cos（κ）;//c2
R[2,2]=Math.Cos（φ）*Math.Cos（ω）;//c3
for（inti=0;i<4;i++）
{
//用共线方程计算x，y的近似值
_X[i]=R[0,0]*（X[i]-Xs）+R[1,0]*（Y[i]-Ys）+R[2,0]*（Z[i]-Zs）;
_Y[i]=R[0,1]*（X[i]-Xs）+R[1,1]*（Y[i]-Ys）+R[2,1]*（Z[i]-Zs）;
_Z[i]=R[0,2]*（X[i]-Xs）+R[1,2]*（Y[i]-Ys）+R[2,2]*（Z[i]-Zs）;
_x[i]=-f*_X[i]/_Z[i];
_y[i]=-f*_Y[i]/_Z[i];
}
for（inti=0;i<4;i++）
{
//计算系数矩阵
A[2*i,0]=（R[0,0]*f+R[0,2]*x[i]）/_Z[i];
A[2*i,1]=（R[1,0]*f+R[1,2]*x[i]）/_Z[i];
A[2*i,2]=（R[2,0]*f+R[2,2]*x[i]）/_Z[i];
A[2*i,3]=y[i]*Math.Sin（ω）-（（x[i]/f）*（x[i]*Math.Cos（κ）-y[i]*Math.Sin（κ））+f*Math.Cos（κ））*Math.Cos（ω）;
A[2*i,4]=-f*Math.Sin（κ）-（x[i]/f）*（x[i]*Math.Sin（κ）+y[i]*Math.Cos（κ））;
A[2*i,5]=y[i];
A[2*i+1,0]=（R[0,1]*f+R[0,2]*y[i]）/_Z[i];
A[2*i+1,1]=（R[1,1]*f+R[1,2]*y[i]）/_Z[i];
A[2*i+1,2]=（R[2,1]*f+R[2,2]*y[i]）/_Z[i];
A[2*i+1,3]=-x[i]*Math.Sin（ω）-（（x[i]/f）*（x[i]*Math.Cos（κ）-y[i]*Math.Sin（κ））-f*Math.Sin（κ））*Math.Cos（ω）;
A[2*i+1,4]=-f*Math.Cos（κ）-（y[i]/f）*（x[i]*Math.Sin（κ）+y[i]*Math.Cos（κ））;
A[2*i+1,5]=-x[i];
//计算常数项
L[2*i,0]=x[i]-_x[i];
L[2*i+1,0]=y[i]-_y[i];
}
////计算A的转置，存在AT中
for（inti=0;i<6;i++）
for（intj=0;j<8;j++）
{
AT[i,j]=A[j,i];
}
////计算A的转置与A的乘积,存在ATA中
for（inti=0;i<6;i++）
for（intj=0;j<6;j++）
{
ATA[i,j]=0;
for（intl=0;l<8;l++）
ATA[i,j]+=AT[i,l]*A[l,j];
}
////////计算ATA的逆矩阵////////
////初始化temp--临时数组
for（intl=0;l<6;l++）
for（intp=0;p<12;p++）
{
temp[l,p]=0;
}
////把ATA各值赋给temp
for（inti=0;i<6;i++）
for（intj=0;j<6;j++）
{
temp[i,j]=ATA[i,j];
}
////在temp中加入初等方阵
for（intl=0;l<6;l++）
temp[l,l+6]=1;
////初等变换
for（intl=0;l<6;l++）
{
if（temp[l,l]!
=1）
{
doublebs=temp[l,l];
temp[l,l]=1;
for（intp=l+1;p<12;p++）
temp[l,p]/=bs;
}
for（intq=0;q<6;q++）
{
if（q!
=l）
{
doublebs=temp[q,l];
for（intp=l;p<12;p++）
temp[q,p]-=bs*temp[l,p];
}
else
continue;
}
}
////得到ATA的逆阵后存在ATAR中
for（inti=0;i<6;i++）
for（intj=0;j<6;j++）
{
ATAR[i,j]=temp[i,j+6];
}
////ATAR*AT存在ATARAT中
for（inti=0;i<6;i++）
for（intj=0;j<8;j++）
{
ATARAT[i,j]=0;
for（intl=0;l<6;l++）
ATARAT[i,j]+=ATAR[i,l]*AT[l,j];
}
////计算ATARAT*L，存在XX中
for（inti=0;i<6;i++）
for（intj=0;j<1;j++）
{
XX[i,j]=0;
for（intl=0;l<8;l++）
XX[i,j]+=ATARAT[i,l]*L[l,0];
}
n++;
if（n>c_N）
{
MessageBox.Show（"请检查阈值和迭代次数！
"+（n-1）,"计算失败",MessageBoxButtons.OK,MessageBoxIcon.Warning）;
break;
}
////计算外方位元素值
Xs+=XX[0,0];
Ys+=XX[1,0];
Zs+=XX[2,0];
φ+=XX[3,0];
ω+=XX[4,0];
κ+=XX[5,0];
}
while（Math.Abs（XX[0,0]）>=d||Math.Abs（XX[1,0]）>=d||Math.Abs（XX[2,0]）>=d||Math.Abs（XX[3,0]）>=1000*d||Math.Abs（XX[4,0]）>=1000*d||Math.Abs（XX[5,0]）>=1000*d）;
////输出显示结果
if（n>c_N）
MessageBox.Show（"超过迭代精度，此数值下可能不收敛！
","提示",MessageBoxButtons.OK,MessageBoxIcon.Information）;
else
{
MessageBox.Show（"计算完成，显示计算结果","提示",MessageBoxButtons.OK,MessageBoxIcon.Information）;
txt计算结果.Text="Xs:
"+Xs.ToString（）+"m"+"\r\n"+"Ys:
"+Ys.ToString（）+"m"+"\r\n"+"Zs:
"+Zs.ToString（）+"m"+"\r\n\r\n"
+"φ:
"+φ.ToString（）+"rad"+"\r\n"+"ω:
"+ω.ToString（）+"rad"+"\r\n"+"κ:
"+κ.ToString（）+"rad"+"\r\n\r\n"
+"dXs:
"+XX[0,0].ToString（）+"m"+"\r\n"+"dYs:
"+XX[1,0].ToString（）+"m"+"\r\n"+"dZs:
"+XX[2,0].ToString（）+"m"+"\r\n"
+"dφ:
"+XX[3,0].ToString（）+"rad"+"\r\n"+"dω:
"+XX[4,0].ToString（）+"rad"+"\r\n"+"dκ:
"+XX[5,0].ToString（）+"rad"+"\r\n\r\n"+"迭代次数：
"+n.ToString（）;
}
////保存结果
if（MessageBox.Show（"是否保存结果数据.","提示",MessageBoxButtons.YesNo,MessageBoxIcon.Question）==DialogResult.Yes）
{
if（ofdSaveFile.ShowDialog（）==DialogResult.OK）
{
StreamWritermystreamwriter=newStreamWriter（ofdSaveFile.FileName）;
mystreamwriter.WriteLine（"计算结果为：
"）;
mystreamwriter.WriteLine（"m:
"+m.ToString（））;
mystreamwriter.WriteLine（"H:
"+H.ToString（）+"m"）;
mystreamwriter.WriteLine（"\r\n"）;
升级会员 