摄影测量学空间后方交会实验报告测绘101徐斌

摄影测量学空间后方交会实验报告测绘101徐斌摄影测量学实验报告实验一、单像空间后方交会学院: 建测学院班级: 测绘101姓名: 徐斌学号: 26一( 实验目的1.深入了解单像空间后方交会的计算过程;2.加强空间后方交会基本公式和误差方程式，法线方程式的记忆;3.通过上机调试程序加强动手能力的培养。

二(实验原理以单幅影像为基础，从该影像所覆盖地面范围内若干控制点和相应点的像坐标量测值出发，根据共线条件方程，求解该影像在航空摄影时刻的相片外方位元素。

三(实验内容1.程序图框图2.实验数据(1)已知航摄仪内方位元素f,153.24mm，Xo,Yo,0。

限差0.1秒(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.313.实验程序Form1.cs 程序using System;using System.Collections.Generic; using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;using System.IO;namespace 后方交会1{public partial class Form1 : Form {public Form1(){InitializeComponent();}public doublef,m,Xs, Ys, Zs,a1, a2, a3, b1, b2, b3, c1, c2, c3, q, w, k;public static int N,s;public double[] x = new double[4];public double[] y = new double[4];public double[] X = new double[4];public double[] Y = new double[4];public double[] Z = new double[4];public double[,] L = new double[N*2, 1];string output="外方位元素\t\n";Matrix XX;string[] a =File.ReadAllLines("d:\\控制点.txt");public void 计算N(){int cnt;using (StreamReader sr = new StreamReader(@"d:\控制点.txt")) {cnt = 0;while (sr.ReadLine() != null){cnt++;}// 这个cnt就是行数}N = cnt/5;}public void 计算初始值(){double temp=0;for (int i = 0; i <N; i++){temp+=X[i];}Xs =temp/ N;double tmp = 0;for (int i = 0; i < N; i++){tmp += Y[i];}Ys =tmp/ N;Zs = m*f;}public void 求改正数(){a1 = Math.Cos(q) * Math.Cos(k) - Math.Sin(q) * Math.Sin(w) * Math.Sin(k);a2 = -Math.Cos(q) * Math.Sin(k) - Math.Sin(q) * Math.Sin(w) * Math.Cos(k);a3 = -Math.Sin(q) * Math.Cos(w);b1 = Math.Cos(w) * Math.Sin(k);b2 = Math.Cos(w) * Math.Cos(k);b3 = -Math.Sin(w);c1 = Math.Sin(q) * Math.Cos(k) + Math.Cos(q) * Math.Sin(w) * Math.Sin(k);c2 = -Math.Sin(q) * Math.Sin(k) + Math.Cos(q) * Math.Sin(w) * Math.Cos(k);c3 = Math.Cos(q) * Math.Cos(w);int p = 0;int j = p;for (p = 0; p < a.Length; p += 5){x[j] = double.Parse(a[p])/1000;y[j] = double.Parse(a[p + 1])/1000;X[j] = double.Parse(a[p + 2]);Y[j] = double.Parse(a[p + 3]);Z[j] = double.Parse(a[p+ 4]);j++;}Matrix A = new Matrix(N*2, 6);int i = 0;int b = i;for (i = 0; i < N; i++){A.m_data[b, 0] = (a1 * f + a3 * x[i]) / (a3 * (X[i] - Xs) + b3 * (Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b, 1] = (b1 * f + b3 * x[i]) / (a3 * (X[i] - Xs) + b3 *(Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b, 2] = (c1 * f + c3 * x[i]) / (a3 * (X[i] - Xs) + b3 *(Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b, 3] = y[i] * Math.Sin(w) - (x[i] / f * (x[i] * Math.Cos(k) - y[i] * Math.Sin(k)) + f * Math.Cos(k)) * Math.Cos(w);A.m_data[b, 4] = -f * Math.Sin(k) - x[i] / f * (x[i] * Math.Sin(k) + y[i] * Math.Cos(k));A.m_data[b, 5] = y[i];A.m_data[b + 1, 0] = (a2 * f + a3 * y[i]) / (a3 * (X[i] - Xs) + b3 * (Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b + 1, 1] = (b2 * f + b3 * y[i]) / (a3 * (X[i] - Xs) + b3 * (Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b + 1, 2] = (c2 * f + c3 * y[i]) / (a3 * (X[i] - Xs) + b3 * (Y[i] - Ys) + c3 * (Z[i] - Zs));A.m_data[b + 1, 3] = -x[i] * Math.Sin(w) - (y[i] / f * (x[i] *Math.Cos(k) - y[i] * Math.Sin(k)) - f * Math.Sin(k)) * Math.Cos(w);A.m_data[b + 1, 4] = -f * Math.Cos(k) - y[i] / f * (x[i] *Math.Sin(k) + y[i] * Math.Cos(k));A.m_data[b + 1, 5] = -x[i];b += 2;}Matrix T = new Matrix(6, N*2);T = A.Transpose();Matrix AT = T * A;Matrix AA = AT.Inverse();Matrix P = AA * T;Matrix L = new Matrix(N*2, 1);int v = 0;int m = v;for (v = 0; v < N; v++){L.m_data[m, 0] = x[v] + f * (a1 * (X[v] - Xs) + b1 * (Y[v] - Ys) + c1 * (Z[v]- Zs)) / (a3 * (X[v] - Xs) + b3 * (Y[v] - Ys) + c3 * (Z[v] - Zs));L.m_data[m + 1, 0] = y[v] + f * (a2 * (X[v] - Xs) + b2 * (Y[v] - Ys) + c2 * (Z[v]- Zs)) / (a3 * (X[v] - Xs) + b3 * (Y[v] - Ys) + c3 * (Z[v] - Zs));m += 2;}XX = P * L;//计算外方位元素Xs += XX[0, 0];Ys += XX[1, 0];Zs += XX[2, 0];q += XX[3, 0];w += XX[4, 0];k += XX[5, 0];}public void 迭代(){求改正数();for (int i = 0; i < 6; i++){for (int j = 0; j < 1; j++){output += XX[i, j] + "\t\n";}}if (Math.Abs(XX[3, 0]) >= 0.000029 && Math.Abs(XX[4, 0]) >= 0.000029 && Math.Abs(XX[5,0]) >= 0.000029){for (int d = 1; d < s; d++){求改正数();for (int i = 0; i < 6; i++){for (int j = 0; j < 1; j++){output += XX[i, j]+ "\t\n";}}if (Math.Abs(XX[3, 0]) <= 0.000029 && Math.Abs(XX[4, 0]) <= 0.000029 &&Math.Abs(XX[5, 0]) <= 0.000029){输出();break;}}if (Math.Abs(XX[3, 0]) >= 0.000029 || Math.Abs(XX[4, 0]) >= 0.000029 ||Math.Abs(XX[5, 0]) >= 0.000029){MessageBox.Show("已达到指定迭代次数,迭代结束!!!");MessageBox.Show(output, "后方交会", MessageBoxButtons.OK,rmation);}}else{输出();}}public void 输出(){_Xs.Text = Xs.ToString();_Ys.Text = Ys.ToString();_Zs.Text = Zs.ToString();_q.Text = q.ToString();_w.Text = w.ToString();_k.Text = k.ToString();}private void button1_Click(object sender, EventArgs e) {try{f = double.Parse(textBox1.Text);}catch{MessageBox.Show("请先输入主距");}s = int.Parse(textBox2.Text);m = int.Parse(textBox3.Text);计算N();int i=0;int j = i;for ( i = 0; i < a.Length; i += 5){x[j] = double.Parse(a[i]);y[j] = double.Parse(a[i + 1]);X[j] = double.Parse(a[i + 2]);Y[j] = double.Parse(a[i + 3]);Z[j] = double.Parse(a[i + 4]);j++;}计算初始值();}private void button2_Click(object sender, EventArgs e){迭代();}private void button3_Click(object sender, EventArgs e){MessageBox.Show("请将点坐标文件保存到d盘根目录并保存为\"控制点 .txt\"\n要求每行一个值，五行表示一个点\n测绘101徐斌");}}}添加Matrix.cs类程序using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace 后方交会1{public class Matrix{//构造方阵public Matrix(int row){m_data = new double[row, row]; }public Matrix(int row, int col) {m_data = new double[row, col]; }//复制构造函数public Matrix(Matrix m){int row = m.Row;int col = m.Col;m_data = new double[row, col]; for (int i = 0; i < row; i++) for (int j = 0; j < col; j++) m_data[i, j] = m.m_data[i, j];}//设为单位阵public void SetUnit(){for (int i = 0; i < m_data.GetLength(0); i++) for (int j = 0; j < m_data.GetLength(1); j++) m_data[i, j] = ((i == j) ? 1 : 0);}//设置元素值public void SetValue(double d){for (int i = 0; i < m_data.GetLength(0); i++) for (int j = 0; j < m_data.GetLength(1); j++) m_data[i, j] = d;}//返中行数public int Row{get{return m_data.GetLength(0);}}//返回列数public int Col{get{return m_data.GetLength(1);}}//重载索引 ,存取数据成员public double this[int row, int col] {get{return m_data[row, col];}set{m_data[row, col] = value;}}//初等变换对调两行:ri<-->rjpublic Matrix Exchange(int i, int j) {double temp;for (int k = 0; k < Col; k++){temp = m_data[i, k];m_data[i, k] = m_data[j, k];m_data[j, k] = temp;}return this;}//初等变换第index 行乘以mulMatrix Multiple(int index, double mul){for (int j = 0; j < Col; j++){m_data[index, j] *= mul;}return this;}//初等变换第src行乘以mul加到第index行Matrix MultipleAdd(int index, int src, double mul) {for (int j = 0; j < Col; j++){m_data[index, j] += m_data[src, j] * mul;}return this;}//transpose 转置public Matrix Transpose(){Matrix ret = new Matrix(Col, Row);for (int i = 0; i < Row; i++)for (int j = 0; j < Col; j++){ret[j, i] = m_data[i, j];}return ret;}//矩阵乘public static Matrix operator *(Matrix lhs, Matrix rhs){if (lhs.Col != rhs.Row) //异常{System.Exception e = new Exception("相乘的两个矩阵的行列数不匹配"); throw e;}Matrix ret = new Matrix(lhs.Row, rhs.Col);double temp;for (int i = 0; i < lhs.Row; i++){for (int j = 0; j < rhs.Col; j++){temp = 0;for (int k = 0; k < lhs.Col; k++){temp += lhs[i, k] * rhs[k, j];}ret[i, j] = temp;}}return ret;}//功能:返回列主元素的行号//参数:row为开始查找的行号//说明:在行号[row,Col)范围内查找第row列中绝对值最大的元素，返回所在行号int Pivot(int row){int index = row;for (int i = row + 1; i < Row; i++){if (m_data[i, row] > m_data[index, row])index = i;}return index;}//逆阵:使用矩阵的初等变换，列主元素消去法public Matrix Inverse(){if (Row != Col) //异常,非方阵{System.Exception e = new Exception("求逆的矩阵不是方阵");throw e;}Matrix tmp = new Matrix(this);Matrix ret = new Matrix(Row); //单位阵ret.SetUnit();int maxIndex;double dMul;for (int i = 0; i < Row; i++){maxIndex = tmp.Pivot(i);if (tmp.m_data[maxIndex, i] == 0){System.Exception e = new Exception("求逆的矩阵的行列式的值等于0,"); throw e;}if (maxIndex != i) //下三角阵中此列的最大值不在当前行，交换{tmp.Exchange(i, maxIndex);ret.Exchange(i, maxIndex);}ret.Multiple(i, 1 / tmp[i, i]); tmp.Multiple(i, 1 / tmp[i, i]); for (int j = i + 1; j < Row; j++) {dMul = -tmp[j, i] / tmp[i, i]; tmp.MultipleAdd(j, i, dMul);ret.MultipleAdd(j, i, dMul);}}for (int i = Row - 1; i > 0; i--) {for (int j = i - 1; j >= 0; j--) {dMul = -tmp[j, i] / tmp[i, i]; tmp.MultipleAdd(j, i, dMul);ret.MultipleAdd(j, i, dMul);}}return ret;}//是方阵吗,public bool IsSquare(){return Row == Col;}//是对称阵吗,public bool IsSymmetric(){if (Row != Col)return false;for (int i = 0; i < Row; i++)for (int j = i + 1; j < Col; j++) if (m_data[i, j] != m_data[j, i]) return false;return true;}//公有数据成员public double[,] m_data;}}4(实验结果四(实验总结此次实验让我深入了解单像空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。