摄影测量实习报告

1、实验报告一：单像空间后方交会程序设计报告内容：实验报告主要是针对实验目的和实验过程进行描述，并对所得实验 结果进行分析和评价，指出实验注意事项和对本实验提出建议。(一)已知数据与作业要求已知数据：XyXYZ1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31内方位元素？?= 153.24，?= ?= 0 ;摄影比例尺分母 m=4000Q要求：已知4对点的影像坐标和地面坐标

2、，计算近似垂直下空间后方交会的 解，得到6个外方位元素的值。(二) 单像空间后方交会计算原理空间后方交会的数学模型是共线条件方程，其具体形式可写为：?(? ? + ?(? ? + ?(? ?(?- ?+ ?(? ? + ?(? ?)?(?- ? + ?(? ? + ?込？2 ?= -? ?(?- ? + ?(?- ? + ?- ?式中(?为像点的像平面坐标，(?？?为像点对应地面点的物方空间坐标，？为相 片主距，？?为外方位线元素，？?= 1,2,3)为相片外方位角元素所组成的 9个方向余弦。对任一控制点，其地面坐标(?？?秽和对应像点坐标(?都是已知的，代入 共线条件方程可以列出两个方程式，

3、因此从理论上讲只要有三个控制点就可以列出 6个方程，从而求解出6个外方位元素。但为了避免粗差，提高测量精度，应有多 余观测，因此，在实际应用中一般需要4个或更多的控制点。(三) 单像空间后方交会的计算过程1. 获取已知数据：从航摄资料中获取像片比例尺 1/m、内方位元素？?,？?，获 取控制点的空间坐标X，丫，Zo2. 量测控制点的像点坐标：量测控制点框标坐标，并经像主点改正得到像点坐 标x，y。3. 确定未知数的初始值：在竖直摄影情况下，角元素的初始值给0,线元素中?0= ? ?Q ?可取四个角上控制点坐标的平均值?=韦？ ?= ?= ? = 7?1 = ? = ORO式中：仃为摄影比例尺分

4、母；n为控制点个数4. 用三个角元素的初始值计算方向余弦，并组成旋转矩阵?=:cos?cos? Sin?sin?Sin?=-COS ?si n ? sin? ?Si n ? ?cos ?=-Sin ?cos? ? ?=cos?Sin ?= ? ? ? ?=cos?cos ? ?2? ?=-Sin ?1?=Sin ?cos? cos?sin ?sin ?=-Sin ?Sin ?+ cos?sin ?cos?=cos?cos?5.用所取未知数的初始值和控制点的地面坐标，代入共线方程，逐点计算像点(x)、 (y)。坐标的近似值,?=?= ?(?-?)+?(?-?+?iX?-?_?(?-?)+?(?-

5、?+?3(?-?)“？?(？?-?)+?(?-?+?2(?-?)(?)=(?)=(?2 =(?2 =(?3 =(?3 =(?4 =(?4 =“？?(？- ? + ?(?- ? + ?(?- ? ?(?- ? + ?(?- ? + ?(?- ? “？(？- ? + ?(?-鋼 + ?(?- ? ?(?- ? + ?(?-漏 + ?(?- ? ?(?- ? + ?(?- ? + ?(?- ? -?(? - ? + ?(?-鬲 + ?孔？- E ?(?- ? + ?(?-鋼 + ?(?- ? _?(? - ? + ?(?-鋼 + ?乱？- ? _?(?- ? + ?(?-錮 + ?(?- ? _?(

6、? - ? + ?(?-钢 + ?乱？- ?) “？(？- ? + ?(?-鋼 + ?(?- ? ?(? - ? + ?(?-钢 + ?(?- ? ?(?- ? + ?(?-窮 + ?(?- ? -?(? - ? + ?(?- ?) + ?- E _?(?- ? + ?(?-笏 + ?(?- ? _?(? - ? + ?(?-鋼 + ?乱？- ?6. 用像点的观测值和由5中计算的近似值，计算每个点的常数项 ？ ?+?+ ?(?初？刃+?(?衿？?)+?寅？务？?)? ?(?-?)+?(?-?+?3<?7?(?2123 4)?(?7?)+?(?*?+?(?7?)(?，3,?(?S-?)+?

7、(?-?+?(?)7. 计算法方程的系数矩阵？?与常数项？?(?-?)+?(?-?)+?3(?-?)? Z= -?2?= 01?=-?4 =:-?(1+祠?5 =?1?=?1 =0?=- ?=-?4 =Z -?5 =-?(1 +Tr)?6 =-?I6X 08 ?2?3?4?5?6?2?3?4?5?6?2?3?4?5?6?2?3?4?5?46?2?3?5?6?2?3?4?5?6?2?3?4?5?6?2?3?4?5?6 ?1?1?1? ? ?礒 000?注意：A矩阵中的（?不是已知的像点坐标，而是将每次迭代后所得到的外方 位元素代入共线条件方程计算出来的，A矩阵在迭代过程中也是变化的8. 求解法方

8、程，得到未知数的改正数。?= （?）1 ?9. 检查计算结果是否收敛：将改正数与限差比较，小于限差则计算终止，否则 用新的近似值重复4-9的计算，直到满足要求为止。?= ?M+?需?=T?"1+ ?=?-1+?勇=?畀-1+?液? =?參-1+ ?层? =?-1+?式中：？?弋表迭代次数。10. 空间后方教会的精度估算。 ?= ?1?+?2 ?3?+?4 ?+?5 ?5 K ? - ? = ?1?22 ?+?3?+?4 ?+?25 ?26 ? - ?由单位权中误差?= ?和6个参数的协因数阵？?= (?)1得出参数？?的中误差为？?= ?丿?爲？（四）算法流程图（五）计算结果MATL

9、A计算结果输出截图:JCs =3.979a44340D520e÷C 4fai =0. W39S5559117g04¥3 =N 7476464S4D27101e4C4Onega OJ 002Ll3670028136ZS =7. 57368833LD23273t÷C3kaa -OJ 065r64B10370,5IJlII -R =1.125401304950231L 2436735942644O- 6377 &56790616O. 9770900f98137C. 0675j406U65390. 003985629934395-0. C6526C4663190.

10、 957715272B03382-0. 0021137265110190. (MQlaISO7341 L-QI 004119272433176Ol 001839750147646OJ 999989823405350Ol 0001596811695 诟OL DOOo"243527C#计算结果输出截图（六）实验心得与体会（七）程序实现代码见附录附录：单像空间后方交会代码 matlab 程序 : clc;clear;format long%单位统一为米xy=1/1000*-86.15 -68.99;-53.40 82.21;-14.78 -76.63;10.46 64.43;% 导入像点坐

11、标XYZ=36589.41 25273.32 2195.17;37631.08 31324.51 728.69;.39100.97 24934.98 2386.50;40426.54 30319.81 757.31; % 导入物点坐标 f=0.15324;x0=0;y0=0;m=40000; %摄影比例尺分母fai=0;omega=0;kappa=0; % 初始外方位角元素Xs=0;Ys=0;for i0=1:4;Xs=Xs+XYZ(i0,1);Ys=Ys+XYZ(i0,2);endZs=m*f;Xs=Xs/4;Ys=Ys/4; %初始外方位线元素R=zeros(3);L=zeros(8,1)

12、;delta=zeros(6,1);A=zeros(8,6);x=zeros(4,1);y=zeros(4,1);t=0;while 1%旋转矩阵 R(1,1)=cos(fai)*cos(kappa)-sin(fai)*sin(omega)*sin(kappa);R(1,2)=-cos(fai)*sin(kappa)-sin(fai)*sin(omega)*cos(kappa); R(1,3)=-sin(fai)*cos(omega);R(2,1)=cos(omega)*sin(kappa);R(2,2)=cos(omega)*cos(kappa);R(2,3)=-sin(omega); R(

13、3,1)=sin(fai)*cos(kappa)+cos(fai)*sin(omega)*sin(kappa); R(3,2)=-sin(fai)*sin(kappa)+cos(fai)*sin(omega)*cos(kappa); R(3,3)=cos(fai)*cos(omega);for i1=1:4 %共线条件方程 x(i1,1)=-f*(R(1,1)*(XYZ(i1,1)-Xs)+R(2,1)*(XYZ(i1,2)-Ys)+R(3,1)*(XYZ(i1,3)-Zs)/. (R(1,3)*(XYZ(i1,1)-Xs)+R(2,3)*(XYZ(i1,2)-Ys)+R(3,3)*(XYZ(

14、i1,3)-Zs);y(i1,1)=-f*(R(1,2)*(XYZ(i1,1)-Xs)+R(2,2)*(XYZ(i1,2)-Ys)+R(3,2)*(XYZ(i1,3)- Zs)/.(R(1,3)*(XYZ(i1,1)-Xs)+R(2,3)*(XYZ(i1,2)-Ys)+R(3,3)*(XYZ(i1,3)-Zs); %误差方程常数项L(2*i1-1,1)=xy(i1,1)-x(i1,1); L(2*i1,1)=xy(i1,2)-y(i1,1);%系数矩阵 A(2*i1-1,1)=-f/(Zs-XYZ(i1,3);A(2*i1-1,2)=0; A(2*i1-1,3)=-x(i1,1)/(Zs-XY

15、Z(i1,3);A(2*i1-1,4)=-f*(1+x(i1,1)*x(i1,1)f2);A(2*i1-1,5)=-x(i1,1)*y(i1,1)/f;A(2*i1-1,6)=y(i1,1);A(2*i1,1)=0; A(2*i1,2)=-f(Zs-XYZ(i1,3);A(2*i1,3)=-y(i1,1)(Zs-XYZ(i1,3);A(2*i1,4)=-x(i1,1)*y(i1,1)f; A(2*i1,5)=-f*(1+y(i1,1)*y(i1,1)f2);A(2*i1,6)=-x(i1,1);enddelta=(A'*A)A'*L; % 改正数Xs=Xs+delta(1,1)

16、;Ys=Ys+delta(2,1);Zs=Zs+delta(3,1);fai=fai+delta(4,1);omega=omega+delta(5,1); kappa=kappa+delta(6,1);t=t+1;if (abs(delta(4,1)<10(-6)&&(abs(delta(5,1)<10(-6)&&( abs(delta(6,1)<10(- 6)breakendendV=A*delta-L;VV=0;for i2=1:8VV=VV+V(i2)*V(i2);endm0=sqrt(VV2); %单位权中误差mi=zeros(6,1)

17、; %未知数中误差 ffcinv=inv(A'*A); %法方程的逆for i3=1:6 mi(i3)=m0*sqrt(ffcinv(i3,i3);end %输出外方位元素，旋转矩阵和未知数中误差XsYsZs fai omega kappa R miC#程序：using System;using System.Collections.Generic;using System.Linq;using System.Text;using System.Threading.Tasks;namespace resectionclass Program/转置static double, Trans

18、pose(double, a)double, b = new doublea.GetLength(1), a.GetLength(0);for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)bj, i = ai, j;return b;/相加static double, addition(double, a, double, b)double, c = new doublea.GetLength(0), a.GetLength(1);if (a.GetLength(0) = b.G

19、etLength(0) && (a.GetLength(1) = b.GetLength(1)for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)ci, j = ai, j + bi, j;return c;elseSystem.Exception e = new Exception(" 两矩阵维数不等，检查数据");throw e;/相减static double, subtraction(double, a, double, b)doubl

20、e, c = new doublea.GetLength(0), a.GetLength(1);if (a.GetLength(0) = b.GetLength(0) && (a.GetLength(1) = b.GetLength(1)for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)ci, j = ai, j - bi, j;return c;elseSystem.Exception e = new Exception(" 两矩阵维数不等，检查数

21、据");throw e;/相乘static double, multiply(double, a, double, b)double, c = new doublea.GetLength(0), b.GetLength(1);if (a.GetLength(1) != b.GetLength(0)System.Exception e = new Exception(" 左矩阵列数与右矩阵行数 不等，检查数据 ");throw e;elsedouble temp;for (int i = 0; i < a.GetLength(0); i+)for (int j

22、 = 0; j < b.GetLength(1); j+)temp = 0;for (int k = 0; k < a.GetLength(1); k+)temp += ai, k * bk, j;ci, j = temp;return c;/行交换static double, swap(int temp, int temp1, double, a)double b=new doublea.GetLength(1);double c=new doublea.GetLength(1);for (int i = 0; i < a.GetLength(1); i+)bi = ate

23、mp, i;ci = atemp1, i;for (int i = 0; i < a.GetLength(1); i+)atemp1, i=bi; atemp, i= ci;return a;/行列式static double Det(double, a)double c = 1;if (a.GetLength(0) != a.GetLength(1)据");System.Exception e = new Exception(" 矩阵行列数不等，检查数 throw e;for (int k = 0; k < a.GetLength(0); k+)for (in

24、t i = k + 1; i < a.GetLength(0); i+) double temp2; if (ai, k != 0) temp2 = ai, k / ak, k; for (int j = 0; j < a.GetLength(1); j+) ai, j -= ak, j * temp2;/上三角矩阵for (int i = 0; i < a.GetLength(0); i+)c *= ai, i;return c;/ 求逆static double, inverse(double, a)double, b = new doublea.GetLength(0)

25、, 2 * a.GetLength(0); for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(0); j+)bi, j = ai, j;for (int j = a.GetLength(0); j < 2 * a.GetLength(0); j+)if (j - a.GetLength(0) = i)bi, j = 1;elsebi, j = 0;/以上步骤为构造 (A,E) 矩阵double temp = Math.Abs(b0, 0);int temp1 = 0;for (int

26、i = 0; i < a.GetLength(0); i+)if (temp <= Math.Abs(bi, 0)temp = bi, 0;temp1 = i;/以上选出第一列中最大值b = swap(0, temp1, b);for (int k = 0; k < b.GetLength(0); k+)for (int i = k + 1; i < b.GetLength(0); i+)double temp2;if (bi, k != 0)temp2 = bi, k / bk, k;for (int j = 0; j < b.GetLength(1); j+)

27、 bi, j -= bk, j * temp2; /上三角矩阵 double c = Det(a); if (c = 0) System.Exception e = new Exception(" 矩阵行列式为零，检查数据");throw e;/检查矩阵行列式是否为 0for (int k = b.GetLength(0) - 1; k >= 0; k-)for (int i = k - 1; i >= 0; i-)double temp3;if (bi, k != 0)temp3 = bi, k / bk, k;for (int j = 0; j < b

28、.GetLength(1); j+) bi, j -= bk, j * temp3;/已处理成对角线矩阵for (int i = 0; i < b.GetLength(0); i+)double temp4 = bi, i;for (int j = 0; j < b.GetLength(1); j+)bi, j /= temp4;/已处理成 (E,A 逆 )的形式for (int i = 0; i < b.GetLength(0); i+)for (int j = b.GetLength(0); j < b.GetLength(1); j+)ai, j - b.GetL

29、ength(0) = bi, j;/此时 a 为原始输入矩阵的逆 return a;static void Main(string args)double, xy = new double, -86.15, -68.99 , -53.40, 82.21 , - 14.78, -76.63 , 10.46, 64.43 ;double, XYZ = new double, 36589.41, 25273.32, 2195.17 , 37631.08, 31324.51, 728.69 , 39100.97, 24934.98, 2386.50 , 40426.54, 30319.81, 757.

30、31 ;double f = 0.15324, m = 40000;double fai = 0, omega = 0, kappa = 0;/ 角元素 double Xs = 0, Ys = 0, Zs = 0;/ 线元素 double, R = new double3, 3;/ 旋转矩阵 double, L = new double2*xy.GetLength(0), 1;/ 常数项 double, A = new double2 * xy.GetLength(0), 6;/ 系数矩阵 double, x = new doublexy.GetLength(0), 1; double, y

31、= new doublexy.GetLength(0), 1; double, Diff = new double6,1;/ 改正数矩阵 double, V = new double2 * xy.GetLength(0),1;/ 误差矩阵 double, FFCinv = new double6, 6;/ 法方程的逆矩阵 double, mi = new double6, 1;/ 未知数中误差 double vv = 0,m0=0;/vv 为残差平方和 m0 为单位权中误差 int t = 0;/ 迭代次数 for (int i = 0; i < 4; i+) for (int j =

32、0; j < 2; j+)xyi, j /= 1000;for (int i = 0; i < 4; i+)Xs = Xs + XYZi, 0;Ys = Ys + XYZi, 1;Zs = m * f;Xs = Xs / 4;Ys = Ys / 4; while (true) /旋转矩阵R0, 0 = Math.Cos(fai) * Math.Cos(kappa) - Math.Sin(fai) * Math.Sin(omega) * Math.Sin(kappa);R0, 1 = -Math.Cos(fai) * Math.Sin(kappa) - Math.Sin(fai) M

33、ath.Sin(omega) * Math.Cos(kappa);R0, 2 = -Math.Sin(fai) * Math.Cos(omega);R1, 0 = Math.Cos(omega) * Math.Sin(kappa);R1, 1 = Math.Cos(omega) * Math.Cos(kappa);R1, 2 = -Math.Sin(omega);R2, 0 = Math.Sin(fai) * Math.Cos(kappa) + Math.Cos(fai) Math.Sin(omega) * Math.Sin(kappa);R2, 1 = -Math.Sin(fai) * Ma

34、th.Sin(kappa) + Math.Cos(fai) Math.Sin(omega) * Math.Cos(kappa);R2, 2 = Math.Cos(fai) * Math.Cos(omega);for (int i = 0; i < 4; i+) /共线条件方程 xi, 0 = -f * (R0, 0 * (XYZi, 0 - Xs) + R1, 0 * (XYZi, 1 - Ys) + R2, 0 * (XYZi, 2 - Zs) /(R0, 2 * (XYZi, 0 - Xs) + R1, 2 * (XYZi, 1 - Ys) + R2, 2 * (XYZi, 2 -

35、Zs);yi, 0 = -f * (R0, 1 * (XYZi, 0 - Xs) + R1, 1 * (XYZi, 1 - Ys) + R2, 1 * (XYZi, 2 - Zs) /(R0, 2 * (XYZi, 0 - Xs) + R1, 2 * (XYZi, 1 - Ys) + R2, 2 * (XYZi, 2 - Zs);/误差方程常数项 L2 * i, 0 = xyi, 0 - xi, 0;L2 * i+1, 0 = xyi, 1 - yi, 0;/系数矩阵A2 * i, 0 = -f / (Zs - XYZi, 2);A2 * i, 1 = 0;A2 * i, 2 = -xi, 0 / (Zs - XYZi, 2);A2 * i, 3 = -f * (1 + xi, 0 * xi, 0 / (f * f);A2 * i, 4 = -xi, 0 * yi, 0 / f;A2 * i, 5 = yi, 0;A2 * i+1, 0 = 0;A2 * i+1, 1 = -f / (Zs - XYZi