单向后方交会实验报告
单像空间后方交会实习报告
摄影测量学单像空间后方交会实习报告一、实习目的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弧度时，迭代结束。
空间后方交会程序设计
单片空间后方交会程序设计实验报告专业：测绘工程班级：081姓名：张发伟学号：09[一]、实习任务：用C或VC++语言实现单片后方交汇的计算。
[二]、实习目的：1.深入理解单片空间后方交会的原理；2.在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程；3．通过上机调试程序加强动手能力的培养。
[三]、实习原理:以单幅影象为基础，从该影象所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素Xs，Ys，Zs,ф,ω,κ.共线条件方程如下：x-x0=-f*[a1(X-Xs)+b1(Y-Ys)+c1(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] y-y0=-f*[a2(X-Xs)+b2(Y-Ys)+c2(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)]x,y为像点的像平面坐标；x0,y0,f为影像的外方位元素；Xs,Ys,Zs为摄站点的物方空间坐标；X,Y,Z为物方点的物方空间坐标；[四]、程序框图：输入原始数据归算像点坐标x,,y计算和确定初值X s 0, Y s 0, Z s 0, φ0,ω0,κ0组成旋转矩阵R计算（x ），（y ）和l x ，l y逐点组成误差方程并法化所有点完否？ 解法方程，求未知数改正数 计算改正后外方位元素未知数改正数<限差否？整理并输出计算结果正常输出迭代次数小于限差否？ 输出中间结果和出错信息非正常结束[五]、实验数据：4个地面控制点的地面坐标及其对应像点的像片坐标：[六]、结果输出：已知条件:像点坐标x,y:-53.4 82.21-14.78 -76.6310.46 64.43153.24 -86.15地点坐标Xa,Ya,Za:37631.1 31324.5 728.6939101 24935 2386.540426.5 30319.8 757.31-68.99 36589.4 25273.3摄影比例尺：1:50000内方位元素:x0=y0=0 f=153.24mm计算结果:旋转矩阵：0.997709 0.0675334 0.00398399-0.0675254 0.997715 -0.00211178-0.0041175 0.00183792 0.99999像点坐标位:(单位：mm)-86.15 -68.99-53.41 82.21-14.78 -76.6310.47 64.43外方位元素：Xs=39795.435 精度为:1.1254813Ys=27476.479 精度为:1.2437625Zs=7572.6929 精度为:0.48380521q=-0.0039840098 精度为:0.00018182003w=0.0021117837 精度为:0.00015959235k=-0.067576934 精度为:7.2440432e-005迭代次数：4Press any key to continue[七]、实习总结：掌握运用程序设计实现单片后方交汇的计算，了解在多余观测情况下，用最小二乘平差方法计算相片外方位元素，及其相关精度。
摄影测量学后方交会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．实验结果四．实验总结此次实验让我深入了解单像空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。
后方交会实验报告
摄影测量学实验报告实验名称：空间后方交会班级：测绘工程12-1一．实验目的掌握运用空间后方交会求解外方元素的过程。
学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用已知公式和计算机编程语言实现空间后方交会的过程，完成外方位元素的求解。
二．实验工具个人电脑，visual basic6.0三．实验源代码Private Sub Command1_Click()Dim c(),kx(),ky(),l(),aa(),a11(),a12(),a13(),a14(),a15(),a16(), a21(),a22(),a23(),a24(),a25()Dim a26(),v(),mm(5)nn=Val(Text2.Text)nnt=nn-1ReDim c(nnt,5),kx(nnt),ky(nnt),l(2*nnt+1,0),aa(2*nnt+1,5), a11(nnt),a12(nnt),a13(nnt),a14(nnt),a15(nnt),a16(nnt),a21(nnt), a22(nnt),a23(nnt),a24(nnt),a25(nnt)ReDim a26(nnt)f=Val(Text1.Text)fei=0#w=0#k=0#If Label2.Caption=""ThenMsgBox"尚未读取文本",,"提示信息"Elsea=Split(Label2.Caption,vbCrLf)For i=0To nnthh=Split(a(i),"")For j=0To5c(i,j)=Val(hh(j))Next jNext im1=(c(0,1)-c(1,1))^2+(c(0,2)-c(1,2))^2 m2=(c(0,3)-c(1,3))^2+(c(0,4)-c(1,4))^2 m=Sqr(m2/m1)zs=m*fxs=0ys=0For i=0To nntxs=xs+c(i,3)ys=xs+c(i,4)Next ixs=xs/nnys=ys/nnDoa1=Cos(fei)*Cos(k)-Sin(fei)*Sin(w)*Sin(k)a2=-Cos(fei)*Sin(k)-Sin(fei)*Sin(w)*Cos(k)a3=-Sin(fei)*Cos(w)b1=Cos(w)*Sin(k)b2=Cos(w)*Cos(k)b3=-Sin(w)c1=Sin(fei)*Cos(k)+Cos(fei)*Sin(w)*Sin(k)c2=-Sin(fei)*Sin(k)+Cos(fei)*Sin(w)*Cos(k)c3=Cos(fei)*Cos(w)For i=0To nntkx(i)=-f*(a1*(c(i,3)-xs)+b1*(c(i,4)-ys)+c1*(c(i,5)-zs)) /(a3*(c(i,3)-xs)+b3*(c(i,4)-ys)+c3*(c(i,5)-zs))ky(i)=-f*(a2*(c(i,3)-xs)+b2*(c(i,4)-ys)+c2*(c(i,5)-zs)) /(a3*(c(i,3)-xs)+b3*(c(i,4)-ys)+c3*(c(i,5)-zs))Next iFor i=0To nntl(i*2,0)=c(i,1)-kx(i)l(i*2+1,0)=c(i,2)-ky(i)Next iFor i=0To nntzb=(a3*(c(i,3)-xs)+b3*(c(i,4)-ys)+c3*(c(i,5)-zs))a11(i)=(a1*f+a3*c(i,1))/zba12(i)=(b1*f+b3*c(i,1))/zba13(i)=(c1*f+c3*c(i,1))/zba21(i)=(a2*f+a3*c(i,2))/zba22(i)=(b2*f+b3*c(i,2))/zba23(i)=(c2*f+c3*c(i,2))/zba14(i)=c(i,2)*Sin(w)-(c(i,1)*(c(i,1)*Cos(k)-c(i,2)*Sin(k)) /f+f*Cos(k))*Cos(w)a15(i)=-f*Sin(k)-c(i,1)/f*(c(i,1)*Sin(k)+c(i,2)*Cos(k))a16(i)=c(i,2)a24(i)=-c(i,1)*Sin(w)-(c(i,2)*(c(i,1)*Cos(k)-c(i,2)*Sin(k)) /f-f*Sin(k))*Cos(w)a25(i)=-f*Cos(k)-c(i,2)/f*(c(i,1)*Sin(k)+c(i,2)*Cos(k)) a26(i)=-c(i,1)Next iFor i=0To nntaa(2*i,0)=a11(i)aa(2*i,1)=a12(i)aa(2*i,2)=a13(i)aa(2*i,3)=a14(i)aa(2*i,4)=a15(i)aa(2*i,5)=a16(i)aa(2*i+1,0)=a21(i)aa(2*i+1,1)=a22(i)aa(2*i+1,2)=a23(i)aa(2*i+1,3)=a24(i)aa(2*i+1,4)=a25(i)aa(2*i+1,5)=a26(i)Next iReDim at(5,2*nnt+1),ata(5,5),nata(5,5),nataat(5,2*nnt+1), detx(5,0)at=zhuanzhi(aa)ata=cheng(at,aa)nata=qiuni(ata)nataat=cheng(nata,at)detx=cheng(nataat,l)xs=xs+detx(0,0)ys=ys+detx(1,0)zs=zs+detx(2,0)fei=fei+detx(3,0)w=w+detx(4,0)k=k+detx(5,0)Loop Until Abs(detx(0,0))<0.000001And Abs(detx(1,0))< 0.000001And Abs(detx(2,0))<0.000001And Abs(detx(3,0))<0.000001And Abs(detx(4,0))<0.000001And Abs(detx(5,0))< 0.000001ReDim v(2*nnt+1,0),ax(2*nnt+1,0)ax=cheng(aa,detx)v=jian(ax,l)vv=0For i=0To2*nnt+1vv=v(i,0)^2+vvNext imo=Sqr(vv/(2*nn-6))For i=0To5mm(i)=mo*Sqr(nata(i,i))Next itxtshow.Text=txtshow.Text&"外方元素为中误差为"&vbCrLf&"Xs"&xs&" "&mm(0)&vbCrLf&"Ys"&ys&" "&mm(1)&vbCrLf&"Zs"&zs&" "&mm(2)&vbCrLf&"Φ"&fei&"" &mm(3)&vbCrLf&"W"&w&"" &mm(4)&vbCrLf&"K"&k&""& mm(5)&vbCrLf&vbCrLftxtshow.Text=txtshow.Text&"报告时间："&Year(Date)&Format(Month(Date),"00")& Format(Day(Date),"00")&vbCrLftxtshow.Text=txtshow.Text&"审批人签字："End IfEnd SubPrivate Sub Command2_Click()CommonDialog1.Filter="文档|*.txt"CommonDialog1.ShowOpenLabel1.Caption=CommonDialog1.FileNameEnd SubPrivate Sub Command3_Click()Label2.Caption=""Dim mylineIf Label1.Caption<>""ThenOpen Label1.Caption For Input As#1Do While Not EOF(1)Line Input#1,mylineLabel2.Caption=Label2.Caption+myline+vbCrLfLoopClose#1End SubPrivate Sub Command4_Click()CommonDialog1.DialogTitle="保存文件" CommonDialog1.Filter="文档|*.doc" CommonDialog1.Action=2If CommonDialog1.FileName<>""ThenOpen CommonDialog1.FileName For Output As#1 Print#1,txtshow.TextClose#1End IfEnd SubPrivate Sub Command5_Click()txtshow.Text=""End SubPrivate Sub Command6_Click()If MsgBox("确认退出程序？",68,"提示信息")=6Then EndEnd IfEnd SubPublic Function cheng(a(),b())ReDim s(UBound(a,1),UBound(b,2))If UBound(a,2)<>UBound(b,1)Then MsgBox"不符合乘法要求"EndElseFor i=0To UBound(s,1)For j=0To UBound(s,2)s(i,j)=0For Q=0To UBound(a,2)s(i,j)=Val(s(i,j))+Val(a(i,Q))*Val(b(Q,j)) Next QNext jNext icheng=sEnd IfEnd FunctionPublic Function zhuanzhi(a())Dim sReDim s(UBound(a,2),UBound(a,1))For i=0To UBound(a,1)For j=0To UBound(a,2)s(j,i)=a(i,j)Next jNext izhuanzhi=sEnd FunctionPublic Function qiuni(m())Dim s()If UBound(m,1)<>UBound(m,2)Then MsgBox("求逆错误")EndElsen=UBound(m,1)Dim e(),m1,m2,m3ReDim e(n,2*n+1)ReDim s(n,2*n+1)For i=0To nFor j=0To ne(i,j)=m(i,j)Next jNext iFor i=0To nFor j=n+1To2*n+1e(i,j)=0Next jNext iFor i=0To ne(i,n+i+1)=1 Next iFor t=0To nIf e(t,t)=0Then For i=t+1To nFor j=0To2*n+1 s(i,j)=e(i,j)e(i,j)=e(t,j)e(t,j)=s(t,j)Next jIf e(t,t)<>0Then GoTo daima1End IfNext iIf e(t,t)=0Then MsgBox("求逆错误") EndGoTo lastlineEnd Ifdaima1:m1=e(t,t)For j=0To2*n+1e(t,j)=e(t,j)/m1Next jFor i=t+1To nm2=e(i,t)For j=0To2*n+1e(i,j)=e(i,j)-m2*e(t,j) Next jNext iNext tFor t=0To n-1For i=t+1To nm3=e(t,i)For j=0To2*n+1e(t,j)=e(t,j)-m3*e(i,j) Next jNext iNext tReDim c(n,n)For i=0To nFor j=0To nc(i,j)=Round(e(i,j+n+1),2)Next jNext iqiuni=cEnd Iflastline:End FunctionPublic Function jian(a(),b())Dim s,tReDim s(UBound(a,1),UBound(a,2))If UBound(a,1)<>UBound(b,1)Or UBound(a,2)<>UBound(b,2) ThenMsgBox"不符合减法要求"EndElseFor i=0To UBound(a,1)For j=0To UBound(a,2)s(i,j)=Val(a(i,j))-Val(b(i,j))Next jNext ijian=sEnd IfEnd FunctionPrivate Sub Form_Load()Form1.Image1.Top=0Form1.Image1.Left=0Form1.Image1.Width=Form1.Width Form1.Image1.Height=Form1.Height End SubPrivate Sub Form_Resize()Form1.Image1.Top=0Form1.Image1.Left=0Form1.Image1.Width=Form1.Width Form1.Image1.Height=Form1.Height End Sub四．实验框图获取已知数据量测控制点的像点坐标确定未知数初始值五．实验数据1-0.08615-0.0689936589.4125273.322195.17 2-0.053400.0822137631.0831324.51728.69 3-0.01478-0.0766339100.9724934.982386.5 40.010460.0644340426.5430319.81757.31六．实验截图七．实验心得此次实验让我更加了解空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。
单像空间后方交会实习报告
框图。
3
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，系电，通力根1保过据护管生高线产中敷工资设艺料技高试术中卷0资不配料仅置试可技卷以术要解是求决指，吊机对顶组电层在气配进设置行备不继进规电行范保空高护载中高与资中带料资负试料荷卷试下问卷高题总中2体2资，配料而置试且时卷可，调保需控障要试各在验类最；管大对路限设习度备题内进到来行位确调。保整在机使管组其路高在敷中正设资常过料工程试况1卷中下安，与全要过，加度并强工且看作尽护下可1都关能可于地以管缩正路小常高故工中障作资高；料中对试资于卷料继连试电接卷保管破护口坏进处范行理围整高，核中或对资者定料对值试某，卷些审弯异核扁常与度高校固中对定资图盒料纸位试，置卷编.工保写况护复进层杂行防设自腐备动跨与处接装理地置，线高尤弯中其曲资要半料避径试免标卷错高调误等试高，方中要案资求，料技编试术写5、卷交重电保底要气护。设设装管备备置线4高、调动敷中电试作设资气高，技料课中并3术试、件资且中卷管中料拒包试路调试绝含验敷试卷动线方设技作槽案技术，、以术来管及避架系免等统不多启必项动要方高式案中，；资为对料解整试决套卷高启突中动然语过停文程机电中。气高因课中此件资，中料电管试力壁卷高薄电中、气资接设料口备试不进卷严行保等调护问试装题工置，作调合并试理且技利进术用行，管过要线关求敷运电设行力技高保术中护。资装线料置缆试做敷卷到设技准原术确则指灵：导活在。。分对对线于于盒调差处试动，过保当程护不中装同高置电中高压资中回料资路试料交卷试叉技卷时术调，问试应题技采，术用作是金为指属调发隔试电板人机进员一行，变隔需压开要器处在组理事在；前发同掌生一握内线图部槽 纸故内资障，料时强、，电设需回备要路制进须造行同厂外时家部切出电断具源习高高题中中电资资源料料，试试线卷卷缆试切敷验除设报从完告而毕与采，相用要关高进技中行术资检资料查料试和，卷检并主测且要处了保理解护。现装场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
后方交会实验报告
摄影测量学实验报告实验名称：空间后方交会班级：测绘工程12-1一．实验目的掌握运用空间后方交会求解外方元素的过程。
学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用已知公式和计算机编程语言实现空间后方交会的过程，完成外方位元素的求解。
二．实验工具个人电脑，visual basic 6.0三．实验源代码Private Sub Command1_Click()Dim c(), kx(), ky(), l(), aa(), a11(), a12(), a13(), a14(), a15(), a16(), a21(), a22(), a23(), a24(), a25()Dim a26(), v(), mm(5)nn = Val(Text2.Text)nnt = nn - 1ReDim c(nnt, 5), kx(nnt), ky(nnt), l(2 * nnt + 1, 0), aa(2 * nnt + 1, 5), a11(nnt), a12(nnt), a13(nnt), a14(nnt), a15(nnt), a16(nnt), a21(nnt), a22(nnt), a23(nnt), a24(nnt), a25(nnt)ReDim a26(nnt)f = Val(Text1.Text)fei = 0#w = 0#k = 0#If Label2.Caption = "" ThenMsgBox "尚未读取文本", , "提示信息"Elsea = Split(Label2.Caption, vbCrLf)For i = 0 To nnthh = Split(a(i), " ")For j = 0 To 5c(i, j) = Val(hh(j))Next jNext im1 = (c(0, 1) - c(1, 1)) ^ 2 + (c(0, 2) - c(1, 2)) ^ 2 m2 = (c(0, 3) - c(1, 3)) ^ 2 + (c(0, 4) - c(1, 4)) ^ 2 m = Sqr(m2 / m1)zs = m * fxs = 0ys = 0For i = 0 To nntxs = xs + c(i, 3)ys = xs + c(i, 4)Next ixs = xs / nnys = ys / nnDoa1 = Cos(fei) * Cos(k) - Sin(fei) * Sin(w) * Sin(k)a2 = -Cos(fei) * Sin(k) - Sin(fei) * Sin(w) * Cos(k)a3 = -Sin(fei) * Cos(w)b1 = Cos(w) * Sin(k)b2 = Cos(w) * Cos(k)b3 = -Sin(w)c1 = Sin(fei) * Cos(k) + Cos(fei) * Sin(w) * Sin(k)c2 = -Sin(fei) * Sin(k) + Cos(fei) * Sin(w) * Cos(k)c3 = Cos(fei) * Cos(w)For i = 0 To nntkx(i) = -f * (a1 * (c(i, 3) - xs) + b1 * (c(i, 4) - ys) + c1 * (c(i, 5) - zs)) / (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))ky(i) = -f * (a2 * (c(i, 3) - xs) + b2 * (c(i, 4) - ys) + c2 * (c(i, 5) - zs)) / (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))Next iFor i = 0 To nntl(i * 2, 0) = c(i, 1) - kx(i)l(i * 2 + 1, 0) = c(i, 2) - ky(i)Next iFor i = 0 To nntzb = (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))a11(i) = (a1 * f + a3 * c(i, 1)) / zba12(i) = (b1 * f + b3 * c(i, 1)) / zba13(i) = (c1 * f + c3 * c(i, 1)) / zba21(i) = (a2 * f + a3 * c(i, 2)) / zba22(i) = (b2 * f + b3 * c(i, 2)) / zba23(i) = (c2 * f + c3 * c(i, 2)) / zba14(i) = c(i, 2) * Sin(w) - (c(i, 1) * (c(i, 1) * Cos(k) - c(i, 2) * Sin(k)) / f + f * Cos(k)) * Cos(w)a15(i) = -f * Sin(k) - c(i, 1) / f * (c(i, 1) * Sin(k) + c(i, 2) * Cos(k))a16(i) = c(i, 2)a24(i) = -c(i, 1) * Sin(w) - (c(i, 2) * (c(i, 1) * Cos(k) - c(i, 2) * Sin(k)) / f - f * Sin(k)) * Cos(w)a25(i) = -f * Cos(k) - c(i, 2) / f * (c(i, 1) * Sin(k) + c(i, 2) * Cos(k))a26(i) = -c(i, 1)Next iFor i = 0 To nntaa(2 * i, 0) = a11(i)aa(2 * i, 1) = a12(i)aa(2 * i, 2) = a13(i)aa(2 * i, 3) = a14(i)aa(2 * i, 4) = a15(i)aa(2 * i, 5) = a16(i)aa(2 * i + 1, 0) = a21(i)aa(2 * i + 1, 1) = a22(i)aa(2 * i + 1, 2) = a23(i)aa(2 * i + 1, 3) = a24(i)aa(2 * i + 1, 4) = a25(i)aa(2 * i + 1, 5) = a26(i)Next iReDim at(5, 2 * nnt + 1), ata(5, 5), nata(5, 5), nataat(5, 2 * nnt + 1), detx(5, 0)at = zhuanzhi(aa)ata = cheng(at, aa)nata = qiuni(ata)nataat = cheng(nata, at)detx = cheng(nataat, l)xs = xs + detx(0, 0)ys = ys + detx(1, 0)zs = zs + detx(2, 0)fei = fei + detx(3, 0)w = w + detx(4, 0)k = k + detx(5, 0)Loop Until Abs(detx(0, 0)) < 0.000001 And Abs(detx(1, 0)) < 0.000001 And Abs(detx(2, 0)) < 0.000001 And Abs(detx(3, 0)) <0.000001 And Abs(detx(4, 0)) < 0.000001 And Abs(detx(5, 0)) < 0.000001ReDim v(2 * nnt + 1, 0), ax(2 * nnt + 1, 0)ax = cheng(aa, detx)v = jian(ax, l)vv = 0For i = 0 To 2 * nnt + 1vv = v(i, 0) ^ 2 + vvNext imo = Sqr(vv / (2 * nn - 6))For i = 0 To 5mm(i) = mo * Sqr(nata(i, i))Next itxtshow.Text = txtshow.Text & "外方元素为中误差为" & vbCrLf & "Xs " & xs & " " & mm(0) & vbCrLf & "Ys " & ys & " " & mm(1) & vbCrLf & "Zs " & zs & " " & mm(2) & vbCrLf & "Φ" & fei & " " & mm(3) & vbCrLf & "W " & w & " " & mm(4) & vbCrLf & "K " & k & " " & mm(5) & vbCrLf & vbCrLftxtshow.Text = txtshow.Text & " 报告时间：" & Year(Date) & Format(Month(Date), "00") & Format(Day(Date), "00") & vbCrLftxtshow.Text = txtshow.Text & " 审批人签字："End IfEnd SubPrivate Sub Command2_Click()CommonDialog1.Filter = "文档|*.txt"CommonDialog1.ShowOpenLabel1.Caption = CommonDialog1.FileNameEnd SubPrivate Sub Command3_Click()Label2.Caption = ""Dim mylineIf Label1.Caption <> "" ThenOpen Label1.Caption For Input As #1Do While Not EOF(1)Line Input #1, mylineLabel2.Caption = Label2.Caption + myline + vbCrLfLoopClose #1End SubPrivate Sub Command4_Click()CommonDialog1.DialogTitle = "保存文件" CommonDialog1.Filter = "文档|*.doc" CommonDialog1.Action = 2If CommonDialog1.FileName <> "" ThenOpen CommonDialog1.FileName For Output As #1Print #1, txtshow.TextClose #1End IfEnd SubPrivate Sub Command5_Click()txtshow.Text = ""End SubPrivate Sub Command6_Click()If MsgBox("确认退出程序？", 68, "提示信息") = 6 Then EndEnd IfEnd SubPublic Function cheng(a(), b())ReDim s(UBound(a, 1), UBound(b, 2))If UBound(a, 2) <> UBound(b, 1) Then MsgBox "不符合乘法要求"EndElseFor i = 0 To UBound(s, 1)For j = 0 To UBound(s, 2)s(i, j) = 0For Q = 0 To UBound(a, 2)s(i, j) = Val(s(i, j)) + Val(a(i, Q)) * Val(b(Q, j)) Next QNext jNext icheng = sEnd IfEnd FunctionPublic Function zhuanzhi(a())Dim sReDim s(UBound(a, 2), UBound(a, 1))For i = 0 To UBound(a, 1)For j = 0 To UBound(a, 2)s(j, i) = a(i, j)Next jNext izhuanzhi = sEnd FunctionPublic Function qiuni(m())Dim s()If UBound(m, 1) <> UBound(m, 2) Then MsgBox ("求逆错误")EndElsen = UBound(m, 1)Dim e(), m1, m2, m3ReDim e(n, 2 * n + 1)ReDim s(n, 2 * n + 1)For i = 0 To nFor j = 0 To ne(i, j) = m(i, j)Next jNext iFor i = 0 To nFor j = n + 1 To 2 * n + 1e(i, j) = 0Next jNext iFor i = 0 To ne(i, n + i + 1) = 1 Next iFor t = 0 To nIf e(t, t) = 0 ThenFor i = t + 1 To nFor j = 0 To 2 * n + 1 s(i, j) = e(i, j)e(i, j) = e(t, j)e(t, j) = s(t, j)Next jIf e(t, t) <> 0 Then GoTo daima1End IfNext iIf e(t, t) = 0 Then MsgBox ("求逆错误") EndGoTo lastlineEnd Ifdaima1:m1 = e(t, t)For j = 0 To 2 * n + 1e(t, j) = e(t, j) / m1Next jFor i = t + 1 To nm2 = e(i, t)For j = 0 To 2 * n + 1e(i, j) = e(i, j) - m2 * e(t, j) Next jNext iNext tFor t = 0 To n - 1For i = t + 1 To nm3 = e(t, i)For j = 0 To 2 * n + 1e(t, j) = e(t, j) - m3 * e(i, j) Next jNext iNext tReDim c(n, n)For i = 0 To nFor j = 0 To nc(i, j) = Round(e(i, j + n + 1), 2)Next jNext iqiuni = cEnd Iflastline:End FunctionPublic Function jian(a(), b())Dim s, tReDim s(UBound(a, 1), UBound(a, 2))If UBound(a, 1) <> UBound(b, 1) Or UBound(a, 2) <> UBound(b, 2) ThenMsgBox "不符合减法要求"EndElseFor i = 0 To UBound(a, 1)For j = 0 To UBound(a, 2)s(i, j) = Val(a(i, j)) - Val(b(i, j))Next jNext ijian = sEnd IfEnd FunctionPrivate Sub Form_Load()Form1.Image1.Top = 0Form1.Image1.Left = 0Form1.Image1.Width = Form1.Width Form1.Image1.Height = Form1.Height End SubPrivate Sub Form_Resize()Form1.Image1.Top = 0Form1.Image1.Left = 0Form1.Image1.Width = Form1.Width Form1.Image1.Height = Form1.Height End Sub四．实验框图获取已知数据量测控制点的像点坐标确定未知数初始值五．实验数据1 -0.08615 -0.06899 36589.41 25273.32 2195.172 -0.05340 0.08221 37631.08 31324.51 728.693 -0.01478 -0.07663 39100.97 24934.98 2386.54 0.01046 0.06443 40426.54 30319.81 757.31 六．实验截图七．实验心得此次实验让我更加了解空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。
试验一 单片空间后方交会程序设计
三、存储数据的变量及数组
针对空间后方交会计算中数据的特点，需定义相 关的变量和数组来存储数据。 存储已知数据的数组：x[]，y[]分别用来存储控制点 所对应的像点坐标；X[]，Y[]，Z[]分别用来存储地 面控制点坐标； 存储旋转矩阵数组：a[]，b[], c []分别用来存储旋 转矩阵的九个方向余弦值 误差方程系数数组： A[]数组用来存储误差方程系数 误差方程常数项数组：l []数组用来存储误差方程常 数项
解法方程，求未知数改正数 计算改正后的外方位元素 否 未知数改正数<限差否？ 整理并输出计算结果 正常结束
输出中间结果 和出错信息 非正常结束
已知航摄仪的内方位元素：f=153.24mm， x0=y0=0.0mm，摄影比例尺为1:50000； 4个地面控制 点的地面坐标及其对应像点的像片坐标如下表：
存储外方位元素改正值的数组：H[]存储某次迭 代计算出的参数改正值 相关变量：Xs0, Ys0, Zs0分别用来存储三个外 方位线元素；t，w，k分别用来存储三个外方 位角元素；f用来存储主距；m存储摄影比例尺 分母
用程序设计语言编写一个完整的单片空间后方交会程序，通 过对提供的试验数据进行计算，输出像片的外方位元素。
二、迭代计算的思想
迭代计算一般需要事先知道参数的初值，初 值选的不好，可能迭代计算收敛速度较慢，甚 至不收敛。
每计算出参数的改正值后，需要把改正值与 限差（迭代收敛的条件）进行比较，直到改正 值小于限差，迭代计算结束。
试验一 单片空间后方交会程序设计
重点： 1.理解单片空间后方交会计算流程 2.掌握迭代计算的思想 3.数据读入与存储方法
输入原始数据 归算像点坐标 x, y
一、空间后方交会 计算程序框图
是 迭 代 次 Ys0 , Z s0 , 0 , 0 , 0 组成旋转矩阵R 计算 ( x), ( y )和 l x , l y 直接生成误差方程B和L 否 所有点完否？
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
班级：测绘一班 学号:20133279日期：2015426Southw^ JIaotong UnjiverSity、计算原理.二、算法流程.三、源程序.四、计算结果.五、结果分析.六、心得体会. 目录131313、计算原理已知条件摄影机主距f=153.24mm, xO=O, yO=O,像片比例尺为1:40000,有四对点的像点坐标与相应的地面坐标如下表。
二、算法流程（1）获取已知数据。
从航摄资料中差取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影坐标。
（2）量测控制点的像点坐标并作系统误差改正。
（3）确定未知数的初始值。
在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即x S 上，Y s0n丄Z0n，Z s mf（4）用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵R31 cos cos sin sin sina2cos sin sin sin cos33 sin cosb1 cos sinb2 cos cosb3 sinC1 sin cos cos sin sin'c sin sin cos sin cosC3 cos cos（5）逐点计算像点坐标的近似值。
禾I」用未知数的近似值和控制点的地面坐标; 带入共线方程式，逐点近似像点坐标的近似值（（6）(7) X）、（y）。
逐点计算误差方程式的系数和常数项，组成误差方程式。
计算法方程的系数矩阵A A和常数项A L L，组成法方程式。
(8) 解法方程，求得外方位元素的改正数dX s、dY s、dZ s、d、d、d 。
(9)用前次迭代取得的近似值，加本次迭代的改正数,xr xL dX S\Y S K Y S K1 dY S K,Z K 计算外方位元素的新值。
z K1dz（K K 1 . K K K 1 . K Kd , d ,（10）将求得的外方位元素改正数与规定的限差比较，负责用新的近似值重复（4）-（9），直到满足要求为止。
若小于限差，则迭代结束。
用共线方程进行空间后方交会的程序框如图所示。
输入原始数据像点坐标计算，系统误差正 确定外方位因素初始值组成旋转矩阵R所有像点完否是C解法方程，求外方位元素改正数计算改正后的外方位元素 外方位元素改正数是否小于限差 输出计算成果，计算并结束、源程序#in elude ”ostream" using n ames pace std;#in elude "fstream" #in clude <ioma nip>const int n=6;void in verse(double c[n ][ n]);temp late<t ypen ame Tl,t ypen ame T2>void transp ose (TI*mat1,T2*mat2,i nt a,i nt b);temp latevty pen ame TI,t ypen ame T2>void multi(T1*mat1,T2*mat2,T2*result, i nt a, i nt b, i nt c); int mai n() {doublex[4][2]={-0.08615,-0.06899,-0.05340,0.08221 ，-0.01478，-0.07663,0.01046,0.06443};doubleX[4][3]={36589.41,25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98, 2386.50,40426.54,30319.81,757.31};int i,j,num=0;//num 为迭代次数double X0[6]={0};//设定未知数(XS,YS,ZS,三个角度)初始值 double f=0.15324;// 摄影机主距 f=153.24mm结束并显示错误信息逐点组成误差方程式并法化否否L迭代次数小于ndouble a=1/40000.0;// 像片比例尺为 1:40000 double R[3][3]={0};// 初始化旋转知阵 R double approx_ x[8]={0};// 用于存放像点估计值 double A[8][6] 二 {0};// 定义了一个系数阵 double AT[6][8] 二 {0};// 定义了 A 的转置知阵 double L[8]={0}:// 定义常数项 const double pi=3.1415926535897932; double Asum[6][6]={0}; double result2[6]={0}; double resultl[6][8]={0}; double sumXYZ[3]={0}; cout.precision(5); cout<< ”已知像点坐标为 :\n ”; for(i=0;i<4;i++)for(j=0;j<2;j++){cout<<fixed: if (j==0)cout<< ”x ”<<i+l<< ”=”<<setw (8)<<x[i][j]<<{cout<< ”Z ”;cout<<i+1;cout<<}}cout<<endl; for(j=0;j<3;j++)for(i=0;i<4;i++) sumXYZ[j]+=X[i][j];for(i=0;i<2;i++)X0[i]=sumXYZ[i]/4;//X0,Y0 初始化X0[i]=1/a*f+sumXYZ[2]/4. 0;// 对 Z0 进行初始化 do{}else{cout<< ”y"<<i+1<< ”=”<<setw(6)<<x[i][j]<<en dl;}}cout<< ”己知地面四个点的坐标为 for(i=0;i<4;i++) for(j=0;j<3;j++) {if (j==0){cout<< ”X";cout<<i+1;cout<<}else if(j==1) {cout<< ”Y";cout<<i+1;cout<<}else:\n ”;<<X[i][j]<< ” ”<<X[i][j)<< ” ””;cout<<X[i][j]<<endl;R[0][0]=cos(X0[3])*cos(X0[5])-sin(X0[3])*sin(X0[4])*sin(X0[5]);R[0][I]=-cos(X0[3])*sin(X0[5])-sin(X0[3])*sin(X0[4])*cos(X0[5]);R[0][2]=-sin(X0[3])*cos(X0[4]);R[1][0]=cos(X0[4])*sin(X0[5]);R[I1[1]=cos(X0[4])*cos(X0[5]);R[1][2]=-sin(X0[4]);R[2][0]=sin(X0[3])*cos(X0[5])+cos(X0[3])*sin(X0[4])*sin(X0[5]);R[2][l]=-sin(X0[3])*sin(X0[5])+cos(X0[3])*sin(X0[4])*cos(X0[5]);R[2][2]=cos(X0[3])*cos(X0[4]);//第一个像点的估计值，第一个点的坐标存放于X [0] [0],X [0] [1],X [0] [2]approx_x[0]=-f*(R[0][0]*(X[0][0]-X0[0])+R[1][0]*(X[0][1]-X0[1])+R[2][0]*(X[0][2]-X0[2]))/( R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2])); approx_x[1]=-f*(R[0][1]*(X[0][0]-X0[0])+R[1][1]*(X[0][1]-X0[1])+R[2][1]*(X[0][2]-X0[2]))/( R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2])); //第二个像点的估计值,第 2 个点的坐标存放于X[1][0],X[1][1],X[1][2] approx_x[2]=-f*(R[0][0]*(X[1][0]-X0[0])+R[1][0]*(X[1][1]-X0[1])+R[2][0]*(X[1][2]-X0[2]))/( R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2])); approx_x[3]=-f*(R[0][1]*(X[1][0]-X0[0])+R[1][1]*(X[1][1]-X0[1])+R[2][1]*(X[1][2]-X0[2]))/( R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2]));//第三个像点的估计值,第 3 个点的坐标存放于X[2][0],X[2][1],X[2][2] approx_x[4]=-f*(R[0][0]*(X[2][0]-X0[0])+R[1][0]*(X[2][1]-X0[1])+R[2][0]*(X[2][2]-X0[2]))/( R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2])); approx_x[5]=-f*(R[0][1]*(X[2][0]-X0[0])+R[1][1]*(X[2][1]-X0[1])+R[2][1]*(X[2][2]-X0[2]))/( R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2])); //第四个像点的估计值,第 4 个点的坐标存放于X[3][0],X[3][1],X[3][2] approx_x[6]=-f*(R[0][0]*(X[3][0]-X0[0])+R[1][0]*(X[3][1]-X0[1])+R[2][0]*(X[3][2]-X0[2]))/( R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2])); approx_x[7]=-f*(R[0][1]*(X[3][0]-X0[0])+R[1][1]*(X[3][1]-X0[1])+R[2][1]*(X[3][2]-X0[2]))/( R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2])); for(i=0;i<4;i++){//第i 个像点估计值放在approx_x[2*(i-1)],approx_x[2*i-1]/*a11*/A[2*i][0]=(R[0][0]*f+R[0][2]*approx_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a12*/A[2*i][1]=(R[1][0]*f+R[1][2]*approx_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1] -X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a13*/A[2*i][2]=(R[2][0]*f+R[2][2]*approx_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1] -X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a21*/A[2*i+1][0]=(R[0][1]*f+R[0][2]*approx_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a22*/A[2*i+1][1]=(R[1][1]*f+R[1][2]*approx_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a23*/A[2*i+1][2]=(R[2][1]*f+R[2][2]*approx_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));/*a14*/A[2*i][3]=approx_x[2*i+1]*sin(X0[4])-(approx_x[2*i]/f*(approx_x[2*i]*cos(X0[5])-approx_x[2*i+1]*sin(X0[5]))+f*cos(X0[5]))*cos(X0[4]);/*a15*/A[2*i][4]=-f*sin(X0[5])-approx_x[2*i]/f*(approx_x[2*i]*sin(X0[5])+approx_x[2*i+1]*co s(X0[5])); /*a16*/A[2*i][5]=approx_x[2*i+1];/*a24*/A[2*i+1][3]=-1*approx_x[2*i]*sin(X0[4])-(approx_x[2*i+1]/f*(approx_x[2*i]*cos(X0[5] )-approx_x[2*i+1]*sin(X0[5]))-f*sin(X0[5]))*cos(X0[4]);/*a25*/A[2*i+1][4]=-1*f*cos(X0[5])-approx_x[2*i+1]/f*(approx_x[2*i]*sin(X0[5])+approx_x[2*i+1]*cos(X0[5]));/*a26*/A[2*i+1][5]=-approx_x[2*i];}//进行常数项的初始化for(i=0;i<4;i++){L[2*i]=x[i][0]-approx_x[2*i];L[2*i+1]=x[i][1]-approx_x[2*i+1];}//A 的转置矩阵for(i=0;i<8;i++)for(j=0;j<6;j++){AT[j][i]=A[i][j];}// 实现 A 与AT 相乘int k=0;for(i=0;i<6;i++)for(j=0;j<6;j++)Asum[i][j]=0;for(i=0;i<6;i++)for(k=0;k<6;k++)for(j=0;j<8;j++)Asum[i][k]+=A T[i][j]*A[j][k];// 得到AT*A 的逆矩阵存放在inverseAsum[6][6] 中inverse(Asum);//实现矩阵Asum[6][6] 与AT[6][8] 的相乘,结果存放在result1[6][8] 中for(i=0;i<6;i++)for(j=0;j<8;j++)result1[i][j]=0;for(i=0;i<6;i++) for(k=0;k<8;k++)for(j=0;j<6;j++)result1[i][k]+=Asum[i][j]*AT[j][k];//实现resultl [6][8]与1[8]的相乘，得到结果放在result2[6]中;for(i=0;i<6;i++)result2[i]=0;for(i=0;i<6;i++)for(j=0;j<8;j++)result2[i]+=result1[i][j]*L[j];num++;for(i=0;i<6;i++){X0[i]=X0[i]+result2[i];}ofstream f7("d:\\A.txt");f7<< std::fixed;cout<<"进行第"<<num<<"次迭代带得到Xs,Ys,Zs, ^,3, K改正数分别为：\n";for(i=0;i<6;i++){ cout<<setw(12)<<result2[i];f7<<setw(12)<<result2[i];}cout<<endl<<endl;f7.close();getchar();}while(abs(result2[3]*206265.0)>6||abs(result2[4]*206265.0)>6||abs(result2[5]*20626 5.0)>6);cout<<"\n 满足限差条件结束循环，最终结果为：\n";cout<<setw(12)<<"Xs"<<setw(12)<<"Ys"<<setw(12)<<"Zs"<<setw(12)<<" p "<<setw(12)<<" 3"<<setw(12)<<" K "<<endl;ofstream f7("d:\\A.txt");f7<< std::fixed;cout.precision(4);for(i=0;i<6;i++){cout<<setw(12)<<X0[i];f7<<setw(16)<<X0[i];}f7.close();//今double XG[6][1];for(i=0;i<6;i++)XG[i][0]=result2[i];double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];multi(A,XG,AXG,8,6,1);for( i=0;i<8;i++)//计算改正数V[i][0]=AXG[i][0]-L[i];transpose (V,VT,1,8);multi(VT,V ,VTV,1,8,1);m0=VTV[0][0]/2;cout<<endl;ofstream f6("d:\\what.txt");// f6<< std::fixed;for(i=0;i<6;i++){for(int j=0;j<6;j++){D[i][j]=m0*Asum[i][j]; cout<<setw(10)<<D[i][j];f6<<setw(15)<<D[i][j];}cout<<endl;f6<<endl;}for(i=0;i<6;i++)cout<<sqrt(D[i][i])<<endl;f6.close();//屏幕输出误差方程系数阵、常数项、改正数//getchar();return 0;}void inverse(double c[n][n]){ int i,j,h,k; double p;double q[n][12];for(i=0;i<n;i++)// 构造高斯矩阵for(j=0;j<n;j++)q[i][j]=c[i][j];for(i=0;i<n;i++)for(j=n;j<12;j++){if(i+6==j)q[i][j]=1;elseq[i][j]=0;}for(h=k=0;k<n-1;k++,h++)// 消去对角线以下的数据for(i=k+1;i<n;i++){ if(q[i][h]==0) continue; p=q[k][h]/q[i][h];for(j=0;j<12;j++){q[i][j]*=p;q[i][j]-=q[k][j];}}for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据for(i=k-1;i>=0;i--) {if(q[i][h]==0)continue;p=q[k][h]/q[i][h];for(j=0;j<12;j++){q[i][j]*=p; q[i][j]-=q[k][j];}}for(i=0;i<n;i++)// 将对角线上数据化为 1{ p=1.0/q[i][i];for(j=0;j<12;j++)q[i][j]*=p;}for(i=0;i<n;i++) // 提取逆矩阵for(j=0;j<n;j++)c[i][j]=q[i][j+6];}template<typename T1,typename T2>void transpose(T1*mat1,T2*mat2,int a,int b){ int i,j;for(i=0;i<b;i++)for(j=0;j<a;j++)mat2[j][i]=mat1[i][j];return;}template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c){for(i=0;i<a;i++) for(j=0;j<c;j++) result[i][j]=0; for(k=0;k<b;k++) result[i][j]+=mat1[i][k]*mat2[k][j];return;程序截图五、结果分析六、心得体会本次实验是关于单张像片的空间后方交会的计算原理及实现过程的课程作业， 其中解题步骤和公式是基本， 熟练运用公式并且成功编程是重点。