前方后方空间交会实验报告
单像空间后方交会实习报告
单像空间后方交会实习报告摘要本报告旨在总结并评估笔者在单像空间后方交会实习中的经验和收获。
首先，报告介绍了单像空间后方交会实习的目的和背景。
接着，报告详细描述了实习期间所进行的实验和操作步骤。
在实习过程中，笔者遇到了一些挑战，但通过团队合作和专业指导取得了成功。
最后，报告总结了实习对于个人职业发展的重要性，并提出了改进实习体验的建议。
1. 引言单像空间后方交会是测量和分析地球或其他星球上的点的空间坐标的方法之一。
该方法通过将来自不同位置的图像投影到一个共同的平面上，并在该平面上对图像进行测量和分析，以确定点的坐标。
本实习旨在将现实生活中的实地测量和图像处理技术相结合，通过实际操作了解和掌握单像空间后方交会的原理和应用。
2. 实习过程本次实习分为三个步骤：图像获取、图像处理和空间坐标计算。
2.1 图像获取首先，为了进行后续的图像处理和分析，我们需要获取一组具有不同视角的图像。
为了实现这个目标，我们选择了一片公共景区进行实地测量。
在测量过程中，我们使用了专业的测量设备和相机，并按照一定的间隔和角度拍摄了一组图像。
这些图像将被用于后续的图像处理和分析。
2.2 图像处理在图像处理阶段，我们使用了专业的图像处理软件对获取到的图像进行处理。
首先，我们使用了相机标定算法对相机内外参数进行校准，以保证后续测量的精度和准确度。
然后，我们对每张图像进行了特征点提取和匹配，以建立图像之间的对应关系。
最后，根据所获得的对应关系，我们重建了图像场景的三维模型，并将其用于后续的空间坐标计算。
2.3 空间坐标计算在空间坐标计算阶段，我们使用了单像空间后方交会的原理，计算了每个图像特征点的空间坐标。
首先，我们将图像场景的三维模型与图像上的特征点进行对应，以确定特征点在三维空间中的位置。
然后，我们利用三角测量原理计算出特征点的三维坐标。
最后，通过对所有图像特征点的计算，我们可以得到目标点的空间坐标。
3. 实习挑战与解决在实习过程中，我们遇到了一些挑战，如图像质量、算法调优和测量误差等。
后方交会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弧度时，迭代结束。
摄影测量学空间后方交会实验报告
摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘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 的编程技能。
在实践中，我遇到了许多问题，包括数据的输入输出、代码的调试和结果的分析等等。
通过和同学的讨论和老师的指导，我不仅解决了这些问题，还对摄影测量学的应用有了更深入的认识。
单像空间后方交会实习报告
框图。
3
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，系电，通力根1保过据护管生高线产中敷工资设艺料技高试术中卷0资不配料仅置试可技卷以术要解是求决指，吊机对顶组电层在气配进设置行备不继进规电行范保空高护载中高与资中带料资负试料荷卷试下问卷高题总中2体2资，配料而置试且时卷可，调保需控障要试各在验类最；管大对路限设习度备题内进到来行位确调。保整在机使管组其路高在敷中正设资常过料工程试况1卷中下安，与全要过，加度并强工且看作尽护下可1都关能可于地以管缩正路小常高故工中障作资高；料中对试资于卷料继连试电接卷保管破护口坏进处范行理围整高，核中或对资者定料对值试某，卷些审弯异核扁常与度高校固中对定资图盒料纸位试，置卷编.工保写况护复进层杂行防设自腐备动跨与处接装理地置，线高尤弯中其曲资要半料避径试免标卷错高调误等试高，方中要案资求，料技编试术写5、卷交重电保底要气护。设设装管备备置线4高、调动敷中电试作设资气高，技料课中并3术试、件资且中卷管中料拒包试路调试绝含验敷试卷动线方设技作槽案技术，、以术来管及避架系免等统不多启必项动要方高式案中，；资为对料解整试决套卷高启突中动然语过停文程机电中。气高因课中此件资，中料电管试力壁卷高薄电中、气资接设料口备试不进卷严行保等调护问试装题工置，作调合并试理且技利进术用行，管过要线关求敷运电设行力技高保术中护。资装线料置缆试做敷卷到设技准原术确则指灵：导活在。。分对对线于于盒调差处试动，过保当程护不中装同高置电中高压资中回料资路试料交卷试叉技卷时术调，问试应题技采，术用作是金为指属调发隔试电板人机进员一行，变隔需压开要器处在组理事在；前发同掌生一握内线图部槽 纸故内资障，料时强、，电设需回备要路制进须造行同厂外时家部切出电断具源习高高题中中电资资源料料，试试线卷卷缆试切敷验除设报从完告而毕与采，相用要关高进技中行术资检资料查料试和，卷检并主测且要处了保理解护。现装场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
空间后方—前方交会的原理
空间后方—前方交会的原理
以空间后方—前方交会的原理为题，我来为大家描述一下。
空间后方—前方交会是一种用于确定目标位置的方法，常用于航空、导航、测绘等领域。
它利用人眼的立体视觉和视差效应，通过观察目标在不同视角下的位置变化，来推断目标的实际位置。
这种方法可以较精确地确定目标的距离和方位，尤其适用于远距离观测。
在进行空间后方—前方交会时，我们首先需要选择两个观测点，它们之间的距离应足够远，以便产生明显的视差效应。
然后，我们分别在这两个观测点上观察目标，并记录下目标在两个观测点的位置。
接下来，我们需要测量观测点之间的距离，并确定观测点与目标之间的夹角。
这些数据将用于计算目标的实际位置。
通过对两个观测点的位置和距离进行几何分析，我们可以得到目标相对于观测点的位移向量。
然后，我们再将这个位移向量与观测点之间的夹角结合起来，就可以计算出目标相对于观测点的实际位置。
空间后方—前方交会的原理基于视差效应，即当我们观察远处的目标时，由于两只眼睛的视角不同，目标在两只眼睛中的位置也会有所不同。
通过比较这两个位置的差异，我们就可以推断出目标的实际位置。
总的来说，空间后方—前方交会是一种利用视差效应来确定目标位
置的方法。
它可以在远距离观测中提供较为准确的测量结果，具有广泛的应用前景。
空间后交-前交程序设计实验报告
空间后交-前交程序设计(实验报告)姓名：班级：学号：时间：空间后交-前交程序设计一、实验目的用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序⑴提交实习报告：程序框图、程序源代码、计算结果、体会⑵计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度二、实验数据f=150.000mm，x0=0，y0=0三、实验思路1.利用空间后方交会求左右像片的外方位元素(1).获取m(于像片中选取两点，于地面摄影测量坐标系中选取同点，分别计算距离，距离比值即为m)，x,y,f,X,Y,Z(2).确定未知数初始值Xs,Ys,Zs,q,w,k(3).计算旋转矩阵R(4).逐点计算像点坐标的近似值（x），（y）(5).组成误差方程式(6).组成法方程式(7).解求外方位元素(8).检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤（3）-（7）2.利用求出的外方位元素进行空间前交，求出待定点地面坐标(1).用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1,R2(2).根据左、右像片的外方位元素，计算摄影基线分量Bx，By，Bz(3).计算像点的像空间辅助坐标（X1,Y1,Z1）和（X2,Y2,Z2）(4).计算点投影系数N1和N2(5).计算未知点的地面摄影测量坐标四、实验过程⑴程序框图函数AandL%求间接平差时需要的系数%%%已知%a=像点坐标x,b=像点坐标y，f内方位元素主距%φ=q,ψ=w,κ=k%像空间坐标系X,Y,Z%地面摄影测量坐标系Xs,Ys,Zsfunction [A1,L1,A2,L2]=AandL(a,b,f,q,w,k,X,Y,Z,Xs,Ys,Zs) %%%%%%%%%%%选择矩阵元素a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a3=-sin(q)*cos(w);b1=cos(w)*sin(k);b2=cos(w)*cos(k);b3=-sin(w);c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c3=cos(q)*cos(w);%%%%%%%共线方程的分子分母X_=a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs);Y_=a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs);Z_=a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs);%%%%%%%近似值x=-f*X_/Z_;y=-f*Y_/Z_;%%%%%%%A组成L组成a11=1/Z_*(a1*f+a3*x);a12=1/Z_*(b1*f+b3*x);a13=1/Z_*(c1*f+c3*x);a21=1/Z_*(a2*f+a3*y);a22=1/Z_*(b2*f+b3*y);a23=1/Z_*(c2*f+c3*y);a14=y*sin(w)-(x/f*(x*cos(k)-y*sin(k))+f*cos(k))*cos(w);a15=-f*sin(k)-x/f*(x*sin(k)+y*cos(k)); a16=y;a24=-x*sin(w)-(y/f*(x*cos(k)-y*sin(k))-f*sin(k))*cos(w);a25=-f*cos(k)-y/f*(x*sin(k)+y*cos(k));a26=-x;lx=a-x;ly=b-y;%%%%%%%%%组成一个矩阵，并返回A1=[a11,a12,a13,a14,a15,a16];A2=[a21,a22,a23,a24,a25,a26];L1=lx;L2=ly;函数deg2dms%%%%%%%%角度转度分秒function y=deg2dms(x)a=floor(x);b=floor((x-a)*60);c=(x-a-b/60)*3600;y=a+(b/100)+(c/10000);函数dms2deg%%%%%度分秒转度function y=dms2deg(x)a=floor(x);b=floor((x-a)*100);c=(x-a-b/100)*10000;y=a+b/60+c/3600;函数ok%%%%%%%%%%%%%%目的是为了保证各取的值的有效值%%xy为n*1,a为1*nfunction result=ok(xy,a)format short gi=size(xy,1);for n=1:io=xy(n)-floor(xy(n,1));o=round(o*(10^a(n)))/(10^a(n));xy(n,1)=floor(xy(n,1))+o;endformat long gresult=xy;函数rad2dmsxy%%%%求度分秒表现形式的三个外方位元素，三个角度function xydms=rad2dmsxy(xy)[a,b,c,d,e,f]=testvar(xy);d=deg2dms(rad2deg(d));e=deg2dms(rad2deg(e));f=deg2dms(rad2deg(f));xydms=[a,b,c,d,e,f]';函数spacehoujiao%%%%%%%空间后交%%% f%%输入p（2*n,1）%%像点坐标x,y,X,Y,Z,均为（n,1）function [xy,m,R]=spacehoujiao(p,x,y,f,X,Y,Z)format long;%%%%%权的矢量化,这是等精度时的，如果非，将函数参数改为P P=diag(p);%%求nj=size(X,2);%%初始化Xs=0;Ys=0;Zs=0;for n=1:jXs=Xs+X(n);Ys=Ys+Y(n);Zs=Zs+Z(n);endSx=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);%%%%两像点之间距离Sd=sqrt((X(2)-X(1))^2+(Y(2)-Y(1))^2);%%%%两地面控制点之间距离m=Sd/Sx; %%%%图像比例系数Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j;m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6);l=rand(2*j,1);%%%%for n=1:j[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n ),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;%%%%%%%%%循环体while 1%%%%%%%%%%%%%%%%[dXs,dYs,dZs,dq,dw,dk]=testvar(det);detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);%%%%%%%%%if((detXs<0.01)&&(detYs<0.01)&&(detZs<0.01)&&(detq<pi/648000)&&(detw<pi/648000)&& (detk<pi/648000))break;elseV=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt((V'*P*V)/(2*j-6));%%m0需要每次的改正数算出来相加%%%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%%%for n=1:j[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n ),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;i=i+1;%%%%end%%%end[dXs,dYs,dZs,dq,dw,dk]=testvar(det);detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);V=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt((V'*P*V)/(2*n-6));%%%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%%%%%%%%%%%%%可以输出迭代次数的i%%%%%%%%%%%%Xs,Ys,Zs,q,w,k,i,dXs,dYs,dZs,dq,dw,dk,detXs,detYs,detZs %%%%%%%%%%%精度mo=m0*sqrt(Q);m=[mo(1,1),mo(2,2),mo(3,3),mo(4,4),mo(5,5),mo(6,6)]';[mXs,mYs,mZs,mq,mw,mk]=testvar(m);%%%%%%%%%输出xy=[Xs,Ys,Zs,q,w,k]';%%输出（6,1）的外方位元素m=[m0,mXs,mYs,mZs,mq,mw,mk]';%%单位误差，各元素中误差R=xyR(xy);%%旋转矩阵函数spaceqianjiao%%空间前交%输入f%输入x1,y1,x2,y2,R1,R2,xy1,xy2 （n,1）%输出X,Y,Z （n,1）function [X,Y,Z]=spaceqianjiao(x1,y1,x2,y2,f,R1,R2,xy1,xy2) i=size(x1,2);[Xs1,Ys1,Zs1,q1,w1,k1]=testvar(xy1);[Xs2,Ys2,Zs2,q2,w2,k2]=testvar(xy2);for n=1:i[X1(n),Y1(n),Z1(n)]=testvar(R1*[x1(n),y1(n),-f]');[X2(n),Y2(n),Z2(n)]=testvar(R2*[x2(n),y2(n),-f]');Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;N1=(Bx*Z2(n)-Bz*X2(n))/(X1(n)*Z2(n)-X2(n)*Z1(n));N2=(Bx*Z1(n)-Bz*X1(n))/(X1(n)*Z2(n)-X2(n)*Z1(n));X(n)=Xs1+N1*X1(n);Z(n)=Zs1+N1*Z1(n);Y(n)=0.5*((Ys1+N1*Y1(n))+(Ys2+N2*Y2(n)));end函数testvar%分割矩阵。
单片空间后方交会实习报告范文
单片空间后方交会实习报告范文程序设计实习报告班级：学号：姓名：实习任务：用C或VC++语言实现单片后方交汇的计算。
实习目的：1、深入理解单片空间后方交会的原理，2、在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
3．通过上机调试程序加强动手能力的培养实习原理以单幅影象为基础，从该影象所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素某，Y，Z,ф,ω,κ.共线条件方程如下：某-某0=-f某[a1(某-某)+b1(Y-Y)+c1(Z-Z)]/[a3(某-某)+b3(Y-Y)+c3(Z-Z)]y-y0=-f某[a2(某-某)+b2(Y-Y)+c2(Z-Z)]/[a3(某-某)+b3(Y-Y)+c3(Z-Z)]某,y为像点的像平面坐标；某0,y0,f为影像的外方位元素；某,Y,Z为摄站点的物方空间坐标；某,Y,Z为物方点的物方空间坐标；.实习主要仪器：电脑（创天中文VC++）程序框图：输入原始数据归算像点坐标某,,y计算和确定初值某0,Y0,Z0,φ0,ω0,κ0组成旋转矩阵R计算（某），（y）和l某，ly逐点组成误差方程并法化迭代次数小于限差否？所有点完否？解法方程，求未知数改正数计算改正后外方位元素未知数改正数尺为1:50000；4个地面控制点的地面坐标及其对应像点的像片坐标：结果输出：已知条件:像点坐标某,y:-53.482.21-14.78-76.6310.4664.43153.24-86.15地点坐标某a,Ya,Za:37631.131324.5728.6939101249352386.540426.530319.8757.31-68.9936589.425273.3内方位元素:某0=y0=0f=153.24mm计算结果:旋转矩阵：0.9977090.06753340.00398399-0.06752540.997715-0.00211178-0.00411750.001837920.99999像点坐标位:(单位：mm)-86.15-68.99-53.4182.21-14.78-76.6310.4764.43单位权中误差：7.2602632e-006外方位元素：某=39795.435精度为:1.1254813Y=27476.479精度为:1.2437625Z=7572.6929精度为:0.48380521q=-0.0039840098精度为:0.00018182003w=0.0021117837精度为:0.00015959235k=-0.067576934精度为:7.2440432e-005迭代次数：4Preanykeytocontinue实习总结：通过这次实习，在最初的程序中遇到了很多地困难，但最终还是克服了它。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
中南大学本科生课程设计(实践)任务书、设计报告 (摄影测量与遥感概论)题目空间后方-前方交会学生姓名指导教师邹峥嵘学院地球科学与信息物理学院专业班级测绘0902班学生学号一、实验目的通过对数字影像空间后交前交的程序设计实验，要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。
利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位元素及其精度，然后通过求得的外方位元素求解未知点的地面摄影测量坐标，达到通过摄影测量量测地面地理数据的目的。
二、实验要求➢用C、VB或者Matlab编写空间后方交会-前方交会计算机程序。
➢提交实验报告：程序框图，程序源代码、计算结果及体会。
➢计算结果：地面点坐标、外方位元素及精度。
➢完成时间：2011年11月17日。
三、实验数据214.618-0.231-76.0060.036349.88-0.782-42.201-1.022486.14-1.346-7.706-2.112548.035-79.962-44.438-79.736f=150.000mm，x0=0，y0=0四、实验思路➢利用后方交会得出两张像片各自的外方位元素1)获取已知数据：从摄影资料中插曲像片比例尺、平均航高、内方位元素以及控制点的地面摄影测量坐标及对应的像点坐标。
2)确定未知数的初始值：在竖直摄影的情况下，胶原素的初始值为0，线元素其中Zs=m*f+，Xs=，Ys=。
3)计算旋转矩阵R。
4)逐点计算像点坐标的近似值：利用共线方程。
5)组成误差方程并法化。
6)解求外方位元素。
7)检查计算是否收敛。
➢利用解求出的外方位元素进行前方交会1)用各自像片的角元素计算出左右像片的旋转矩阵R1和R2。
2)根据左右像片的外方位元素计算摄影基线分量Bx,By，Bz。
3)逐点计算像点的空间辅助坐标。
4) 计算投影系数。
5) 计算未知点的地面摄影测量坐标。
6) 重复以上步骤完成所有点的地面坐标的计算。
五、 实验过程➢ 程序流程框图不满足限差则重复计算➢程序中的主要函数设计子函数（矩阵求积multiply，计算函数Resection，矩阵转置transpose，矩阵求逆inMerse1，输出函数shuchu，左片的外方位元素求解函数zuobian。
右片的外方位元素求解函数youbian。
）➢程序源代码#include "stdio.h"#include "math.h"double Xs1,Xs2,Ys1,Ys2,Zs1,Zs2,p01,p02,w01,w02,k01,k02;//求矩阵a的转置矩阵b，a为m行、n列void transpose(double *a, double *b, int m, int n);//矩阵a乘以矩阵b，结果存储在c中，a为m×n大小，b为n×l大小void multiply(double *a, double *b, double *c, int m, int n, int l);//求矩阵a的逆int inMerse1(double *a, int n);//输出m行、n列的矩阵avoid shuchu(double *a, int m, int n);//计算并输出左片的外方位元素void zuobian();//计算并输出右片的外方位元素void youbian();void zuobian(){FILE *fp = NULL;FILE *fp1 = NULL;if((fp=fopen("F:\image.txt","r")) == NULL){printf("Open file error!");return;}if((fp1=fopen("F:\ground.txt","r")) == NULL){printf("Open file error!");return;}//像点坐标和地面点坐标double imagecontrol[4][2]={0.0};double groundcontrol[4][3]={0.0};//摄影比例尺分母double m = 9943;double f=0.15;long i,j,k;for(i=0; i<4; i++){for(j=0; j<2; j++){fscanf(fp, "%lf", &imagecontrol[i][j]);imagecontrol[i][j] /= 1000.0;}for(k=0; k<3; k++){fscanf(fp1, "%lf", &groundcontrol[i][k]);}}fclose(fp);fclose(fp1);//计算外方位元素初始值for( i=0;i<4;i++){Xs1+=groundcontrol[i][0];Ys1+=groundcontrol[i][1];Zs1+=groundcontrol[i][2];}Xs1/=4.0;Ys1/=4.0;Zs1/=4.0;Zs1+=m*f;double R[3][3]={0.0};double L3=0.0,L1=0.0,L2=0.0;double L[8][1]={0.0},x=0.0,y=0.0;double A[8][6]={0.0},A T[6][8]={0.0},ATA[6][6]={0.0},B[6][8]={0.0};double V[6][1]={0.0};int n=0;do{//计算旋转矩阵R[0][0]=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R[0][1]=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R[0][2]=(-1)*sin(p01)*cos(w01);R[1][0]=cos(w01)*sin(k01);R[1][1]=cos(w01)*cos(k01);R[1][2]=(-1)*sin(w01);R[2][0]=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R[2][1]=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R[2][2]=cos(p01)*cos(w01);for(i=0,j=0;j<4;i+=2,j++){//计算像点坐标的近似值L1=R[0][0]*(groundcontrol[j][0]-Xs1)+R[1][0]*(groundcontrol[j][1]-Ys1)+R[2][0]*(ground control[j][2]-Zs1);L2=R[0][1]*(groundcontrol[j][0]-Xs1)+R[1][1]*(groundcontrol[j][1]-Ys1)+R[2][1]*(ground control[j][2]-Zs1);L3=R[0][2]*(groundcontrol[j][0]-Xs1)+R[1][2]*(groundcontrol[j][1]-Ys1)+R[2][2]*(ground control[j][2]-Zs1);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3;//计算常数项L[2*j][0]=imagecontrol[j][0]-x;L[2*j+1][0]=imagecontrol[j][1]-y;//计算系数矩阵A[i][0]=(R[0][0]*f+R[0][2]*imagecontrol[j][0])/L3;A[i][1]=(R[1][0]*f+R[1][2]*imagecontrol[j][0])/L3;A[i][2]=(R[2][0]*f+R[2][2]*imagecontrol[j][0])/L3;A[i][3]=imagecontrol[j][1]*sin(w01)-((imagecontrol[j][0]/f)*(imagecontrol[j][0]*cos(k01)-i magecontrol[j][1]*sin(k01))+f*cos(k01))*cos(w01);A[i][4]=(-1)*f*sin(k01)-(imagecontrol[j][0]/f)*(imagecontrol[j][0]*sin(k01)+imagecontrol[j] [1]*cos(k01));A[i][5]=imagecontrol[j][1];A[i+1][0]=(R[0][1]*f+R[0][2]*imagecontrol[j][1])/L3;A[i+1][1]=(R[1][1]*f+R[1][2]*imagecontrol[j][1])/L3;A[i+1][2]=(R[2][1]*f+R[2][2]*imagecontrol[j][1])/L3;A[i+1][3]=(-1)*imagecontrol[j][0]*sin(w01)-((imagecontrol[j][1]/f)*(imagecontrol[j][0]*cos (k01)-imagecontrol[j][1]*sin(k01))-f*sin(k01))*cos(w01);A[i+1][4]=(-1)*f*cos(k01)-(imagecontrol[j][1]/f)*(imagecontrol[j][0]*sin(k01)+imagecontro l[j][1]*cos(k01));A[i+1][5]=(-1)*imagecontrol[j][0];}transpose(&A[0][0],&AT[0][0],8,6);multiply(&AT[0][0],&A[0][0],&ATA[0][0],6,8,6);inMerse1(*ATA,6);multiply(&ATA[0][0],&AT[0][0],&B[0][0],6,6,8);multiply(&B[0][0],&L[0][0],&V[0][0],6,8,1);Xs1+=V[0][0];Ys1+=V[1][0];Zs1+=V[2][0];p01+=V[3][0];w01+=V[4][0];k01+=V[5][0];n++;}while(fabs(V[3][0])>=0.00001||fabs(V[4][0])>=0.00001||fabs(V[5][0])>=0.00001);R[0][0]=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R[0][1]=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R[0][2]=(-1)*sin(p01)*cos(w01);R[1][0]=cos(w01)*sin(k01);R[1][1]=cos(w01)*cos(k01);R[1][2]=(-1)*sin(w01);R[2][0]=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R[2][1]=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R[2][2]=cos(p01)*cos(w01);//进行未知数的精度评定double A V[8][1],X[8][1],XT[1][8],XTX[1][1],mo,D[6][6],mi[6];multiply(&A[0][0],&V[0][0],&A V[0][0],8,6,1);for(i=0;i<8;i++)X[i][0]=A V[i][0]-L[i][0];transpose(&X[0][0],&XT[0][0],8,1);multiply(&XT[0][0],&X[0][0],&XTX[0][0],1,8,1);mo=XTX[0][0]/2;for(i=0;i<6;i++)for(j=0;j<6;j++)D[i][j]=mo*ATA[i][j];for(i=0;i<6;i++){mi[i]=sqrt(D[i][i]);}printf("左片结果为：\n\n");printf("旋转矩阵R为：\n");shuchu(&R[0][0],3,3);printf("外方为元素为：\n");printf("Xs1=%lf\n",Xs1);printf("Ys1=%lf\n",Ys1);printf("Zs1=%lf\n",Zs1);printf("p01=%lf\n",p01);printf("w01=%lf\n",w01);printf("k01=%lf\n",k01);printf("各外方位元素精度为：\n");printf("m1=%lf\nm2=%lf\nm3=%lf\nm4=%lf\nm5=%lf\nm6=%lf\n",mi[0],mi[1],mi[2],mi[3],mi[ 4],mi[5]);printf("迭代次数为：%d\n\n\n\n",n);fclose(fp);}void youbian(){FILE *fp = NULL;FILE *fp1 = NULL;if((fp=fopen("F:\image2.txt","r")) == NULL){printf("Open file error!");return;}if((fp1=fopen("F:\ground.txt","r")) == NULL){printf("Open file error!");return;}//像点坐标和地面点坐标double imagecontrol[4][2]={0.0};double groundcontrol[4][3]={0.0};//摄影比例尺分母double m = 10177;double f=0.15;long i,j,k;for(i=0; i<4; i++){for(j=0; j<2; j++){fscanf(fp, "%lf", &imagecontrol[i][j]);imagecontrol[i][j] /= 1000.0;}for(k=0; k<3; k++){fscanf(fp1, "%lf", &groundcontrol[i][k]);}}fclose(fp);fclose(fp1);//计算外方位元素初始值for( i=0;i<4;i++){Xs2+=groundcontrol[i][0];Ys2+=groundcontrol[i][1];Zs2+=groundcontrol[i][2];}Xs2/=4.0;Ys2/=4.0;Zs2/=4.0;Zs2+=m*f;double R[3][3]={0.0};double L3=0.0,L1=0.0,L2=0.0;double L[8][1]={0.0},x=0.0,y=0.0;double A[8][6]={0.0},A T[6][8]={0.0},ATA[6][6]={0.0},B[6][8]={0.0};double V[6][1]={0.0};int n=0;do{R[0][0]=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R[0][1]=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R[0][2]=(-1)*sin(p02)*cos(w02);R[1][0]=cos(w02)*sin(k02);R[1][1]=cos(w02)*cos(k02);R[1][2]=(-1)*sin(w02);R[2][0]=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R[2][1]=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R[2][2]=cos(p02)*cos(w02);for(i=0,j=0;j<4;i+=2,j++){//计算像点坐标的近似值L1=R[0][0]*(groundcontrol[j][0]-Xs2)+R[1][0]*(groundcontrol[j][1]-Ys2)+R[2][0]*(ground control[j][2]-Zs2);L2=R[0][1]*(groundcontrol[j][0]-Xs2)+R[1][1]*(groundcontrol[j][1]-Ys2)+R[2][1]*(ground control[j][2]-Zs2);L3=R[0][2]*(groundcontrol[j][0]-Xs2)+R[1][2]*(groundcontrol[j][1]-Ys2)+R[2][2]*(ground control[j][2]-Zs2);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3;//计算常数项矩阵L[2*j][0]=imagecontrol[j][0]-x;L[2*j+1][0]=imagecontrol[j][1]-y;//计算系数矩阵A[i][0]=(R[0][0]*f+R[0][2]*imagecontrol[j][0])/L3;A[i][1]=(R[1][0]*f+R[1][2]*imagecontrol[j][0])/L3;A[i][2]=(R[2][0]*f+R[2][2]*imagecontrol[j][0])/L3;A[i][3]=imagecontrol[j][1]*sin(w02)-((imagecontrol[j][0]/f)*(imagecontrol[j][0]*cos(k02)-i magecontrol[j][1]*sin(k02))+f*cos(k02))*cos(w02);A[i][4]=(-1)*f*sin(k02)-(imagecontrol[j][0]/f)*(imagecontrol[j][0]*sin(k02)+imagecontrol[j] [1]*cos(k02));A[i][5]=imagecontrol[j][1];A[i+1][0]=(R[0][1]*f+R[0][2]*imagecontrol[j][1])/L3;A[i+1][1]=(R[1][1]*f+R[1][2]*imagecontrol[j][1])/L3;A[i+1][2]=(R[2][1]*f+R[2][2]*imagecontrol[j][1])/L3;A[i+1][3]=(-1)*imagecontrol[j][0]*sin(w02)-((imagecontrol[j][1]/f)*(imagecontrol[j][0]*cos (k02)-imagecontrol[j][1]*sin(k02))-f*sin(k02))*cos(w02);A[i+1][4]=(-1)*f*cos(k02)-(imagecontrol[j][1]/f)*(imagecontrol[j][0]*sin(k02)+imagecontro l[j][1]*cos(k02));A[i+1][5]=(-1)*imagecontrol[j][0];}transpose(&A[0][0],&AT[0][0],8,6);multiply(&AT[0][0],&A[0][0],&ATA[0][0],6,8,6);inMerse1(*ATA,6);multiply(&ATA[0][0],&AT[0][0],&B[0][0],6,6,8);multiply(&B[0][0],&L[0][0],&V[0][0],6,8,1);Xs2+=V[0][0];Ys2+=V[1][0];Zs2+=V[2][0];p02+=V[3][0];w02+=V[4][0];k02+=V[5][0];n++;}while(fabs(V[3][0])>=0.00001||fabs(V[4][0])>=0.00001||fabs(V[5][0])>=0.00001);R[0][0]=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R[0][1]=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R[0][2]=(-1)*sin(p02)*cos(w02);R[1][0]=cos(w02)*sin(k02);R[1][1]=cos(w02)*cos(k02);R[1][2]=(-1)*sin(w02);R[2][0]=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R[2][1]=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R[2][2]=cos(p02)*cos(w02);//进行未知数的精度评定double A V[8][1],X[8][1],XT[1][8],XTX[1][1],mo,D[6][6],mi[6];multiply(&A[0][0],&V[0][0],&A V[0][0],8,6,1);for(i=0;i<8;i++)X[i][0]=A V[i][0]-L[i][0];transpose(&X[0][0],&XT[0][0],8,1);multiply(&XT[0][0],&X[0][0],&XTX[0][0],1,8,1);mo=XTX[0][0]/2;for(i=0;i<6;i++)for(j=0;j<6;j++)D[i][j]=mo*ATA[i][j];for(i=0;i<6;i++){mi[i]=sqrt(D[i][i]);}printf("右片结果为：\n\n");printf("旋转矩阵R为：\n");shuchu(&R[0][0],3,3);printf("外方为元素为：\n");printf("Xs2=%lf\n",Xs2);printf("Ys2=%lf\n",Ys2);printf("Zs2=%lf\n",Zs2);printf("p02=%lf\n",p02);printf("w02=%lf\n",w02);printf("k02=%lf\n",k02);printf("各外方位元素精度为：\n");printf("m1=%lf\nm2=%lf\nm3=%lf\nm4=%lf\nm5=%lf\nm6=%lf\n",mi[0],mi[1],mi[2],mi[3],mi[ 4],mi[5]);printf("迭代次数为：%d\n",n);fclose(fp);}//求矩阵a的转置矩阵b，a为m行、n列void transpose(double *a, double *b, int m, int n){int i,j;for(i=0;i<m;i++){for(j=0;j<n;j++){b[j*m+i] = a[i*n+j];}}}//矩阵a乘以矩阵b，结果存储在c中，a为m×n大小，b为n×l大小void multiply(double *a, double *b, double *c, int m, int n, int l){int i,j,k;double t;for(i=0;i<m;i++){for(j=0;j<l;j++){t=0;for(k=0;k<n;k++){t += a[i*n+k]*b[k*l+j];}c[i*l+j]=t;}}}//求矩阵a的逆int inMerse1(double *a, int n){int *is, *js, i, j, k, l, u, M;double d,p;is = new int[n];js = new int[n];for(k=0; k<=n-1; k++){d=0.0;for(i=k; i<=n-1; i++){for(j=k; j<=n-1; j++){l=i*n+j; p=fabs(a[l]);if (p>d) { d=p; is[k]=i; js[k]=j;}}}if(d+1.0==1.0){delete []is;delete []js;printf("Error, not inMerse!\n");return 0;}if(is[k] != k)for(j=0; j<=n-1; j++){u=k*n+j;M=is[k]*n+j;p=a[u];a[u]=a[M];a[M]=p;}if(js[k]!=k)for(i=0; i<=n-1; i++){u=i*n+k;M=i*n+js[k];p=a[u];a[u]=a[M];a[M]=p;}l=k*n+k;a[l]=1.0/a[l];for(j=0; j<=n-1; j++)if(j!=k){u=k*n+j; a[u]=a[u]*a[l];}for(i=0; i<=n-1; i++)if(i!=k)for(j=0; j<=n-1; j++)if(j!=k){u=i*n+j;a[u]=a[u]-a[i*n+k]*a[k*n+j];}for(i=0; i<=n-1; i++)if(i!=k){u=i*n+k;a[u]=-a[u]*a[l];}}for(k=n-1; k>=0; k--){if (js[k]!=k){for(j=0; j<=n-1; j++){u=k*n+j; M=js[k]*n+j;p=a[u]; a[u]=a[M]; a[M]=p;}}if(is[k]!=k){for(i=0; i<=n-1; i++){u=i*n+k; M=i*n+is[k];p=a[u]; a[u]=a[M]; a[M]=p;}}}delete []is;delete []js;return 1;}//输出m行、n列的矩阵avoid shuchu(double *a, int m, int n){int i,j;for(i=0;i<m;i++){for(j=0;j<n;j++){printf("%9lf ",a[i*n+j]);}printf("\n");}}void main(){zuobian();youbian();double a11,a12,a13,a21,a22,a23,b11,b12,b13,b21,b22,b23,c11,c12,c13,c21,c22,c23;double x1,y1,x2,y2,X1,Y1,Z1,X2,Y2,Z2,N1,N2,Xt,Yt,Zt,Bx,By,Bz;double f=0.15;scanf("x1=%lf,y1=%lf,x2=%lf,y2=%lf",&x1,&y1,&x2,&y2);//scanf("y1=%lf",&y1);//printf("x2=");//scanf("x2=%lf",&x2);//旋转矩阵a11=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);a12=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);a13=(-1)*sin(p01)*cos(w01);a21=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);a22=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);a23=(-1)*sin(p02)*cos(w02);b11=cos(w01)*sin(k01);b12=cos(w01)*cos(k01);b13=(-1)*sin(w01);b21=cos(w02)*sin(k02);b22=cos(w02)*cos(k02);b23=(-1)*sin(w02);c11=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);c12=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);c13=cos(p01)*cos(w01);c21=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);c22=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);c23=cos(p02)*cos(w02);//计算各待定点的像空间辅助坐标X1=a11*x1/1000+a12*y1/1000-a13*f;Y1=b11*x1/1000+b12*y1/1000-b13*f;Z1=c11*x1/1000+c12*y1/1000-c13*f;X2=a21*x2/1000+a22*y2/1000-a23*f;Y2=b21*x2/1000+b22*y2/1000-b23*f;Z2=c21*x2/1000+c22*y2/1000-c23*f;//计算摄影基线分量Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;//计算点投影系数N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1);N2=(Bx*Z1-Bz*X1)/(X1*Z2-X2*Z1);//计算待定点的地面摄影测量坐标Xt=((N1*X1+Xs1)+(N2*X2+Xs2))/2;Yt=((N1*Y1+Ys1)+(N2*Y2+Ys2))/2;Zt=((N1*Z1+Zs1)+(N2*Z2+Zs2))/2-1500;printf("Xt=%lf\nYt=%lf\nZt=%lf\n",Xt,Yt,Zt);}➢实验结果左片右片地面摄影测量坐标点号x y x y X Y Z GCP116.01279.963-73.9378.7065083.2055852.099527.925 GCP288.5681.134-5.25278.1845780.025906.365571.549 GCP313.362-79.37-79.122-78.8795210.8794258.446461.81 GCP482.24-80.027-9.887-80.0895909.2644314.283455.484 151.75880.555-39.95378.4635431.489 5879.367 549.723 214.618-0.231-76.0060.0365147.388 5055.549 485.001 349.88-0.782-42.201-1.0225495.786 5082.726 506.677 486.14-1.346-7.706-2.1125844.173 5109.861 528.420 548.035-79.962-44.438-79.7365559.940 4268.185 463.523六、实验总结从编写思路及步骤，然后是梳理框架，要用到的子函数，子函数参数，怎么编写，有不懂的东西网上参考资料；再是自己动手编程，发现编程要求很细腻，出不得一点差错，程序编写出来后，有很多的错误，要求逐一进行改正，修改；最后才得以运行。