摄影测量实习报告
——单张影像空间后方交会程序设计
一、 实习目的
深入理解单张影像空间后方交会呃原理，体会再有多余观测情况下，用最小二乘平差方法编程实现解求影像外方为元素的过程。
利用Visual C++或者Matlab（或其他熟悉的计算机语言）边学一个完整的单张影像空间后方交会程序，通过对提供的试验数据进行计算，输出相片的外方为元素并进行评定精度。
通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的主动能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。
二、 实习环境
1、 硬件环境：Windows 操作系统；
2、 软件环境：VC++活Matlab 或其他计算机语言。
三、 实习内容
利用一定数量的地面控制点，根据共线条件方程求解相片外方为元素并进行精度评定。
四、 实习原理
1．共线方程
2．精度评定
其中
五、实习数据
2．模拟像片一对：左片号23 右片号24
2．像片比例尺： 1/30000
3．航摄机主距：f=150mm
4．每张像片有4个控制点
5．点位略图
6．各片像点坐标及其地面坐标
七、验证数据
1．已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。
2．已知4对点的影像坐标和地面坐标
八、程序设计如下
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
using System.IO;
using System.Data.OleDb;
namespace 单像空间后方交会
{
public partial class 后方交会 : Form
{
public 后方交会()
{
InitializeComponent();
}
//定义已知数据，x[]、y[]存储像点坐标，X[]、Y[]、Z[]存储地面点坐标
double []x = new double[4], y = new double[4],X = new double[4], Y = new double[4], Z = new double[4];
double f, d;//f为主距,d为迭代精度
private void openDataToolStripMenuItem_Click(object sender, EventArgs e)
{
////打开*.txt文件数据
if (ofdOpenFile.ShowDialog() == DialogResult.OK)
{
StreamReader mystreamreader = new StreamReader(ofdOpenFile.FileName);//将文件载入到mystreamreader
string Line;//读取的每行信息将记录在该字符串中
double[] Temp = new double[22];//将每一行的信息存放在该一维数组中
int temp = 0;
//将*.text文件中的信息读如到Temp[]一维数组中
while ((Line = mystreamreader.ReadLine()) != null)
{
Temp[temp++] = Convert.ToDouble(Line);
}
mystreamreader.Close();
for (int i = 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];
//将已知数据添加到列表视图中
ListViewItem a;
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";
}
}
private void btnCalculate_Click(object sender, EventArgs e)
{
/////////////////初始计算///////////////////
////已知数据
//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 };
//double f = 0.15324, d = 0.000001;
const int c_N = 100;//迭代次数
int n=0;//迭代次数统计变量
double H,m=0;//H为航高，m为比例尺
//// 求m的值
int t = 0;
double sum;
double[] mid = new double[6];//存储比例尺m的值
for (int i = 0; i < 4; i++)//求六对线段的长度比
{
for (int j = 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 (int i = 0; i < t; i++)
{
sum = sum + mid[i];
}
m = sum / t;//比例尺为各段长度比的均值
H = f * m;//计算航高
lbl比例尺m.Text = "比例尺m=" + m.ToString();//显示比例尺
lblH.Text = "航高H=" + H.ToString()+"m";//显示航高
////定义外方位元素，并附初值
double Xs , 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 = new double[4];
double[] _y = new double[4];
////定义共线方程中的分子分母项，便于计算
double[] _X = new double[4];
double[] _Y = new double[4];
double[] _Z = new double[4];
////定义旋转矩阵R
double[,] R = new double[3, 3];
////定义解方程所用矩阵
double[,] A = new double[8, 6];//A系数阵
double[,] AT = new double[6, 8];//A系数阵转置
double[,] ATA = new double[6, 6];//A的转置与A的乘积
double[,] temp = new double[6, 12];//临时矩阵
double[,] ATAR = new double[6, 6];//A的转置与A的乘积的逆
double[,] ATARAT = new double[6, 8];//A的转置与A的乘积的逆的转置
double[,] L = new double[8, 1];//误差方程中的常数项
double[,] XX = new double[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 (int i = 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 (int i = 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 (int i = 0; i < 6; i++)
for (int j = 0; j < 8; j++)
{
AT[i, j] = A[j, i];
}
////计算A的转置与A的乘积,存在ATA中
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
{
ATA[i, j] = 0;
for (int l = 0; l < 8; l++)
ATA[i, j] += AT[i, l] * A[l, j];
}
////////计算ATA的逆矩阵////////
////初始化temp--临时数组
for (int l = 0; l < 6; l++)
for (int p = 0; p < 12; p++)
{
temp[l, p] = 0;
}
////把ATA各值赋给temp
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
{
temp[i, j] = ATA[i, j];
}
////在temp中加入初等方阵
for (int l = 0; l < 6; l++)
temp[l, l + 6] = 1;
////初等变换
for (int l = 0; l < 6; l++)
{
if (temp[l, l] != 1)
{
double bs = temp[l, l];
temp[l, l] = 1;
for (int p = l + 1; p < 12; p++)
temp[l, p] /= bs;
}
for (int q = 0; q < 6; q++)
{
if (q != l)
{
double bs = temp[q, l];
for (int p = l; p < 12; p++)
temp[q, p] -= bs * temp[l, p];
}
else
continue;
}
}
////得到ATA的逆阵后存在ATAR中
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
{
ATAR[i, j] = temp[i, j + 6];
}
////ATAR * AT存在ATARAT中
for (int i = 0; i < 6; i++)
for (int j = 0; j < 8; j++)
{
ATARAT[i, j] = 0;
for (int l = 0; l < 6; l++)
ATARAT[i, j] += ATAR[i, l] * AT[l, j];
}
////计算ATARAT * L，存在XX中
for (int i = 0; i < 6; i++)
for (int j = 0; j < 1; j++)
{
XX[i, j] = 0;
for (int l = 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)
{
StreamWriter mystreamwriter = new StreamWriter(ofdSaveFile.FileName);
mystreamwriter.WriteLine("计算结果为：" );
mystreamwriter.WriteLine("m:" + m.ToString());
mystreamwriter.WriteLine("H:" + H.ToString() + "m");
mystreamwriter.WriteLine("\r\n");
mystreamwriter.WriteLine("Xs:" + Xs.ToString() + "m" );
mystreamwriter.WriteLine("Ys:" + Ys.ToString() + "m");
mystreamwriter.WriteLine("Zs:" + Zs.ToString() + "m");
mystreamwriter.WriteLine("\r\n");
mystreamwriter.WriteLine("φ:" + φ.ToString() + "rad");
mystreamwriter.WriteLine("ω:" + ω.ToString() + "rad");
mystreamwriter.WriteLine("κ:" + κ.ToString() + "rad");
mystreamwriter.WriteLine("\r\n");
mystreamwriter.WriteLine("dXs:" + XX[0, 0].ToString() + "m");
mystreamwriter.WriteLine("dYs:" + XX[1, 0].ToString() + "m");
mystreamwriter.WriteLine("dZs:" + XX[2, 0].ToString() + "m");
mystreamwriter.WriteLine("dφ:" + XX[3, 0].ToString() + "rad");
mystreamwriter.WriteLine("dω:" + XX[4, 0].ToString() + "rad");
mystreamwriter.WriteLine("dκ:" + XX[5, 0].ToString() + "rad");
mystreamwriter.WriteLine("迭代次数：" + n.ToString());
mystreamwriter.Close();
}
}
else
return;
}
private void editorToolStripMenuItem_Click(object sender, EventArgs e)
{
MessageBox.Show("地信 XXX","个人信息");
}
}
}
九、运行结果
十、实习总结
通过这次摄影测量实习，是我学到了许多课本上学不到的知识和实践经验，也初步学会了关于单张影像空间后方交会的程序设计，丰富了自己的专业知识，让自己对摄影测量又有了性得了解；同时，在实习中也发现了自身的许多不足，使相关的专业知识和技巧得到了进一步提高。
摄影测量实习报告
——单张影像空间后方交会程序设计
姓名：蒋瑞
学号：
班级：
学院：矿业学院
指导教师：