摄影测量学空间后方交会实验报告(共13页)

1、摄影测量学实验报告实验一、单像空间后方交会学院： 建测学院 班级： 测绘082 姓名： 肖澎 学号： 15 一 实验目的1.深入了解单像空间后方交会的计算过程；2.加强空间后方交会基本公式和误差方程式，法线方程式的记忆；3.通过上机调试程序加强动手能力的培养。二实验原理以单幅影像为基础，从该影像所覆盖地面范围内若干控制点和相应点的像坐标量测值出发，根据共线条件方程，求解该影像在航空摄影时刻的相片外方位元素。三 实验内容1.程序图框图2.实验数据（1）已知航摄仪内方位元素f153.24mm，XoYo0。限差0.1秒（2）已知4对点的影像坐标和地面坐标：影像坐标地面坐标x（mm）y（mm）X(m)

2、Y（m）Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.313.实验程序using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace ConsoleApplication3 class Program static void Main() /输

3、入比例尺，主距，参与平参点的个数 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"); str

4、ing n1 = Console.ReadLine(); int n = (int)Convert.ToSingle(n1); /像点坐标的输入代码 double arr1 = new double2 * n; /1.像点x坐标的输入 for (int i = 0; i < n; i+) Console.WriteLine("请输入已进行系统误差改正的像点坐标的x0值：r", i+1); string u = Console.ReadLine(); for (int j = 0; j < n; j += 2) arr1j = (double)Convert.To

5、Single(u); /2.像点y坐标的输入 for (int i = 0; i < n; i+) Console.WriteLine("请输入已进行系统误差改正的像点坐标的y0值：r", i+1); string v = Console.ReadLine(); for (int j = 1; j < n; j += 2) arr1j = (double)Convert.ToSingle(v); /控制点的坐标输入代码 double, arr2 = new doublen, 3; /1.控制点X坐标的输入 for (int j = 0; j < n; j+

6、) Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的X0值：r", j+1); string u = Console.ReadLine(); arr2j , 0 = (double)Convert.ToSingle(u); /2.控制点Y坐标的输入 for (int k = 0; k < n; k+) Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的Y0值：r", k+1); string v = Console.ReadLine(); arr2k , 1 = (double)Conver

7、t.ToSingle(v); /3.控制点Z坐标的输入 for (int p =0; p < n; p+) Console.WriteLine("请输入控制点在地面摄影测量坐标系的坐标的Z0值：r", p+1); string w = Console.ReadLine(); arr2p , 2 = (double)Convert.ToSingle(w); /确定外方位元素的初始值 /1.确定Xs的初始值： double Xs0 = 0; double sumx = 0; for (int j = 0; j < n; j+) double h = arr2j, 0

8、; sumx += h; Xs0 = sumx / n; /2.确定Ys的初始值： double Ys0 = 0; double sumy = 0; for (int j = 0; j < n; j+) double h = arr2j, 1; sumy += h; Ys0 = sumy / n; /3.确定Zs的初始值： double Zs0 = 0; double sumz = 0; for (int j = 0; j <= n - 1; j+) double h = arr2j, 2; sumz += h; Zs0 = sumz / n; double 0 = 0; doub

9、le 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 double3, 3; for (int i = 0; i < 3; i+) arr3i, i = 1; double a1 = arr30, 0; double a2 = arr30, 1; double a3 = ar

