摄影测量学后方交会matlab实习报告

合集下载

matlab_实习报告

matlab_实习报告

matlab_实习报告在大学的学习生涯中,实习是一个非常重要的环节,它能够让我们将理论知识与实际应用相结合,提升自己的专业技能和综合素质。

本次实习,我选择了使用 MATLAB 这个强大的工具进行实践操作,通过一段时间的学习和实践,我收获颇丰。

一、实习目的本次实习的主要目的是深入了解和掌握MATLAB 软件的使用方法,能够运用其解决实际问题,并提高自己的编程能力和逻辑思维能力。

同时,通过实际项目的操作,培养自己的团队协作精神和解决问题的能力,为今后的学习和工作打下坚实的基础。

二、实习单位及岗位介绍我实习的单位是_____,在实习期间,我主要负责利用 MATLAB 进行数据分析和算法实现的相关工作。

三、实习内容及过程(一)基础学习在实习的初期,我首先对 MATLAB 的基本语法和操作进行了系统的学习。

了解了变量的定义、数据类型、矩阵运算、函数的编写等基础知识。

通过大量的练习和实例,我逐渐熟悉了 MATLAB 的编程环境,能够熟练地编写简单的程序。

例如,在学习矩阵运算时,我通过编写程序实现了矩阵的加法、乘法、求逆等操作,深刻理解了矩阵运算在数学和工程中的重要应用。

(二)项目实践在掌握了基础知识后,我开始参与实际的项目。

其中一个项目是对一组数据进行分析和处理,以提取有用的信息。

首先,我使用MATLAB 读取数据文件,并对数据进行预处理,包括去除噪声、缺失值处理等。

然后,运用统计学方法对数据进行分析,计算均值、方差、相关性等统计量。

最后,通过绘图函数将分析结果以直观的图表形式展示出来,以便更好地理解数据的特征和趋势。

在这个过程中,我遇到了很多问题。

例如,数据的格式不一致导致读取错误,算法的复杂度过高导致运行时间过长等。

通过查阅资料、请教同事和不断地调试,我最终解决了这些问题,顺利完成了项目任务。

(三)算法实现除了数据分析,我还参与了算法的实现工作。

在一个图像识别的项目中,需要使用机器学习算法对图像进行分类。

摄影测量学实习报告-实习报告.doc

摄影测量学实习报告-实习报告.doc

摄影测量学实习报告-实习报告第一篇:摄影测量学实习报告摄影测量学实习报告为期两周的摄影测量学实习今天正式结束了,虽然两周时间并不长,但是对于我来说,学到的东西远不能用时间来衡量。

在这两周里,我们完成了全数字摄影测量系统实习、数字影像分割程序编制、立体影像匹配程序编制等内容,这些东西让我们的两周很充实,很有意义。

其实刚开始时一直怀疑摄影测量学实习有什么意义,到了今天,我才发现这是有意义的。

因为通过本次实习,我们可以将课堂理论与实践相结合,使我们深入掌握摄影测量学基本概念和原理,加强摄影测量学的基本技能训练,并且培养了我们的分析问题和解决实际问题的能力。

通过使用数字摄影测量工作站,我们可以了解数字摄影测量的内定向、相对定向、绝对定向、测图过程及方法;通过开发数字影像分割程序和立体影像匹配程序,使自己掌握数字摄影测量基本方法与实现技术,为今后从事有关应用遥感技术应用和数字摄影测量打下坚实基础。

所以,就算现在觉得没什么用,但是也为将来奠定了很好的基础。

正因为如此,在这两周中我们都很认真的在学习并且完成实习任务。

其实说是两周,但时间真的更短,毕竟赶上了元旦假期,联欢晚会等一系列活动。

所以如何在短暂的时间里,更出色的完成任务,是我们必须考虑的。

记得实习动员的时候,老师花了很长时间又给我们讲了一次这次实习对我们的重要性,这很触动我们,毕竟老师的苦口婆心我们都看在眼里。

不光如此,老师又耐心的把实习要求,实习任务,实习步骤讲解了一遍,让我们大致明白了这次实习从何入手,这让本来很迷茫的我们瞬间找到了方向,也为我们接下来的工作提供了便利。

动员结束的日子,我们便进入机房,正式开始了实习。

首先我们结束了全数字摄影测量系统,这款软件是我们从来没有接触过的,所以刚开始的时候很陌生,不知道怎么用,也不知道能用来做什么。

还好,我们有老师的细致讲解,并且借助帮助向导可以解决我们很多问题。

所以在这个实习中,我们没有遇到太多困难。

让我印象深刻的是,我在做我们小组的绝对定向时,总是提示同名点数不够,就因为此,很难往下一步进行。

摄影测量学单像空间后方交会编程实习报告(精品资料).doc

摄影测量学单像空间后方交会编程实习报告(精品资料).doc

【最新整理,下载后即可编辑】摄影测量学单像空间后方交会编程实习报告班级:130x姓名:xx学号:2013302590xxx指导老师:李欣一、实习目的通过对提供的数据进行计算,输出像片的外方位元素并评定精度。

