单像空间后方交会实习报告.doc

摄影测量学单像空间后方交会实习报告一、 实习目的1. 掌握空间后方交会的定义和实现算法(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，。(2)算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数Xs，Ys，Zs，。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。2. 了解摄影测量平差的基本过程(1) 获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标X，Y，Z。(2) 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。(3) 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，、的初值都设为0。或者的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。(4) 计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。(5) 逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。(6) 逐点计算误差方程式的系数和常数项，组成误差方程式。(7) 计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。(8) 解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。(9) 检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对，的改正数，给予限差，通常为0.000001弧度，当3个改正数均小于0.000001弧度时，迭代结束。否则用新的近似值重复(4)（8）步骤的计算，直到满足要求为止。3. 通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。4. 实习过程：4.1学习单张像片空间后方交会的基本理论，掌握其基本思想。如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。而单像空间后方交会就是用于测定像片的外方位元素的，它的基本思想是：以单幅影像为基础，从影像所覆盖的地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,p,w,k.由于空间后方交会所采用的数学模型共线方程是非线性函数，为了便于外方位元素的解求，首先将其线性化。4.2在纸上绘出空间后方交会的计算机程序框图。为了能够在宏观上指导我们编写程序，我们需要在草稿纸上绘出程序框图。输入原始数据归算像点坐标x，y计算和确定初值Xs0，Ys0，Zs0，0，0，0组成旋转矩阵R计算(x),(y)和lx，ly逐点组成误差方程式并法化所有点完否？计算改正后的外方位元素解法方程，求未知数改正数未知数改正数限差否？整理并输出计算结果正常结束否迭代次数小于限差否？否是否输出中间结果和出错信息非正常结束二、 源代码using System;usingSystem.Collections.Generic;usingSystem.Linq;usingSystem.Text;namespace 单像空间的后方交会public class calculate private double j, k, l, Xs, Ys, Zs;/六个外方位元素 private double f = 28.1539;/主距struct point /像点和地面点坐标 public double x, y, X, Y, Z; private point p = new point4;/存取控制点的坐标 private double R = new double9;/旋转矩阵 private double a = new double8;/近似值坐标 private double L = new double8;/误差方程常数项 private double, A = new double8, 6;/误差方程系数项private int count = 0;public void Y(double q) for (int n = 0; n 4; n+) int m = n * 5;Console.WriteLine(请输入第0控制点的坐标, n + 1); qm = Convert.ToDouble(Console.ReadLine();qm + 1 = Convert.ToDouble(Console.ReadLine();qm + 2 = Convert.ToDouble(Console.ReadLine();qm + 3 = Convert.ToDouble(Console.ReadLine();qm + 4 = Convert.ToDouble(Console.ReadLine();pn.x = qm;pn.y = qm + 1;pn.X = qm + 2;pn.Y = qm + 3;pn.Z = qm + 4; double ave = 0, sum = 0;/ 求比例尺分母，用来求外方位元素for (int n = 0; n 3; n+) for (int m = n + 1; m 4; m+) sum += Math.Sqrt(Math.Pow(pn.X - pm.X, 2) + Math.Pow(pn.Y - pm.Y, 2) / Math.Sqrt(Math.Pow(pn.x - pm.x, 2) + Math.Pow(pn.y - pm.y, 2); ave = sum / 6; double j = 0.054882; k = 0.057034; l = -0.036175;/六个外方位元素的近似值doubleXs = 500215.49;doubleYs = 4185301.89;doubleZs = 1475.56; private double sin(double m) returnMath.Sin(m); private double cos(double m) returnMath.Cos(m); public void calX()/计算旋转矩阵 R0 = cos(j) * cos(k) - sin(j) * sin(l) * sin(k);R1 = -cos(j) * sin(k) - sin(j) * sin(l) * cos(k);R2 = -sin(j) * cos(l);R3 = cos(l) * sin(k);R4 = cos(l) * cos(k);R5 = -sin(l);R6 = sin(j) * cos(k) + cos(j) * sin(l) * sin(k);R7 = -sin(j) * sin(k) + cos(j) * sin(l) * cos(k);R8 = cos(j) * cos(l); public void calJSZ()/计算像点坐标近似值 for (int n = 0; n 4; n+) a2 * n = -f * (R0 * (pn.X - Xs) + R3 * (pn.Y - Ys) + R6 * (pn.Z - Zs) / (R2 * (pn.X - Xs) + R5 * (pn.Y - Ys) + R8 * (pn.Z - Zs); a2 * n + 1 = -f * (R1 * (pn.X - Xs) + R4 * (pn.Y - Ys) + R7 * (pn.Z - Zs) / (R2 * (pn.X - Xs) + R5 * (pn.Y - Ys) + R8 * (pn.Z - Zs); public void calXSJZandCSX()/计算系数矩阵和常数项 for (int n = 0; n 4; n+)/计算常数项 L2 * n = pn.x - a2 * n;L2 * n + 1 = pn.y - a2 * n + 1; for (int n = 0; n 4; n+)/计算系数矩阵 double z = R2 * (pn.X - Xs) + R5 * (pn.Y - Ys) + R8 * (pn.Z - Zs);A2 * n, 0 = (R0 * f + R2 * pn.x) / z;A2 * n, 1 = (R3 * f + R5 * pn.x) / z;A2 * n, 2 = (R6 * f + R8 * pn.x) / z;A2 * n, 3 = pn.y * sin(l) - (pn.x * (pn.x * cos(k) - pn.y * sin(k) / f + f * cos(k) * cos(l);A2 * n, 4 = -f * sin(k) - pn.x / f * (pn.x * sin(k) + pn.y * cos(k);A2 * n, 5 = pn.y;A2 * n + 1, 0 = (R1 * f + R2 * pn.y) / z;A2 * n + 1, 1 = (R4 * f + R5 * pn.y) / z;A2 * n + 1, 2 = (R7 * f + R8 * pn.y) / z;A2 * n + 1, 3 = pn.x * sin(l) - (pn.y * (pn.x * cos(k) - pn.y * sin(k) / f - f * sin(k) * cos(l);A2 * n + 1, 4 = -f * cos(k) - pn.y / f * (pn.x * sin(k) + pn.y * cos(k);A2 * n, 5 = -pn.x; public double calJSGZS()/计算改正数 double, AT = new double6, 8;/A转置 double, temp = new double6, 6;/A转置与A的乘积 double X = new double6;/改正数 double ATL = new double6;/A转置与L相乘的积int n, m, s; for (n = 0; n 8; n+)/求A的转置矩阵ATfor (m = 0; m 6; m+) ATm, n = ATn, m; for (n = 0; n 6; n+)/求A转置与A的乘积for (m = 0; m 6; m+) tempn, m = 0;for (s = 0; s 8; s+) tempn, m += ATn, s * As, m; MA(temp);/求逆 for (n = 0; n 6; n+)/求A的转置与L的乘积 for (m = 0; m 8; m+) ATLn += ATn, m * Lm; for (n = 0; n 6; n+)/求改正数X for (m = 0; m = 4) calX();calJSZ();calXSJZandCSX();Console.WriteLine(第0次迭代的结果, count);calJSGZS();count+; private void MA(double, c) int i, j, h, m;constint n = 6;double l;double, q = new doublen, 12; for (i = 0; i n; i+)/求高斯矩阵 for (j = 0; j 12; j+) qi, j = ci, j; for (i = 0; i n; i+)/ for (j = 0; j 12; j+)/构造单位阵 if (i + 6 = j)qi, j = 1;elseqi, j = 0; for (h = 0, m = 0; m n - 1; m+, h+)/消去对角线以下的数据 for (i = m + 1; i n; i+) if (qi, h = 0d)continue; l = qm, h / qi, h;for (j = 0; j 0; m-, h-)/消去对角线以上的数据 for (i = m - 1;