下载文档原格式
单像空间后方交会名词解释
单像空间后方交会是摄影测量学中的一个重要概念，它是指利用单个影像进行地物测量和定位的方法。
在单像空间后方交会中，通过对单张影像进行分析，可以确定地面上物体的位置和形状。
这个过程涉及到对影像中的特征点进行识别和匹配，然后利用相机内外参数以及影像上的像点坐标来计算地物的三维坐标。
单像空间后方交会的过程包括以下几个步骤，首先是对影像进行预处理，包括去畸变、影像配准等操作；然后是特征点的提取和匹配，这一步是通过计算机视觉算法来实现的，可以利用角点、边缘等特征来进行匹配；接下来是相机内外参数的标定，这一步是为了将像素坐标转换为实际世界坐标而进行的；最后是利用已知的相机参数和像点坐标来计算地物的三维坐标。
单像空间后方交会在航空摄影、遥感影像解译和地图制图等领域有着广泛的应用。
它可以通过对单张影像的处理，实现对地物的测量和定位，为地理信息系统和地图制图提供了重要的数据基础。
同时，随着计算机视觉和图像处理技术的不断发展，单像空间后方交会的精度和效率也在不断提高，为各种应用领域提供了更加可靠和精确的地物信息。
单像空间后方交会实习报告一、实习目的单像空间后方交会是摄影测量中确定像片外方位元素的重要方法。
通过本次实习，旨在深入理解单像空间后方交会的基本原理和计算过程，熟练掌握相关软件的操作，提高对摄影测量数据处理的实践能力，并培养解决实际问题的思维和方法。
二、实习原理单像空间后方交会的目的是利用像片上的像点坐标以及相应的地面控制点坐标，通过数学模型求解像片的外方位元素（三个线元素 Xs、Ys、Zs 和三个角元素φ、ω、κ）。
其基本原理基于共线条件方程，即摄影中心、像点和相应的地面点位于同一条直线上。
共线条件方程可以表示为：＼＼begin{align}x x_0&＝ f\frac{a_1(X X_s) ＋ b_1(Y Y_s) ＋ c_1(Z Z_s)｝｛a_3(X X_s) ＋ b_3(Y Y_s) ＋ c_3(Z Z_s)｝＼＼y y_0&＝ f\frac{a_2(X X_s) ＋ b_2(Y Y_s) ＋ c_2(Z Z_s)｝｛a_3(X X_s) ＋ b_3(Y Y_s) ＋ c_3(Z Z_s)｝＼end{align}＼其中，＼（（x,y)＼）为像点坐标，＼（（x_0,y_0)＼）为主点坐标，＼（f\）为摄影机焦距，＼（（X,Y,Z)＼）为地面点的物方空间坐标，＼（（X_s,Y_s,Z_s)＼）为摄影中心的物方空间坐标，＼（（a_1,b_1,c_1)，（a_2,b_2,c_2)，（a_3,b_3,c_3)＼）为由角元素φ、ω、κ 构成的旋转矩阵的元素。
三、实习数据本次实习使用了一组航空像片，像片比例尺为 1:5000，焦距为152mm，像主点坐标为\（（x_0,y_0)＝（5000mm,5000mm)＼）。
同时，提供了 6 个均匀分布在像片范围内的地面控制点的物方空间坐标和像点坐标。
四、实习步骤1、数据准备整理地面控制点的物方空间坐标和像点坐标，确保数据的准确性。
输入像片的基本参数，如像主点坐标、焦距等。
单像空间后方交会测绘学院 成晓倩1 概述1.1 定义利用一定数量的地面控制点和对应像点坐标求解单张像片外方位元素的方法称为空间后方交会。
1.2 所需控制点个数与分布共线条件方程的一般形式为：⎪⎪⎩⎪⎪⎨⎧-+-+--+-+--=--+-+--+-+--=-)()()()()()()()()()()()(33322203331110S S S S S S S S S S S S Z Z c Y Y b X X a Z Z c Y Y b X X a f y y Z Z c Y Y b X X a Z Z c Y Y b X X a f x x （1）式中包含有六个外方位元素，即κωϕ、、、、、S S S Z Y X ，只有确定了这六个外方位元素的值，才能利用共线条件方程真正确定一张像片的任一像点与对应地面点的坐标关系。
个数：对任一控制点，我们已知其地面坐标)(i i i Z Y X 、、和对应像点坐标)(i i y x 、，代入共线条件方程可以列出两个方程式，因此，只少需要3个控制点才能解算出六个外方位元素。
在实际应用中，为了避免粗差，应有多余检查点，因此，一般需要4～6个控制点。
分布：为了最有效地控制整张像片，控制点应均匀分布于像片边缘，如下图所示。
由于共线条件方程是非线性的，直接答解十分困难，所以首先将共线方程改化为线性形式，然后再答解最为简单的线性方程组。
2 空间后方交会的基本思路分布合理 分布合理 分布不合理2.1 共线条件方程线性化的基本思路在共线条件方程中，令)()()()()()()()()(333222111S S S S S S S S S Z Z c Y Y b X X a Z Z Z c Y Y b X X a Y Z Z c Y Y b X X a X -+-+-=-+-+-=-+-+-= （2） 则共线方程变为⎪⎪⎩⎪⎪⎨⎧-=--=-ZY fy y Z Xf x x 00 （3） 对上式两侧同乘Z ，并移至方程同侧，则有⎩⎨⎧=-+=-+0)(0)(00Z y y Y f Z x x X f （4） 令⎩⎨⎧-+=-+=Zy y Y f Fy Zx x X f Fx )()(00 （5） 由于上式是共线方程的变形，因此，Fy Fx 、是κωϕ、、、、、S S S Z Y X 的函数。
单像空间后方交会原理非常简单，稍稍懂一点线性代数和微分学知识就很容易理解。
但它的意义很大，所有的基于视觉技术的定位（测量）都是从它扩展来的，没有它就没有后面的复杂算法的演进，以面向更复杂的场景以及需求。
这块内容一般来说就是在已知地面上若干点的地面坐标以后，反求该相应摄影光束的外参数，当用以作摄影机的标定时，还可借以同时求出摄影的内参数。
具体的说就是一个多元线性模型，根据不同的需求可以基于它进一步构建各种平差模型。
摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘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．实验结果四．实验总结此次实验让我深入了解单像空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。
单像空间后方交会原理你知道单像空间后方交会吗？这可是摄影测量里一个超有趣的概念呢！咱们先来说说啥是单像空间后方交会。
想象一下，你拿着相机拍了一张照片，这张照片里有好多好多的景物。
那单像空间后方交会呢，就是通过这一张照片里的信息，去算出拍摄这张照片的时候，相机在空间里的位置和姿态。
比如说，照片里有一座山，还有一条河，还有几棵大树。
那咱们怎么通过这些东西来知道相机当时在哪，朝哪个方向呢？这就用到单像空间后方交会啦！这当中有几个关键的东西哦。
一个是控制点，就好像是我们的“小帮手”。
这些控制点是我们事先知道它们在空间里准确位置的点。
比如说，有个特别明显的大石头，我们知道它在地球上的坐标是多少。
然后呢，还有像片的内方位元素。
这就像是相机的“小秘密”，比如说相机的焦距啦等等。
那怎么通过这些来算出相机的位置和姿态呢？这就像是一个解谜的过程！咱们得先把照片上控制点的像点坐标找出来，这就像是在照片里给这些控制点“定位”。
然后呢，根据一些数学公式和算法，把这些坐标啊、内方位元素啊、控制点的空间坐标啊等等都放到一起，就像是把一堆拼图的碎片拼起来。
这个过程可不容易哦，得算好多好多的数学式子。
但是别担心，咱们聪明的科学家们早就想出了办法，有各种软件和工具能帮咱们完成这些复杂的计算。
你可能会想，这有啥用啊？用处可大啦！比如说，我们要做地图，要对一个地方进行测量，单像空间后方交会就能帮我们得到相机的位置和姿态，这样就能更准确地知道照片里的东西在实际空间里的位置啦。
而且哦，现在科技越来越发达，单像空间后方交会的精度也越来越高。
这就像是我们的眼睛越来越厉害，能看得更清楚，更准确！想象一下，如果没有单像空间后方交会，那我们看到的照片就只是一张好看的图片，没办法知道那么多背后的信息。
但是有了它，一张照片就像是一个装满了秘密的宝盒，我们可以一点点地解开，发现更多有趣的东西。
怎么样，是不是觉得单像空间后方交会很神奇很有趣呀？希望我讲得能让你明白这个有点复杂但又超级酷的原理！。
简述单像空间后方交会的程序设计步骤
单像空间后方交会是一种用于测量摄影点在三维空间中位置的方法。
以下是简述的程序设计步骤：
1.读取摄影测量数据：首先，从摄影测量设备（如相机）中读取图像和相关的内参数据，包括相机的焦距、像点大小等。
2.图像处理：对读取的图像进行预处理。
可能需要进行去畸变操作，校正图像的畸变。
3.特征提取：从图像中提取关键点或特征点。
这些特征点可以是角点、边缘、斑点等。
提取出的特征点用于后方交会计算。
4.求解相机位姿：使用特征点的像素坐标和已知内参数，通过解非线性方程组的方法，计算相机在三维空间中的位姿（即相机的位置和方向）。
5.求解三维点坐标：对于每个特征点，使用单像模型，将像素坐标投影到相机坐标系中。
然后，通过解线性方程组的方法，计算特征点在三维空间中的坐标。
6.误差检测与优化：计算测量误差，并进行误差检测。
可以使用一些优化算法，如最小二乘法，来优化相机位姿和三维点坐标。
7.输出测量结果：将计算得到的三维点坐标输出，可以是数字格式或者可视化结果。
以上是单像空间后方交会的基本程序设计步骤。
每个步骤可能会有不同的具体实现，根据具体的应用场景和需求进行设计和调整。
单像空间后方交会和双像解析空间后方-前方交会的算法程序实现遥感科学与技术摘要：如果已知每张像片的6个外方位元素，就能确定被摄物体与航摄像片的关系。
因此，利用单像空间后方交会的方法，可以迅速的算出每张像片的6个外方位元素。
而前方交会的计算，可以算出像片上点对应于地面点的三维坐标。
基于这两点，利用计算机强大的运算能力，可以代替人脑快速的完成复杂的计算过程。
关键词：后方交会，前方交会，外方位元素，C++编程0.引言：单张像片空间后方交会是摄影测量基本问题之一,是由若干控制点及其相应像点坐标求解摄站参数(X S,Y S,ZS,ψ、ω、κ)。
单像空间后方交会主要有三种方法:基于共线条件方程的平差解法、角锥法、基于直接线性变换的解法。
而本文将介绍第一种方法，基于共线条件方程反求象片的外方位元素。
而空间前方交会先以单张像片为单位进行空间后方交会，分别求出两张像片的外方位元素，再根据待定点的一对像点坐标，用空间前方交会的方法求解待定点的地面坐标。
可以说，这种求解地面点的坐标的方法是以单张像片空间后方交会为基础的，因此，单张像片空间后方交会成为解决这两个问题以及算法程序实现的关键。
1.单像空间后方交会的算法程序实现：（1）空间后方交会的基本原理：对于遥感影像,如何获取像片的外方位元素,一直是摄影测量工作者探讨的问题,其方法有:利用雷达(Radar)、全球定位系统(GPS)、惯性导航系统(I N S)以及星像摄影机来获取像片的外方位元素;也可以利用一定数量的地面控制点,根据共线方程,反求像片的外方位元素,这种方法称为单像空间后方交会(如图1所示)。
图中,地面坐标X i、Yi、Zi和对应的像点坐标x i、yi是已知的,外方位元素XS、Y S、ZS(摄站点坐标),ψ、ω、κ(像片姿态角)是待求的。
（2）空间后方交会数学模型：空间后方交会的数学模型是共线方程, 即中心投影的构像方程：式中X、Y、Z是地面某点在地面摄影测量坐标系中的坐标,x,y是该地面点在像片上的构像点的像片坐标,对于空间后方交会而言它们是已知的,还有主距f是已知的。
实验一近景单张像片空间后方交会一、实验目的通过单张像片空间后方交会算法编程，掌握空间后方交会的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。
二、实验内容利用C++编程平台，通过给定的单片像点的像点坐标、内方位元素和控制点物方空间坐标，计算出像片的外方位元素。
三、实验步骤：1、编程流程图：2、编程代码：#include "iostream"#include"stdio.h"#include "stdlib.h"#include<math.h>#define N 10using namespace std;void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2)//矩阵相乘{int i,j,k;for(i=0;i<i_1;i++)for(j=0;j<j_2;j++){result[i*j_2+j]=0.0;for(k=0;k<j_12;k++)result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];}return;}/* ######################矩阵求逆##################*/int invers_matrix(double *m1,int n){int *is,*js;int i,j,k,l,u,v;double temp,max_v;is=(int *)malloc(n*sizeof(int));js=(int *)malloc(n*sizeof(int));if(is==NULL||js==NULL){printf("out of memory!\n");return(0);}for(k=0;k<n;k++){max_v=0.0;for(i=k;i<n;i++)for(j=k;j<n;j++){temp=fabs(m1[i*n+j]);if(temp>max_v){max_v=temp; is[k]=i; js[k]=j;}}if(max_v==0.0){free(is); free(js);printf("invers is not availble!\n");return(0);}if(is[k]!=k)for(j=0;j<n;j++){u=k*n+j; v=is[k]*n+j;temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}if(js[k]!=k)for(i=0;i<n;i++){u=i*n+k; v=i*n+js[k];temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}l=k*n+k;m1[l]=1.0/m1[l];for(j=0;j<n;j++)if(j!=k){u=k*n+j;m1[u]*=m1[l];}for(i=0;i<n;i++)if(i!=k)for(j=0;j<n;j++)if(j!=k){u=i*n+j;m1[u]-=m1[i*n+k]*m1[k*n+j];}for(i=0;i<n;i++)if(i!=k){u=i*n+k;m1[u]*=-m1[l];}}for(k=n-1;k>=0;k--){if(js[k]!=k)for(j=0;j<n;j++){u=k*n+j; v=js[k]*n+j;temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}if(is[k]!=k)for(i=0;i<n;i++){u=i*n+k; v=i*n+is[k];temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;}}free(is); free(js);return(1);}/*------------------------------矩阵转置----------------------------------*/void transpose(double *m1,double *m2,int m,int n){int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++)m2[j*m+i]=m1[i*n+j];return;}void main(){double Xs,Ys,Zs,q,w,k;//外方位元素//double a[3],b[3],c[3];//旋转矩阵//double x0,y0,f;//内方位元素//double x[N],y[N];double X[N],Y[N],Z[N];double x1[N],y1[N];double m;//中误差//double L[2*N];double XX[6];double A[2*N][6];double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];//外方位元素初始置//int i,n=0;//迭代次数//double sum=0,m0;/*for(i=0;i<N;i++){printf("请输入第%d个点的地面坐标：",i+1); //----------输入点地面坐标-scanf_s("%lf%lf%lf",&X[i],&Y[i],&Z[i]);} */X[0]=1617.712 ;Y[0]=1883.870 ;Z[0]=-4917.308;X[1]=1615.572 ;Y[1]= 1687.514 ; Z[1]= -4913.750;X[2]=1614.640 ;Y[2]= 1315.066 ;Z[2]=-4916.724;X[3]=1615.784 ;Y[3]= 1001.358 ;Z[3]= -4916.691;X[4]=1617.197 ;Y[4]= 690.040 ;Z[4]= -4917.070;X[5]=1615.933 ;Y[5]= 387.492;Z[5]= -4917.234;X[6]=1618.946 ; Y[6]= 145.491;Z[6]= -4917.807;X[7]=1616.753 ;Y[7]= -186.640 ;Z[7]= -4918.531;X[8]=1617.714 ;Y[8]= -482.264;Z[8]= -4918.561;X[9]=1617.859 ;Y[9]= -680.767;Z[9]= -4918.799;/*for(i=0;i<N;i++){printf("请输入第%d个点的像片坐标：",i+1); //---------输入点像片坐标scanf_s("%lf%lf",&x[i],&y[i]);}cout<<endl;*/x[0]=1105.312313;y[0]=3367.979242;x[1]=1111.483438;y[1]=3203.998009;x[2]=1128.936009;y[2]=2886.146374;x[3]=1145.729428;y[3]=2615.385301;x[4]=1163.579521;y[4]=2343.603584;x[5]=1179.483899;y[5]=2076.672738;x[6]=1196.664025;y[6]=1861.971885;x[7]=1216.048311;y[7]=1565.046087;x[8]=1236.196347;y[8]=1299.275413;x[9]=1249.998721;y[9]=1120.52425;/*-----------------设定外方位元素初始值--------------*/x0=48.6401;y0=11.4765;f=4548.5949;Xs=3391.658;Ys=-135.645;Zs=97.250;q=-0.0102;w=0.0620;k=-0.0887;XX[3]=1;/*------------------迭代计算--------------------------*/while((XX[3]>0.00001 || XX[4]>0.00001 || XX[5]>0.00001)&&n<100) {/*----------------旋转矩阵R-----------------------*/a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);/*-----------------像点坐标计算值------------------*/for(i=0;i<N;i++){X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs);Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs);Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);x1[i]=x0-f*X0[i]/Z0[i];y1[i]=y0-f*Y0[i]/Z0[i];}/*-------------误差方程中各偏导数的值--------------*/for(i=0;i<N;i++){A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i];A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i];A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin( k))/f+f*cos(k))*cos(w);A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;A[2*i][5]=y[i]-y0;L[2*i]=x[i]-x1[i];A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];A[1+2*i][1]=((b[1]*f+b[2]*(y[i]-y0)))/Z0[i];A[1+2*i][2]=((c[1]*f+c[2]*(y[i]-y0)))/Z0[i];A[1+2*i][3]=-(x[i]-x0)*sin(w)-((y[i]-y0)*((x[i]-x0)*cos(k)-(y[i]-y0)*sin(k))/f-f*sin(k))*cos(w);A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k ))/f;A[1+2*i][5]=-x[i]+x0;L[1+2*i]=y[i]-y1[i];}/*-------------------解法方程--------------------*/transpose(&A[0][0],&At[0][0],2*N,6);mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);invers_matrix(&result1[0][0],6);mult(&At[0][0],L,&result2[0][0],6,2*N,1);mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);Xs+=XX[0];Ys+=XX[1];Zs+=XX[2];q+=XX[3];w+=XX[4];k+=XX[5];n++;}/*----------------旋转矩阵R-----------------------*/a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);cout<<"迭代次数为:"<<n<<endl;printf("\n像片外方位元素的解\n");cout<<"航向顷角q:"<<q<<endl;cout<<"旁向倾角w:"<<w<<endl;cout<<"像片旋角k:"<<k<<endl;printf("Xs=%lf\t\tYs=%lf\t\tZs=%lf\n",Xs,Ys,Zs);cout<<endl;printf("旋转矩阵R:\n");for(i=0;i<3;i++)printf("%lf\t",a[i]);printf("\n");for(i=0;i<3;i++)printf("%lf\t",b[i]);printf("\n");for(i=0;i<3;i++)printf("%lf\t",c[i]);printf("\n");/*-------------------计算单位权中误差---------------*/ for(i=0;i<2*N;i++)sum+=L[i]*L[i];m0=sqrt(sum/(2*N-6));cout<<"单位权中误差m0="<<m0<<endl;cout<<"测量精度:"<<endl;cout<<"Xs的精度为:"<<m0*sqrt(result1[0][0])<<endl;cout<<"Ys的精度为:"<<m0*sqrt(result1[1][1])<<endl;cout<<"Zs的精度为:"<<m0*sqrt(result1[2][2])<<endl;cout<<"q的精度为:"<<m0*sqrt(result1[3][3])<<endl;cout<<"w的精度为:"<<m0*sqrt(result1[4][4])<<endl;cout<<"k的精度为:"<<m0*sqrt(result1[5][5])<<endl;system("pause");}四、实验结果分析五、实验心得通过这次实验，使我掌握了在C++及Microsoft visual studio 环境下单张像片空间后方交会算法编程及实现，掌握空间后方交会的基本原理和基本步骤，并且锻炼了我编程的能力，看到自己编出的程序成功实现后，感到很有成就感，这也是我有了更大的自信，也为以后的工作打下了更好的基础，感谢老师的这次实习。