深入理解单像空间后方交会的思想,体会在有多余观测情况下,用最小二乘平差方法编程实现解求影像外方位元素的过程。

通过尝试编程实现加强编程处理问题的能力和对实习内容的理解,通过对实验结果的分析,增强综合运用所学知识解决实际问题的能力。

了解摄影测量平差的基本过程,掌握空间后方交会的定义和实现算法。

二、实习内容根据学习的单像空间后方交会的知识,用程序设计语言(C++或C语言)编写一个完整的单像空间后方交会程序,通过对提供的数据进行计算,输出像片的外方位元素并评定精度。

三、实习数据已知航摄仪的内方位元素:fk =153.24mm,x=y=0,摄影比例尺为1:15000;4个地面控制点的地面坐标及其对应像点的像片坐标:四、实习原理如果我们知道每幅影像的6个外方位元素,就能确定被摄物体与航摄影像的关系。

因此,如何获取影像的外方位元素,一直是摄影测量工作者所探讨的问题。

可采取的方法有:利用雷达、全球定位系统(GPS)、惯性导航系统(INS)以及星相摄影机来获取影像的外方位元素;也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标,根据共线条件方程,反求该影像的外方位元素,这种方法称为单幅影像的空间后方交会。

单像空间后方交会的基本思想是:以单幅影像为基础,从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发,根据共线条件方程,解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。

五、 实习流程1. 获取已知数据。