10、r30, 2; double b1 = arr31, 0; double b2 = arr31, 1; double b3 = arr31, 2; double c1 = arr32, 0; double c2 = arr32, 1; double c3 = arr32, 2; /*利用线元素的初始值和控制点的地面坐标，代入共线方程（3-5-2）， * 逐点计算像点坐标的近似值*/ /1.定义存放像点近似值的数组 double arr4 = new double2 * n;/-近似值矩阵 /2.逐点像点坐标计算近似值 /a.计算像点的x坐标近似值（x） for (int i = 0; i &l

11、t; 2 * n; i += 2) for (int j = 0; j < n; j+) arr4i = -f * (a1 * (arr2j, 0 - Xs0) + b1 * (arr2j, 1 - Ys0) + c1 * (arr2j, 2 - Zs0) / (a3 * (arr2j, 0 - Xs0) + b3 * (arr2j, 1 - Ys0) + c3 * (arr2j, 2 - Zs0); /b.计算像点的y坐标近似值（y） for (int i = 1; i < 2 * n; i += 2) for (int j = 0; j < n; j+) arr4i =

12、-f * (a2 * (arr2j, 0 - Xs0) + b2 * (arr2j, 1 - Ys0) + c2 * (arr2j, 2 - Zs0) / (a3 * (arr2j, 0 - Xs0) + b3 * (arr2j, 1 - Ys0) + c3 * (arr2j, 2 - Zs0); /逐点计算误差方程式的系数和常数项，组成误差方程： double, arr5 = new double2 * n, 6; /-系数矩阵（A） /1.计算dXs的系数 for (int i = 0; i < 2 * n; i += 2) arr5i, 0 = -1 / m; /-f/H = -1

13、/m /2.计算dYs的系数 for (int i = 1; i < 2 * n; i += 2) arr5i, 1 = -1 / m; /-f/H = -1/m /3.a.计算误差方程式Vx中dZs的系数 for (int i = 0; i < 2 * n; i += 2) arr5i, 2 = -arr1i / m * f; /3.b.计算误差方程式Vy中dZs的系数 for (int i = 1; i < 2 * n; i += 2) arr5i, 2 = -arr1i / m * f; /4.a.计算误差方程式Vx中d的系数 for (int i = 0; i <

14、; 2 * n; i += 2) arr5i, 3 = -f * (1 + arr1i * arr1i / f * f); /4.a.计算误差方程式Vy中d的系数 for (int i = 1; i < 2 * n; i += 2) arr5i, 3 = -arr1i - 1 * arr1i / f; /5.a.计算误差方程式Vx中d的系数 for (int i = 0; i < 2 * n; i += 2) arr5i, 4 = -arr1i * arr1i + 1 / f; /5.b.计算误差方程式Vy中d的系数 for (int i = 1; i < 2 * n; i

15、+= 2) arr5i, 4 = -f * (1 + arr1i * arr1i / f * f); /6.a.计算误差方程式Vx中dk的系数 for (int i = 0; i < 2 * n; i += 2) arr5i, 5 = arr1i + 1; /6.b.计算误差方程式Vy中dk的系数 for (int i = 1; i < 2 * n; i += 2) arr5i, 5 = -arr1i - 1; /定义外方位元素组成的数组 double arr6 = new double6;/-外方位元素改正数矩阵（X） /定义常数项元素组成的数组 double arr7 = ne

16、w double2 * n;/-常数矩阵（L） /计算lx的值 for (int i = 0; i < 2 * n; i += 2) arr7i = arr1i - arr4i; /将近似值矩阵的元素代入 /计算ly的值 for (int i = 1; i <= 2 * (n - 1); i += 2) arr7i = arr1i - arr4i; /将近似值矩阵的元素代入 /* 对于所有像点的坐标观测值，一般认为是等精度量测，所以权阵P为单位阵. 所以X=(ATA)-1ATL */ /1.计算AT double, arr5T = new double6, 2 * n; for (

17、int i = 0; i < 6; i+) for (int j = 0; j < 2 * n; j+) arr5Ti, j = arr5j, i; /A的转置与A的乘积,存放在arr5AA中 double, arr5AA = new double6, 6; for (int i = 0; i < 6; i+) for (int j = 0; j < 6; j+) arr5AAi, j = 0; for (int l = 0; l < 2 * n; l+) arr5AAi, j += arr5Ti, l * arr5l, j; nijuzhen(arr5AA);

18、/arr5AA经过求逆后变成原矩阵的逆矩阵 /arr5AA * arr5T存在arr5AARAT double, arr5AARAT = new double6, 2 * n; for (int i = 0; i < 6; i+) for (int j = 0; j < 2 * n; j+) arr5AARATi, j = 0; for (int p = 0; p < 6; p+) arr5AARATi, j += arr5AAi, p * arr5Tp, j; /计算arr5AARAT x L，存在arrX中 double arrX = new double6; for (

19、int i = 0; i < 6; i+) for (int j = 0; j < 1; j+) arrXi = 0; for (int vv = 0; vv < 6; vv+) arrXi += arr5AARATi, vv * arr7vv; /计算外方位元素值 double Xs, Ys, Zs, , , K; Xs = Xs0 + arrX0; Ys = Ys0 + arrX1; Zs = Zs0 + arrX2; = 0 + arrX3; = 0 + arrX4; K = K0 + arrX5; for (int i = 0; i <= 2; i+) Xs +

20、= arrX0; Ys += arrX1; Zs += arrX2; += arrX3; += arrX4; K += arrX5; 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 double6, 6; int i, j, k; int row = 0; int col = 0; double max, te

21、mp; int p = new int6; for (i = 0; i < 6; i+) pi = i; Bi, 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(ai, j); if (max < temp) max = temp; row = i; col = j; /交换行列，将主元调整到k行k列上 if (row != k) for (j = 0; j < 6; j+) temp = arow, j; arow, j = ak, j; ak, j = temp; temp = Brow, j; Brow, j = Bk, j; Bk, j = temp; i = prow; prow = pk; pk = i; if (col != k) for (i = 0; i < 6; i+) temp =