从摄影资料中查取影像比例尺1/m ,平均摄影距离(航空摄影的航高、内方位元素x 0,y 0,f ;获取控制点的空间坐标X t ,Y t ,Z t 。

2. 量测控制点的像点坐标并进行必要的影像坐标系统误差改正,得到像点坐标。

后方交会MATLAB程序实习报告

后方交会MATLAB程序实习报告
Xs=Xs0
Ys=Ys0
Zs=Zs0
fai=fai0
omig=omig0
ka=ka0
If abs(XX(4))<0.0000291&&abs(XX(5))<0.0000291&&abs(XX(6))<0.0000291
break
end
End
R=[a1,a2,a3;b1,b2,b3;c1,c2,c3];
a2=-cos(fai0)*sin(ka0)-sin(fai0)*sin(omig0)*cos(ka0);
a3=-sin(fai0)*cos(omig0);
b1=cos(omig0)*sin(ka0);
b2=cos(omig0)*cos(ka0);
b3=-sin(omig0);
c1=sin(fai0)*cos(ka0)+cos(fai0)*sin(omig0)*sin(ka0);
a16=y(h);
a21=(a2*f+a3*y(h))/Q;
a22=(b2*f+b3*y(h))/Q;
a23=(c2*f+c3*y(h))/Q;
a24=-x(h)*sin(omig0)-(y(h)/f*(x(h)*cos(ka0)-y(h)*sin(ka0))-f*sin(ka0))*cos(omig0);
④计算旋转矩阵R。利用角元素近似值计算方向余弦值,组成R阵。
⑤逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值(x),(y)。
⑥逐点计算误差方程式的系数和常数项,组成误差方程式。
⑦计算法方程的系数矩阵ATA与常数项ATL,组成法方程式。
⑧解求外方位元素。根据法方程,解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。

后方交会MATLAB程序实习报告

后方交会MATLAB程序实习报告

《摄 影 测 量 学》单像空间后方交会实习报告班 级: XXXX姓 名: X X X学 号:XXXXXXXXXXXXX指导教师: X X X一、实习目的1、掌握空间后方交会的定义和实现算法;2、了解摄影测量平差的基本过程;3、熟练MATLAB 等程序编写。

二、实习原理利用至少三个已知地面控制点的坐标),,(A A A Z Y X A 、),,(B B B Z Y X B 、),,(C C C Z Y X C ,与其影像上对应的三个像点的影像坐标),(a a y x a 、),(b b y x b 、),(c c y x c ,根据共线方程,反求该像片的外方位元素κωϕ、、、、、S S S Z Y X 。

共线条件方程式:将共线式线性化并取一次小值项得:三、解算过程①获取已知数据。

包括影像比例尺1/m,平均摄影距离(航空摄影的航高)H,内方位元素x0、y0、f,控制点的空间坐标X、Y、Z。

②量测控制点的像点坐标并进行必要的影像坐标系统误差改正,得到像点坐标。

③确定未知数的初始值。

单像空间后方交会必须给出待定参数的初始值,在竖直航空摄影且地面控制点大体对称分布的情况下,Xs0和Ys0为均值,Zs0为航高,φ、ω、κ的初值都设为0。

或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。

④计算旋转矩阵R。

利用角元素近似值计算方向余弦值,组成R阵。

⑤逐点计算像点坐标的近似值。

利用未知数的近似值按共线条件式计算控制点像点坐标的近似值(x),(y)。

⑥逐点计算误差方程式的系数和常数项,组成误差方程式。

⑦计算法方程的系数矩阵ATA与常数项ATL,组成法方程式。

⑧解求外方位元素。

根据法方程,解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。

⑨检查计算是否收敛。

将所求得的外方位元素的改正数与规定的限差比较,通常对φ,ω,κ的改正数△φ,△ω,△κ给予限差,通常为0.000001弧度,当3个改正数均小于0.000001弧度时,迭代结束。

摄影测量学后方交会matlab实习报告

摄影测量学后方交会matlab实习报告

摄影测量原理单张影像后方交会实习目录一实习目的 (3)二实习原理 (3)1. 间接平差 (3)2. 共线方程 (3)3. 单向空间后方交会 (4)三计算流程 (4)1. 求解步骤 (4)2.计算机框图 (4)四程序实现 (5)五结果分析 (6)1.外方位元素 (6)2.误差 (6)3.旋转矩阵R (7)六实习体会 (7)1. 平台的选择 (7)2.问题的解决 (7)3.心得体会 (8)七代码展示 (8)一实习目的为了增强同学们对后方交会公式的理解,培养同学们对迭代循环编程的熟悉感,本次摄影测量课间实习内容定为用C语言或其他程序编写单片空间后方交会程序,最终输出像点坐标、地面坐标、单位权中误差、外方位元素及其精度。

已知四对点的影像坐标和地面坐标如下。

内方位元素fk=153.24mm,x0=y0=0。

本次实习,我使用了matlab2014进行后方交会程序实现。

结果与参考答案一致,精度良好。

二实习原理题干中有四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标,由此可列出8个误差方程,存在2个多余观测(n=2)。

故可利用间接平差的最小二乘法则求解。

由于共线方程是非线性函数模型,为了方便计算,需要将其“线性化”。

但如果仅取泰勒级数展开式的一次项,未知数的近似值改正是不精确的。

因此必须采用迭代趋近法计算,直到外方位元素的改正值小于限差。

1.间接平差间接平差为平差计算最常用的方法。

在确定多个未知量的最或然值时,选择它们之间不存在任何条件关系的独立量作为未知量组成用未知量表达测量的函数关系、列出误差方程式,按最小二乘法原理求得未知量的最或然值的平差方法。

在一个间接平差问题中,当所选的独立参数X个数与必要观测值t个数相等时,可将每个观测值表达成这t个参数的函数,组成观测方程。

函数模型为:L = BX + d。

2.共线方程共线方程是中心投影构像的数学基础,也是各种摄影测量处理方法的重要理论基础。

式中:x,y 为像点的像平面坐标;x0,y0,f 为影像的内方位元素;XS,YS,ZS 为摄站点的物方空间坐标;XA,YA,ZA 为物方点的物方空间坐标;ai,bi,ci (i = 1,2,3)为影像的3 个外方位角元素组成的9 个方向余弦。

摄影测量学空间后方交会实验报告

摄影测量学空间后方交会实验报告

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

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

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

限差0.1秒(2)已知4对点的影像坐标和地面坐标:3.实验程序using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace ConsoleApplication3{class Program{static void Main(){//输入比例尺,主距,参与平参点的个数Console.WriteLine("请输入比例尺分母m:\r");string m1 = Console.ReadLine();double m = (double)Convert.ToSingle(m1);Console.WriteLine("请输入主距f:\r");string f1 = Console.ReadLine();double f = (double)Convert.ToSingle(f1);Console.WriteLine("请输入参与平差控制点的个数n:\r");string n1 = Console.ReadLine();int n = (int)Convert.ToSingle(n1);//像点坐标的输入代码double[] arr1 = new double[2 * n];//1.像点x坐标的输入for (int i = 0; i < n; i++){Console.WriteLine("请输入已进行系统误差改正的像点坐标的x{0}值:\r", i+1);string u = Console.ReadLine();for (int j = 0; j < n; j += 2){arr1[j] = (double)Convert.ToSingle(u);}}//2.像点y坐标的输入for (int i = 0; i < n; i++){Console.WriteLine("请输入已进行系统误差改正的像点坐标的y{0}值:\r", i+1);string v = Console.ReadLine();for (int j = 1; j < n; j += 2){arr1[j] = (double)Convert.ToSingle(v);}}//控制点的坐标输入代码double[,] arr2 = new double[n, 3];//1.控制点X坐标的输入for (int j = 0; j < n; j++){Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的X{0}值:\r", j+1);string u = Console.ReadLine();arr2[j , 0] = (double)Convert.ToSingle(u);}//2.控制点Y坐标的输入for (int k = 0; k < n; k++){Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的Y{0}值:\r", k+1);string v = Console.ReadLine();arr2[k , 1] = (double)Convert.ToSingle(v);}//3.控制点Z坐标的输入for (int p =0; p < n; p++){Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的Z{0}值:\r", p+1);string w = Console.ReadLine();arr2[p , 2] = (double)Convert.ToSingle(w);}//确定外方位元素的初始值//1.确定Xs的初始值:double Xs0 = 0;double sumx = 0;for (int j = 0; j < n; j++){double h = arr2[j, 0];sumx += h;}Xs0 = sumx / n;//2.确定Ys的初始值:double Ys0 = 0;double sumy = 0;for (int j = 0; j < n; j++){double h = arr2[j, 1];sumy += h;}Ys0 = sumy / n;//3.确定Zs的初始值:double Zs0 = 0;double sumz = 0;for (int j = 0; j <= n - 1; j++){double h = arr2[j, 2];sumz += h;}Zs0 = sumz / n;doubleΦ0 = 0;doubleΨ0 = 0;double K0 = 0;Console.WriteLine("Xs0,Ys0,Zs0,Φ0,Ψ0,K0的值分别是:{0},{1},{2},{3},{4},{5}", Xs0, Ys0, Zs0, 0, 0, 0);//用三个角元素的初始值按(3-4-5)计算各方向余弦值,组成旋转矩阵,此时的旋转矩阵为单位矩阵I:double[,] arr3 = new double[3, 3];for (int i = 0; i < 3; i++)arr3[i, i] = 1;}double a1 = arr3[0, 0]; double a2 = arr3[0, 1]; double a3 = arr3[0, 2];double b1 = arr3[1, 0]; double b2 = arr3[1, 1]; double b3 = arr3[1, 2];double c1 = arr3[2, 0]; double c2 = arr3[2, 1]; double c3 = arr3[2, 2];/*利用线元素的初始值和控制点的地面坐标,代入共线方程(3-5-2),* 逐点计算像点坐标的近似值*///1.定义存放像点近似值的数组double[] arr4 = new double[2 * n];//----------近似值矩阵//2.逐点像点坐标计算近似值//a.计算像点的x坐标近似值(x)for (int i = 0; i < 2 * n; i += 2){for (int j = 0; j < n; j++){arr4[i] = -f * (a1 * (arr2[j, 0] - Xs0) + b1 * (arr2[j, 1] - Ys0) + c1 * (arr2[j, 2] - Zs0)) / (a3 * (arr2[j, 0] - Xs0) + b3 * (arr2[j, 1] - Ys0) + c3 * (arr2[j, 2] - Zs0)); }}//b.计算像点的y坐标近似值(y)for (int i = 1; i < 2 * n; i += 2){for (int j = 0; j < n; j++){arr4[i] = -f * (a2 * (arr2[j, 0] - Xs0) + b2 * (arr2[j, 1] - Ys0) + c2 * (arr2[j, 2] - Zs0)) / (a3 * (arr2[j, 0] - Xs0) + b3 * (arr2[j, 1] - Ys0) + c3 * (arr2[j, 2] - Zs0)); }}//逐点计算误差方程式的系数和常数项,组成误差方程:double[,] arr5 = new double[2 * n, 6]; //------------系数矩阵(A)//1.计算dXs的系数for (int i = 0; i < 2 * n; i += 2){arr5[i, 0] = -1 / m; //-f/H == -1/m}//2.计算dYs的系数for (int i = 1; i < 2 * n; i += 2){arr5[i, 1] = -1 / m; //-f/H == -1/m}//3.a.计算误差方程式Vx中dZs的系数for (int i = 0; i < 2 * n; i += 2)arr5[i, 2] = -arr1[i] / m * f;}//3.b.计算误差方程式Vy中dZs的系数for (int i = 1; i < 2 * n; i += 2){arr5[i, 2] = -arr1[i] / m * f;}//4.a.计算误差方程式Vx中dΦ的系数for (int i = 0; i < 2 * n; i += 2){arr5[i, 3] = -f * (1 + arr1[i] * arr1[i] / f * f);}//4.a.计算误差方程式Vy中dΦ的系数for (int i = 1; i < 2 * n; i += 2){arr5[i, 3] = -arr1[i - 1] * arr1[i] / f;}//5.a.计算误差方程式Vx中dΨ的系数for (int i = 0; i < 2 * n; i += 2){arr5[i, 4] = -arr1[i] * arr1[i + 1] / f;}//5.b.计算误差方程式Vy中dΨ的系数for (int i = 1; i < 2 * n; i += 2){arr5[i, 4] = -f * (1 + arr1[i] * arr1[i] / f * f);}//6.a.计算误差方程式Vx中dk的系数for (int i = 0; i < 2 * n; i += 2){arr5[i, 5] = arr1[i + 1];}//6.b.计算误差方程式Vy中dk的系数for (int i = 1; i < 2 * n; i += 2){arr5[i, 5] = -arr1[i - 1];}//定义外方位元素组成的数组double[] arr6 = new double[6];//--------------------外方位元素改正数矩阵(X)//定义常数项元素组成的数组double[] arr7 = new double[2 * n];//-----------------常数矩阵(L)//计算lx的值for (int i = 0; i < 2 * n; i += 2)arr7[i] = arr1[i] - arr4[i]; //将近似值矩阵的元素代入}//计算ly的值for (int i = 1; i <= 2 * (n - 1); i += 2){arr7[i] = arr1[i] - arr4[i]; //将近似值矩阵的元素代入}/* 对于所有像点的坐标观测值,一般认为是等精度量测,所以权阵P为单位阵.所以X=(ATA)-1ATL *///1.计算ATdouble[,] arr5T = new double[6, 2 * n];for (int i = 0; i < 6; i++){for (int j = 0; j < 2 * n; j++){arr5T[i, j] = arr5[j, i];}}//A的转置与A的乘积,存放在arr5AA中double[,] arr5AA = new double[6, 6];for (int i = 0; i < 6; i++){for (int j = 0; j < 6; j++){arr5AA[i, j] = 0;for (int l = 0; l < 2 * n; l++){arr5AA[i, j] += arr5T[i, l] * arr5[l, j];}}}nijuzhen(arr5AA);//arr5AA经过求逆后变成原矩阵的逆矩阵//arr5AA * arr5T存在arr5AARATdouble[,] arr5AARAT = new double[6, 2 * n];for (int i = 0; i < 6; i++){for (int j = 0; j < 2 * n; j++){arr5AARAT[i, j] = 0;for (int p = 0; p < 6; p++){arr5AARAT[i, j] += arr5AA[i, p] * arr5T[p, j];}}}//计算arr5AARAT x L,存在arrX中double[] arrX = new double[6];for (int i = 0; i < 6; i++){for (int j = 0; j < 1; j++){arrX[i] = 0;for (int vv = 0; vv < 6; vv++){arrX[i] += arr5AARAT[i, vv] * arr7[vv];}}}//计算外方位元素值double Xs, Ys, Zs, Φ, Ψ, K;Xs = Xs0 + arrX[0];Ys = Ys0 + arrX[1];Zs = Zs0 + arrX[2];Φ = Φ0 + arrX[3];Ψ = Ψ0 + arrX[4];K = K0 + arrX[5];for (int i = 0; i <= 2; i++){Xs += arrX[0];Ys += arrX[1];Zs += arrX[2];Φ += arrX[3];Ψ += arrX[4];K += arrX[5];}Console.WriteLine("Xs,Ys,Zs,Φ,Ψ,K的值分别是:{0},{1},{2},{3},{4},{5}", Xs0, Ys0, Zs0, Φ, Ψ, K);Console.Read();}//求arr5AA的逆矩public static double[,] nijuzhen(double[,] a) {double[,] B = new double[6, 6];int i, j, k;int row = 0;int col = 0;double max, temp;int[] p = new int[6];for (i = 0; i < 6; i++){p[i] = i;B[i, i] = 1;}for (k = 0; k < 6; k++){//找主元max = 0; row = col = i;for (i = k; i < 6; i++){for (j = k; j < 6; j++){temp = Math.Abs(a[i, j]);if (max < temp){max = temp;row = i;col = j;}}}//交换行列,将主元调整到k行k列上if (row != k){for (j = 0; j < 6; j++){temp = a[row, j];a[row, j] = a[k, j];a[k, j] = temp;temp = B[row, j];B[row, j] = B[k, j];B[k, j] = temp;i = p[row]; p[row] = p[k]; p[k] = i; }if (col != k){for (i = 0; i < 6; i++){temp = a[i, col];a[i, col] = a[i, k];a[i, k] = temp;}}//处理for (j = k + 1; j < 6; j++){a[k, j] /= a[k, k];}for (j = 0; j < 6; j++){B[k, j] /= a[k, k];a[k, k] = 1;}for (j = k + 1; j < 6; j++){for (i = 0; j < k; i++){a[i, j] -= a[i, k] * a[k, j];}for (i = k + 1; i < 6; i++){a[i, j] -= a[i, k] * a[k, j];}}for (j = 0; j < 6; j++){for (i = 0; i < k; i++){B[i, j] -= a[i, k] * B[k, j];}for (i = k + 1; i < 6; i++){B[i, j] -= a[i, k] * B[k, j];}for (i = 0; i < 6; i++) {a[i, k] = 0;a[k, k] = 1;}}//恢复行列次序for (j = 0; j < 6; j++){for (i = 0; i < 6; i++) {a[p[i], j] = B[i, j]; }}for (i = 0; i < 6; i++){for (j = 0; j < 6; j++) {a[i, j] = a[i, j];}}return a;}4.实验结果四.实验总结此次实验让我深入了解单像空间后方交会的计算过程,加强了对空间后方交会基本公式和误差方程式,法线方程式的记忆。

摄影测量学单像空间后方交会编程实习报告

摄影测量学单像空间后方交会编程实习报告

摄影测量学单像空间后方交会编程实习报告实习背景在本次实习中,我们学习了摄影测量学单像空间后方交会的编程实现。

这是一种通过计算影像中各点的坐标来确定被摄物的三维坐标的方法,应用广泛于测绘、地理信息、建筑等领域。

本次实习采用 MATLAB 软件进行编程,目的是将理论知识应用到实际操作中,让我们更深入理解摄影测量学单像空间后方交会的原理和应用。

实习内容理论部分首先,我们在工作室进行了理论部分的学习。

老师讲解了单像空间后方交会的原理,以及如何通过影像坐标、相机外方位元素、像点坐标和像平面坐标等参数来计算被摄物的三维坐标。

在理论部分的学习过程中,我们通过公式的推导和实例分析,更加深入地理解了单像空间后方交会的原理。

实践部分实践部分是本次实习的重头戏。

我们利用 MATLAB 软件进行了单像空间后方交会的编程实现,具体步骤如下:1.输入相机外方位元素通过读取文本文件,将相机外方位元素(相机在拍摄时的姿态、位置等参数)输入到 MATLAB 中。

2.输入影像坐标通过读取文本文件,将影像中的像点坐标输入到 MATLAB 中。

3.计算像平面坐标利用相机内定标参数,将像点坐标转化为像平面坐标。

4.计算被摄物三维坐标根据单像空间后方交会的原理,利用相机外方位元素、像平面坐标和像点坐标等参数,计算被摄物的三维坐标。

5.输出结果将计算结果输出到文本文件中,以便后续的数据处理和分析。

在实际操作中,我们首先编写了 MATLAB 脚本文件,根据上述步骤逐步实现了单像空间后方交会的计算过程。

然后,我们利用自己拍摄的实际照片进行实验,将相机外方位元素和像点坐标输入到程序中,最终得到了被摄物的三维坐标结果。

实习收获通过本次实习,我从理论到实践,更深入地理解了摄影测量学单像空间后方交会的原理和应用,同时也掌握了 MATLAB 的编程技能。

在实践中,我遇到了许多问题,包括数据的输入输出、代码的调试和结果的分析等等。

通过和同学的讨论和老师的指导,我不仅解决了这些问题,还对摄影测量学的应用有了更深入的认识。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

摄影测量原理单张影像后方交会实习目录一实习目的 (3)二实习原理 (3)1. 间接平差 (3)2. 共线方程 (3)3. 单向空间后方交会 (4)三计算流程 (4)1. 求解步骤 (4)2.计算机框图 (4)四程序实现 (5)五结果分析 (6)1.外方位元素 (6)2.误差 (6)3.旋转矩阵R (7)六实习体会 (7)1. 平台的选择 (7)2.问题的解决 (7)3.心得体会 (8)七代码展示 (8)一实习目的为了增强同学们对后方交会公式的理解,培养同学们对迭代循环编程的熟悉感,本次摄影测量课间实习内容定为用C语言或其他程序编写单片空间后方交会程序,最终输出像点坐标、地面坐标、单位权中误差、外方位元素及其精度。

已知四对点的影像坐标和地面坐标如下。

内方位元素fk=153.24mm,x0=y0=0。

本次实习,我使用了matlab2014进行后方交会程序实现。

结果与参考答案一致,精度良好。

二实习原理题干中有四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标,由此可列出8个误差方程,存在2个多余观测(n=2)。

故可利用间接平差的最小二乘法则求解。

由于共线方程是非线性函数模型,为了方便计算,需要将其“线性化”。

但如果仅取泰勒级数展开式的一次项,未知数的近似值改正是不精确的。

因此必须采用迭代趋近法计算,直到外方位元素的改正值小于限差。

1.间接平差间接平差为平差计算最常用的方法。

在确定多个未知量的最或然值时,选择它们之间不存在任何条件关系的独立量作为未知量组成用未知量表达测量的函数关系、列出误差方程式,按最小二乘法原理求得未知量的最或然值的平差方法。

在一个间接平差问题中,当所选的独立参数X个数与必要观测值t个数相等时,可将每个观测值表达成这t个参数的函数,组成观测方程。

函数模型为:L = BX + d。

2.共线方程共线方程是中心投影构像的数学基础,也是各种摄影测量处理方法的重要理论基础。

式中:x,y 为像点的像平面坐标;x0,y0,f 为影像的内方位元素;XS,YS,ZS 为摄站点的物方空间坐标;XA,YA,ZA 为物方点的物方空间坐标;ai,bi,ci (i = 1,2,3)为影像的3 个外方位角元素组成的9 个方向余弦。

3.单向空间后方交会利用至少三个已知地面控制点的坐标A、B、Z与其影像上对应的三个像点的影像坐标a、b、c,根据共线方程,反求该像点的外方位元素。

这种解算方法以单张相片为基础,亦称单像空间后方交会。

三计算流程1.求解步骤(1)获取已知数据。

(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正,得到像点坐标。

(3)确定未知数的初始值。

在竖直航空摄影且地面控制点大体对称分布的情况下,可按如下方法确定初始值:Z S0 = H = m*fX S0 = 1/n*(X1+X2+…+X n)Y S0 = 1/n*(Y1+Y2+…+Y n)q0 = w0 = 0k0可在航迹图上找出,或根据控制点坐标通过坐标正反变换求出。

(4)利用角元素近似值计算方向余弦值,组成旋转矩阵R。

(5)利用未知数的近似值按共线方程条件式计算控制点像点坐标的近似值x0,y0。

(6)逐点计算误差方程式的系数和常数项,组成误差方程式。

(7)解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。

(8)规定限差限制迭代次数,检查计算是否收敛。

2.计算机框图四程序实现1.定义各种参数名并赋值。

由于题目是近似垂直观测,角元素q,w,k的初值我统一赋了0。

从最终结果也可以看出三个值都非常接近0。

2.读取坐标新建名为“坐标.txt”的文本文件,按顺序输入每个点的影像坐标和地面坐标。

使用matlab读取文件语句textread将x,y,X,Y,Z分别存储到对应矩阵中。

3.从计算R阵开始,直到外方位元素改正,需要进行数次迭代循环。

结束循环的条件为线元素最终改正值绝对值小于limit1(0.001),角元素最终改正值绝对值小于limit1(0.001),两者必须同时成立。

4.进入循环后,先按照书上公式,用程序中外方位元素赋的初值计算旋转矩阵R。

5.按照公式计算共线方程。

共线方程左侧应为x-x0/y-y0,但本题中x0、y0都为0,故略去一步,直接记录左侧数值为x/y。

存储在矩阵的顺序与c矩阵相同,方便二者作差求出L阵。

6.利用一个i(1~4)的循环,计算出系数阵A。

在此之前需要计算Zbar。

7.计算常数项L,进而推导出外方位元素改正值dx。

8.改正外方位元素大小,判断是否满足循环条件。

若不满足则输出结果。

9.循环结束后,分别输出外方位元素值、外方位元素误差和R阵,保存在txt文件中。

五结果分析1.外方位元素结果如下;迭代次数为5次。

结果与参考答案一致,迭代次数也符合要求,说明外方位元素初值选取的比较好。

本张像片近似垂直摄影。

2.外方位元素误差都比较小,同时角元素精度远高于线元素精度。

3.旋转矩阵R各元素大小四舍五入后与参考答案一致。

六实习体会七代码展示%单幅影像后方交会%线元素10^-3,角元素10^-6close all; clear all; clc;disp('后方交会求解')%屏幕输出[x,y,X,Y,Z]=textread('坐标.txt','%f %f %f %f %f');%读取文件x=x/1000;y=y/1000;%修正单位fk=153.24/1000;%主距修正单位m=50000;%比例尺系数limit1=0.001;limit2=0.000001;%改正值精度Xs=sum(X)/4;Ys=sum(Y)/4;Zs=fk*m+sum(Z)/4;%待定参数初始值q=0;w=0;k=0;%外方位角元素R=zeros(3,3);%旋转矩阵dx=ones(6,1);%改正值H=m*fk;%航高round=0;%迭代次数c=[x(1);y(1);x(2);y(2);x(3);y(3);x(4);y(4)];%方便L计算PI=3.1415926;%循环条件:线元素改正值<=limit1且角元素改正值<=limit2while(abs(dx(1,1))>limit1||abs(dx(2,1))>limit1||abs(dx(3,1))>limit1||a bs(dx(4,1))>limit2||abs(dx(5,1))>limit2||abs(dx(6,1))>limit2)%迭代计算旋转矩阵RR(1,1)=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);R(1,2)=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);R(1,3)=-sin(q)*cos(w);R(2,1)=cos(w)*sin(k);R(2,2)=cos(w)*cos(k);R(2,3)=-sin(w);R(3,1)=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);R(3,2)=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);R(3,3)=cos(q)*cos(w);xy_0=zeros(8,1);%共线方程。

奇数为x-x0,偶数为y-y0xy_0(1,1)=-fk*(R(1,1)*(X(1)-Xs)+R(2,1)*(Y(1)-Ys)+R(3,1)*(Z(1)-Zs))/(R( 1,3)*(X(1)-Xs)+R(2,3)*(Y(1)-Ys)+R(3,3)*(Z(1)-Zs));xy_0(2,1)=-fk*(R(1,2)*(X(1)-Xs)+R(2,2)*(Y(1)-Ys)+R(3,2)*(Z(1)-Zs))/(R( 1,3)*(X(1)-Xs)+R(2,3)*(Y(1)-Ys)+R(3,3)*(Z(1)-Zs));xy_0(3,1)=-fk*(R(1,1)*(X(2)-Xs)+R(2,1)*(Y(2)-Ys)+R(3,1)*(Z(2)-Zs))/(R( 1,3)*(X(2)-Xs)+R(2,3)*(Y(2)-Ys)+R(3,3)*(Z(2)-Zs));xy_0(4,1)=-fk*(R(1,2)*(X(2)-Xs)+R(2,2)*(Y(2)-Ys)+R(3,2)*(Z(2)-Zs))/(R( 1,3)*(X(2)-Xs)+R(2,3)*(Y(2)-Ys)+R(3,3)*(Z(2)-Zs));xy_0(5,1)=-fk*(R(1,1)*(X(3)-Xs)+R(2,1)*(Y(3)-Ys)+R(3,1)*(Z(3)-Zs))/(R( 1,3)*(X(3)-Xs)+R(2,3)*(Y(3)-Ys)+R(3,3)*(Z(3)-Zs));xy_0(6,1)=-fk*(R(1,2)*(X(3)-Xs)+R(2,2)*(Y(3)-Ys)+R(3,2)*(Z(3)-Zs))/(R( 1,3)*(X(3)-Xs)+R(2,3)*(Y(3)-Ys)+R(3,3)*(Z(3)-Zs));xy_0(7,1)=-fk*(R(1,1)*(X(4)-Xs)+R(2,1)*(Y(4)-Ys)+R(3,1)*(Z(4)-Zs))/(R( 1,3)*(X(4)-Xs)+R(2,3)*(Y(4)-Ys)+R(3,3)*(Z(4)-Zs));xy_0(8,1)=-fk*(R(1,2)*(X(4)-Xs)+R(2,2)*(Y(4)-Ys)+R(3,2)*(Z(4)-Zs))/(R( 1,3)*(X(4)-Xs)+R(2,3)*(Y(4)-Ys)+R(3,3)*(Z(4)-Zs));A=zeros(8,6);%系数阵for i=1:4Xbar=R(1,1)*(X(i)-Xs)+R(2,1)*(Y(i)-Ys)+R(3,1)*(Z(i)-Zs);Ybar=R(1,2)*(X(i)-Xs)+R(2,2)*(Y(i)-Ys)+R(3,2)*(Z(i)-Zs);Zbar=R(1,3)*(X(i)-Xs)+R(2,3)*(Y(i)-Ys)+R(3,3)*(Z(i)-Zs);A(2*i-1,1)=1/Zbar*(R(1,1)*fk+R(1,3)*xy_0(2*i-1,1));A(2*i-1,2)=1/Zbar*(R(2,1)*fk+R(2,3)*xy_0(2*i-1,1));A(2*i-1,3)=1/Zbar*(R(3,1)*fk+R(3,3)*xy_0(2*i-1,1));A(2*i,1)=1/Zbar*(R(1,2)*fk+R(1,3)*xy_0(2*i,1));A(2*i,2)=1/Zbar*(R(2,2)*fk+R(2,3)*xy_0(2*i,1));A(2*i,3)=1/Zbar*(R(3,2)*fk+R(3,3)*xy_0(2*i,1));A(2*i-1,4)=xy_0(2*i,1).*sin(w)-(xy_0(2*i-1,1)./fk.*(xy_0(2*i-1,1).*cos (k)-xy_0(2*i,1).*sin(k))+fk*cos(k))*cos(w);A(2*i-1,5)=-fk*sin(k)-xy_0(2*i-1,1)./fk.*(xy_0(2*i-1,1).*sin(k)+xy_0(2 *i,1).*cos(k));A(2*i-1,6)=xy_0(2*i,1);A(2*i,4)=-xy_0(2*i-1,1).*sin(w)-(xy_0(2*i,1)./fk.*(xy_0(2*i-1,1).*cos( k)-xy_0(2*i,1).*sin(k))-fk*sin(k))*cos(w);A(2*i,5)=-fk*cos(k)-xy_0(2*i,1)./fk.*(xy_0(2*i-1,1).*sin(k)+xy_0(2*i,1 ).*cos(k));A(2*i,6)=-xy_0(i*2-1,1);endround=round+1;L=c-xy_0;dx=((A')*A)\(A')*L;%改正外方位元素Xs=Xs+dx(1,1);Ys=Ys+dx(2,1);Zs=Zs+dx(3,1);q=q+dx(4,1);w=w+dx(5,1);k=k+dx(6,1);end%写入文件fid1=fopen('后方交会.txt','w');%输出外方位元素fprintf(fid1,'Xs = %f\r\n',Xs);fprintf(fid1,'Ys = %f\r\n',Ys);fprintf(fid1,'Zs = %f\r\n\r\n',Zs);fprintf(fid1,'q = %f\r\n',q);fprintf(fid1,'w = %f\r\n',w);fprintf(fid1,'k = %f\r\n\r\n',k);fprintf(fid1,'round = %d\r\n',round);fclose(fid1);fid2=fopen('R阵.txt','w');%输出矩阵Rfprintf(fid2,'\r\n');for m=1:3for n=1:3fprintf(fid2,'%f ',R(m,n));endfprintf(fid2,'\r\n');endfclose(fid2);fid3=fopen('精度评定.txt','w');%输出中误差fprintf(fid3,'dXs = ±%f(mm)\r\n',1000*abs(dx(1,1)));fprintf(fid3,'dYs = ±%f(mm)\r\n',1000*abs(dx(2,1)));fprintf(fid3,'dZs = ±%f(mm)\r\n',1000*abs(dx(3,1)));fprintf(fid3,'\r\ndq = ±%f(″)\r\n',abs(dx(4,1))*180*3600/PI);fprintf(fid3,'dw = ±%f(″)\r\n',abs(dx(5,1))*180*3600/PI); fprintf(fid3,'dk = ±%f(″)\r\n',abs(dx(6,1))*180*3600/PI); fclose(fid3);。

相关文档
最新文档