空间后方交会程序
空间后方交会的解算
空间后方交会的解算一． 空间后方交会的目的摄影测量主要利用摄影的方法获取地面的信息，主要是是点位信息，属性信息，因此要对此进行空间定位和建模，并首先确定模型的参数，这就是空间后方交会的目的，用以求出模型外方位元素。
二． 空间后方交会的原理空间后方交会的原理是共线方程。
共线方程是依据相似三角形原理给出的，其形式如下111333222333()()()()()()()()()()()()A S A S A S A S A S A S AS A S A S A S A S A S a X X b Y Y c Z Z x f a X X a Y Y a Z Z a X X b Y Y c Z Z y f a X X a Y Y a Z Z -+-+-=--+-+--+-+-=--+-+-上式成为中心投影的构线方程，我们可以根据几个已知点，来计算方程的参数，一般需要六个方程，或者要三个点，为提高精度，可存在多余观测，然后利用最小二乘求其最小二乘解。
将公式利用泰勒公式线性化，取至一次项，得到其系数矩阵A ；引入改正数（残差）V ，则可将其写成矩阵形式：V AX L =-其中111333222333[,]()()()()()()()()()()()()()()Tx y A S A S A S x A S A S A S A S A S A S y A S A S A S L l l a X X b Y Y c Z Z l x x x fa X X a Y Y a Z Z a X Xb Y Yc Z Z l y y y fa X X a Y Y a Z Z =-+-+-=-=+-+-+--+-+-=-=+-+-+- 则1()T T X A A A L -=X 为外方位元素的近似改正数，由于采用泰勒展开取至一次项，为减少误差，要将的出的值作为近似值进行迭代，知道小于规定的误差三． 空间后方交会解算过程1. 已知条件近似垂直摄影00253.24mmx y 0f ===2. 解算程序流程图MATLAB 程序format long;s1=xlsread('data.xls');%读取数据a1=s1(1:4,1:2);%影像坐标b1=s1(1:4,3:5);%地面摄影测量坐标a2=s1.*10^-3;%影像坐标单位转化j1=a2(1,:)-a2(2,:);j2=j1(1,1)^2+j1(1,2)^2;lengh_a1=sqrt(j2); %相片某一长度j1=b1(1,:)-b1(1,:);j2=j1(1,1)^2+j1(1,2)^2;lengh_b1=sqrt(j2); %地面对应的长度m=lengh_b1/lengh_a1;%求出比例尺n0=0;p0=0;q0=0;x0=mean(b1(:,1));y0=mean(b1(:,2));f=153.24*10^-3;z0=m*f;x001={x0,x0,x0,x0};X0=cell2mat(x001)';y001={y0,y0,y0,y0};Y0=cell2mat(y001)';z001={z0,z0,z0,z0};Z0=cell2mat(z001)';%初始化外方位元素的值aa1=cos(n0)*cos(q0)-sin(n0)*sin(p0)*sin(q0);aa2=-sin(q0)*cos(n0)-sin(n0)*sin(p0)*cos(q0);aa3=-sin(n0)*cos(p0);bb1=sin(q0)*cos(p0);bb2=cos(q0)*cos(p0);bb3=-sin(p0);cc1=sin(n0)*cos(q0)+sin(p0)*cos(n0)*sin(q0);cc2=-sin(n0)*sin(q0)+sin(p0)*cos(q0)*cos(n0);cc3=cos(n0)*cos(p0);%计算改正数XX1=aa1.*(b1(:,1)-X0)+bb1.*(b1(:,2)-Y0)+cc1.*(b1(:,3)-Z0); XX2=aa2.*(b1(:,1)-X0)+bb2.*(b1(:,2)-Y0)+cc2.*(b1(:,3)-Z0); XX3=aa3.*(b1(:,1)-X0)+bb3.*(b1(:,2)-Y0)+cc3.*(b1(:,3)-Z0); lx=a1(:,1)+f.*(XX1./XX3);ly=a1(:,2)+f.*(XX2./XX3);l={lx',ly'};L=cell2mat(l)';%方程系数A=[-3.969*10^-5 0 2.231*10^-5 -0.2 -0.04 -0.06899;0 -3.969*10^-5 1.787*10^-5 -0.04 -0.18 0.08615;-2.88*10^-5 0 1*10^-5 -0.17 0.03 0.08211;0 -2.88*10^-5 -1.54*10^-5 0.03 -0.2 0.0534;-4.14*10^-5 0 4*10^-6 -0.15 -7.4*10^-3 -0.07663;0 -4.14*10^-5 2.07*10^-5 -7.4*10^-3 -0.19 0.01478;-2.89*10^-5 0 -1.98*10^-6 -0.15 -4.4*10^-3 0.06443;0 -2.89*10^-5 -1.22*10^-5 -4.4*10^-3 -0.18 0.01046];%L=[-1.28 3.78 -3.02 -1.45 -4.25 4.98 -4.72 -0.385]'.*10^-2; %第一次迭代X=inv(A'*A)*A'*L;3.结果X=1492.41127406195-554.4015671761941425.68660973544-0.0383847815608609 0.00911624039769785 -0.105416434087641S=1492.41127406195-554.401567176194 1425.68660973544 38436.9616152184 27963.1641162404-0.105416434087641。
空间后方交会程序设计
单片空间后方交会程序设计专业：测绘工程班级：测绘101姓名：张雯学号：2010220200512 用C 或VC++语言实现单片后方交汇的计算。
以单幅影像为基础，从该影象所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素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.00018182003 w=0.0021117837 精度为:0.00015959235 k=-0.067576934 精度为:7.2440432e-005 迭代次数：4Press any key to continue。
单像空间后方交会
像空间辅助坐标系(u,v,w)
坐标原点位于S，但坐标轴不一定与像平 面坐标轴平行，按需要定义。
像空坐标系与像空辅助坐标系之关系
物方坐标系
地面测量坐标系（Xt,Yt,Zt)
义。
地面摄影测量坐标系（X,Y,Z,)
原点位于地面某一已知点，坐标轴按需要定
地面测量坐标系与摄影测量坐标系之关系
确定像片相对S 的位置。 --焦距 --像主点 在像平面坐标系中 的坐标 例
外方位元素
1、确定S在物方空 间坐标系中位置的元 素（直线元素）。 Xs,Ys,Zs 例 Xs=1140.2m Ys=2003.5m Zs=1035.7m
பைடு நூலகம்2、确定像片在
物方空间坐标系中 位置的元素（角 元素）。 1） 角元素
像方坐标系与物方 坐标系之关系
共线方程线性化:
前式具体化:
即有
'
2
'
2
(5-9a)具体化:
写成
即
综合上述推导,有共线方程的线性形式:
式中
二.解算中的具体公式
利用(a)式解求外方位元素时,有6个未知数,须用像 片及地面3个点的3对已知的（X,Y,Z）、（x,y)组6个 方程.实用中为提高精度常取多余点多余观测,为此要按 最小二乘平差计算.则平差算式如下:
分）
单像空间后方交会（第五章部
根据单张航测像片上一定数量的已 知点(像片坐标和地面坐标已知),计算该 像片的外方位元素(摄影中心S的坐标 Xs,Ys,Zs,像片的角元素 ).
知道外方位 元素,可用来恢 复像片在摄影时 的空间位置,重 建像片与被摄地 面之间的相互关 系
内方位元素
( X1 , Y1 , Z1 )
后方交会法流程
后方交会法流程后方交会法呀，这可有点小意思呢。
后方交会法是一种在测量里挺有用的方法哦。
咱先说一下啥是后方交会法的基本原理吧。
它就是通过在未知点上设站，然后观测三个或者更多的已知控制点，根据这些观测数据来计算出这个未知点的坐标。
这就好比你在一个陌生的地方，你能看到周围几个你认识的标志性建筑，然后根据你看这些建筑的角度啊之类的信息，就能知道自己在啥位置啦。
那它的具体操作流程是啥样的呢？一、前期准备工作。
咱们得先准备好测量仪器呀，像全站仪这些，这就像是战士上战场要带好武器一样重要。
而且要保证仪器的状态良好，要是仪器出了问题，那后面测量出来的数据可就不准啦。
同时呢，还要拿到那些已知控制点的坐标资料，这可是关键的参考信息，没有这些，就像没了地图的旅行者，完全不知道方向呢。
二、设站。
把全站仪安置在这个未知点上，这个未知点的选择也有点小讲究哦。
要选择一个视野开阔的地方，这样才能很好地观测到那些已知控制点。
要是周围都是遮挡物，那可就没法好好测量啦。
在设站的时候呢，要把仪器调平，这一步可不能马虎，如果仪器没有调平，就像盖房子没有打好地基，后面的数据肯定是错得离谱。
三、观测已知控制点。
接下来就是观测那些已知控制点啦。
要仔细地测量每个控制点的水平角和垂直角之类的数据。
测量的时候要认真哦，多测几次取个平均值就更好啦，这样可以减少误差。
就像你做数学题，多检查几遍答案就更准确啦。
而且观测的时候要按照一定的顺序，这样不容易乱，也方便后面整理数据。
四、数据记录与整理。
把观测到的数据认真地记录下来，这可关系到最后的计算结果呢。
要是记录错了一个数字，那就像做饭放错了盐一样，整个味道就不对啦。
记录完之后呢，要检查一下有没有遗漏或者错误的地方。
然后根据这些数据开始计算未知点的坐标啦。
这个计算过程可能有点复杂，不过现在有很多软件可以帮助我们计算，但是咱也得知道计算的原理呀，不能完全依赖软件。
五、精度检查。
计算出坐标之后呢，可不能就这么完事儿了。
单像空间后方交会程序《JAVA版》
10. 将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。否则用新的近似值重 复 4~9 满足要求为止。 11. 由误差方程的计算式：
0 vx a11X s a12 Ys a13Z s a14 a15 a16 x x 0 v y a21X s a22 Ys a23Z s a24 a25 a26 y y
1. Java 程序 以下第一部分为主程序
import java.util.*; import java.math.*;
5. 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标，带入以下共线方程式，
x f y f a1 ( X A X S ) b1 (YA YS ) c1 ( Z A Z S ) a3 ( X A X S ) b3 (YA YS ) c3 ( Z A Z S ) a2 ( X A X S ) b2 (YA YS ) c2 ( Z A Z S ) a3 ( X A X S ) b3 (YA YS ) c3 ( Z A Z S )
以单像空间后方交会方法，求解该像片的外方位元素。 二、 计算原理
1. 获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为 地面摄测坐标。 2. 测量控制点的像点坐标并作系统误差改正。 3. 确定未知数的初始值。 在竖直摄影且地面控制点大体对称分布的情况下， 按如下方法确定初始值， 即
其他项类似可得。
x2 ) f2
y2 ) f2
由常数项计算公式：
a ( X X S ) b1 (Yi YS ) c1 ( Z i Z S ) xi f 1 i a3 ( X i X S ) b3 (Yi YS ) c3 ( Z i Z S ) lxi Li l yi y f a2 ( X i X S ) b2 (Yi YS ) c2 ( Z i Z S ) i a3 ( X i X S ) b3 (Yi YS ) c3 ( Z i Z S ) (i 1, 2,3, 4)
编写单片空间后方交会程序实习一
实验一 编写单片空间后方交会程序
一、实验目的
用程序设计语言（Visual C++或者C 语言、matlab 语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。
本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
二、实验内容与要求
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。
三、实验原理
根据摄影过程的几何反转原理。
四、实验数据准备
﹡ 已知航摄仪的内方位元素0,24.15300===y x mm f （摄影比例尺为1:50000）； ﹡ 4个地面控制点的地面坐标及其对应像点的像片坐标：
1．每人必须独立上机操作；
2．实验前认真复习单片空间后方交会的有关内容。
六、实验报告
原理、实验步骤及程序框图
上机调试程序并打印结果。
摄影测量学空间后方交会实验报告
摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘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.评定精度四、实验意义及能力培养正文：单像空间后方交会是一种基于摄影测量原理的算法，通过计算影像的外方位元素，实现对影像的定位和测量。
在已知地面上若干点的地面坐标和对应像点的像片坐标的情况下，利用共线条件方程求解像片外方位元素。
以下详细介绍单像空间后方交会的程序设计步骤：一、单像空间后方交会概念阐述单像空间后方交会是指在影像覆盖范围内，根据一定数量的分布合理的地面控制点（已知其像点和地面点的坐标），利用共线条件方程求解像片外方位元素的过程。
它是摄影测量中一种基础的算法，为后续复杂算法的演进提供了基础。
二、单像空间后方交会程序设计目的单像空间后方交会程序设计的目的是让学生深入理解单片空间后方交会的原理，通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
三、单像空间后方交会程序设计步骤1.数据准备：已知航摄仪的内方位元素，摄影比例尺，以及地面控制点的地面坐标和对应像点的像片坐标。
2.建立中心投影几何模型：根据针孔相机模型，建立中心投影的几何模型，包括物空间坐标系、像空间坐标系和摄影坐标系。
3.求解像片外方位元素：利用共线条件方程，通过最小二乘平差方法求解像片的外方位元素。
4.评定精度：根据实验数据，评定求解得到的像片外方位元素的精度。
四、实验意义及能力培养单像空间后方交会实验有助于学生掌握摄影测量基本原理，了解单像空间后方交会的计算过程，提高动手实践能力。
通过实验，学生可以深入理解线性代数和微分学在摄影测量中的应用，为后续学习提供更扎实的基础。
此外，实验还可以培养学生的解决问题的能力和综合运用所学知识的能力，为未来从事相关领域工作打下坚实基础。
综上所述，单像空间后方交会是一种重要的摄影测量方法，通过程序设计实现对影像的定位和测量。
摄影测量实验报告(空间后方交会—前方交会)
空间后方交会—空间前方交会程序编程实验一．实验目的要求掌握运用空间后方交会-空间前方交会求解地面点的空间位置.学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。
然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程.二．仪器用具计算机、编程软件(MATLAB）三．实验数据实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。
此四对坐标运用最小二乘法求解左右像片的外方位元素，即完成了空间后方的过程.另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。
实验数据如下：内方位元素：f=152。
000mm，x0=0，y0=0 四．实验框图此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0。
00003，相当于0。
1'的角度值)为止。
在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值,每次迭代所得的改正数都应加到上一次的初始值之中。
在空间后方交会中运用的数学模型为共线方程确定Xs，Ys，Zs的初始值时，对于左片可取地面左边两个GCP的坐标的平均值作为左片Xs 和Ys的初始值，取右边两个GCP的坐标平均值作为右片Xs 和Ys的初始值。
Zs可取地面所有GCP的Z坐标的平均值再加上航高.空间前方交会的数学模型为：五．实验源代码function Main_KJQHFJH（）global R g1 g2 m G a c b1 b2;m=10000；a=5;c=4；feval（＠shuru)；％调用shuru(）shurujcp（)函数完成像点及feval（＠shurujcp)；%CCP有关数据的输入XYZ=feval（＠MQZqianfangjh）; %调用MQZqianfangjh（）函数完成空间前方、%%%％%% 单位权中误差％%%％%后方交会计算解得外方位元素global V1 V2；%由于以上三个函数定义在外部文件中故需VV=[］; ％用feval（）完成调用过程for i=1:2＊cVV(i)=V1（i)；VV(2*i+1）=V2(i)；endm0=sqrt（VV＊（VV’)/（2＊c-6））；disp('单位权中误差m0为正负：’）；disp(m0）; ％计算单位权中误差并将其输出显示输入GCP像点坐标及地面摄影测量坐标系坐标的函数和输入所求点像点坐标函数：function shurujcp（)global c m;m=input（’摄影比例尺:')；%输入GCP像点坐标数据函数并分别将其c=input（'GCP的总数=')；％存入到不同的矩阵之中disp（'GCP左片像框标坐标：')；global g1;g1=zeros（c，2）;i=1;while i<=cm=input('x='）；n=input（'y='）;g1(i，1)=m;g1（i，2)=n；i=i+1；enddisp('GCP右片像框标坐标：’);global g2；g2=zeros(c,2);i=1;while i〈=cm=input（'x=’）;n=input('y=’）;g2（i，1)=m；g2（i，2）=n;i=i+1;end％%%%％%%%%%%%%％％%%％％%％％％％%％%%％％%%％％%％％%%％% function shuru（)global a；a=input（'计算总像对点数='); %完成想计算所需的像平面坐标global b1；％坐标输入,存入不同的矩阵中b1=zeros(a,2)；disp('左片像点坐标：'）i=1；while i〈=am=input（'x=’）；n=input(’y=’);b1(i，1)=m；b1（i,2)=n;i=i+1;end%%global b2;b2=zeros（a,2）；disp（’右片像点坐标：')i=1;while i〈=am=input('x=’）;n=input（'y=’）；b2(i,1）=m；b2(i,2)=n；i=i+1;end％%global c;c=input(’GCP的总数='）；disp('GCP摄影测量系坐标：’)global G;G=zeros(3,c)；i=1；while i〈=cm=input（’X='）;n=input(’Y='）;v=input（’Z=')；G(i，1）=m;G（i，2)=n;G（i,3)=v；i=i+1;end %%％％％％%%％%%%%%%％％％％％％％％%%％％%％%%%％％%%%空间前方交会和后方交会函数:function XYZ=MQZqianfangjh（）global R1 R2 a f b1 b2 Ra Rb；global X1 X2;R1=Ra；R2=Rb;R1=zeros（3,3）;R2=zeros（3,3)；global g1 g2 G V1 V2 V WF c QXX QXX1 QXX2；xs0=（G（1，1)+G(3,1))/2；ys0=（G（1，2）+G(3，2)）/2;［Xs1，Ys1，Zs1,q1,w1，k1 R]=houfangjh(g1，xs0,ys0）；％对左片调用后方交会函数R1=R；V1=V；WF1=WF;QXX1=QXX;save '左片外方位元素为。
摄影测量学单像空间后方交会程序设计作业
using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace 单像空间后方交会{class Program{static void Main(string[] args){int x0, y0, i, j; double f, m;Console.Write("请输入像片比例尺：");m = double.Parse(Console.ReadLine());Console.Write("请输入像片的内方位元素x0:");//均以毫米为单位x0 = int.Parse(Console.ReadLine());Console.Write("请输入像片的内方位元素y0:");y0 = int.Parse(Console.ReadLine());Console.Write("请输入摄影机主距f:");f = double.Parse(Console.ReadLine());Console.WriteLine();//输入坐标数据double[,] zuobiao = new double[4, 5];for (i = 0; i < 4; i++){for (j = 0; j < 5; j++){if (j < 3){Console.Write("请输入第{0}个点的第{1}个地面坐标：", i + 1, j + 1);zuobiao[i, j] =double.Parse(Console.ReadLine());}else{Console.Write("请输入第{0}个点的第{1}个像点坐标：", i + 1, j - 2);zuobiao[i, j] =double.Parse(Console.ReadLine());}} Console.WriteLine();}//归算像点坐标for (i = 0; i < 4; i++){for (j = 3; j < 5; j++){if (j == 3)zuobiao[i, j] = zuobiao[i, j] - x0;elsezuobiao[i, j] = zuobiao[i, j] - y0;}}//计算和确定初值double zs0 = m * f, xs0 = 0, ys0 = 0;for (i = 0; i < 4; i++){xs0 = xs0 + zuobiao[i, 0];ys0 = ys0 + zuobiao[i, 1];}xs0 = xs0 / 4;ys0 = ys0 / 4;//逐点计算误差方程系数double[,] xishu = new double[8, 6];for (i = 0; i < 8; i += 2){double x, y;x = zuobiao[i / 2, 3]; y = zuobiao[i / 2, 4];xishu[i, 0] = xishu[i + 1, 1] = -1 / m; xishu[i, 1] = xishu[i + 1, 0] = 0; xishu[i, 2] = -x / (m * f); xishu[i, 3] = -f * (1 + x * x / (f * f));xishu[i, 4] = xishu[i + 1, 3] = -x * y / f; xishu[i, 5] = y; xishu[i + 1, 2] = -y / (m * f); xishu[i + 1, 4] = -f * (1 + y * y / (f * f)); xishu[i + 1, 5] = -x;}//计算逆阵double[,] dMatrix =matrixChe(matrixTrans(xishu), xishu);double[,] dReturn = ReverseMatrix(dMatrix, 6);Console.WriteLine("逆矩阵为：");if (dReturn != null){matrixOut(dReturn);}//求解过程double phi0 = 0, omega0 = 0, kappa0 = 0; int q = 0;double[,] r = new double[3, 3];double[,] jinsi = new double[4, 2];double[] chazhi = new double[8];double[] jieguo = new double[6];double[,] zhong = matrixChe(dReturn,matrixTrans(xishu));do{ //计算旋转矩阵rr[0, 0] = Math.Cos(phi0) * Math.Cos(kappa0) - Math.Sin(phi0) * Math.Sin(omega0) * Math.Sin(kappa0);r[0, 1] = -Math.Cos(phi0) * Math.Sin(kappa0) - Math.Sin(phi0) * Math.Sin(omega0) * Math.Cos(kappa0);r[0, 2] = -Math.Sin(phi0) * Math.Cos(omega0);r[1, 0] = Math.Cos(omega0) * Math.Sin(kappa0);r[1, 1] = Math.Cos(omega0) * Math.Cos(kappa0);r[1, 2] = -Math.Sin(omega0);r[2, 0] = Math.Sin(phi0) * Math.Cos(kappa0) + Math.Cos(phi0) * Math.Sin(omega0) * Math.Sin(kappa0);r[2, 1] = -Math.Sin(phi0) * Math.Sin(kappa0) + Math.Cos(phi0) * Math.Sin(omega0) * Math.Cos(kappa0);r[2, 2] = Math.Cos(phi0) * Math.Cos(omega0);//计算x,y的近似值for (i = 0; i < 4; i++){jinsi[i, 0] = -f * (r[0, 0] * (zuobiao[i, 0] - xs0) + r[1, 0] * (zuobiao[i, 1] - ys0) + r[2, 0] * (zuobiao[i, 2] - zs0)) / (r[0, 2] * (zuobiao[i, 0] - xs0) + r[1, 2] * (zuobiao[i, 1] - ys0) + r[2, 2] * (zuobiao[i, 2] - zs0));jinsi[i, 1] = -f * (r[0, 1] * (zuobiao[i, 0] - xs0) + r[1, 1] * (zuobiao[i, 1] - ys0) + r[2, 1] * (zuobiao[i, 2] - zs0)) / (r[0, 2] * (zuobiao[i, 0] - xs0) + r[1, 2] * (zuobiao[i, 1] - ys0) + r[2, 2] * (zuobiao[i, 2] - zs0));}for (i = 0; i < 8; i += 2){chazhi[i] = zuobiao[i / 2, 3] - jinsi[i / 2, 0];chazhi[i + 1] = zuobiao[i / 2, 4] - jinsi[i / 2, 1];}for (i = 0; i < zhong.GetLength(0); i++){double k = 0;for (j = 0; j < zhong.GetLength(1); j++){k = k + zhong[i, j] * chazhi[j];}jieguo[i] = k;}//求新的近似值xs0 += jieguo[0]; ys0 += jieguo[1]; zs0 += jieguo[2];phi0 += jieguo[3]; omega0 += jieguo[4]; kappa0 += jieguo[5];q++;if (q > 1000)break;} while ((Math.Abs(jieguo[0]) > 0.020 ||Math.Abs(jieguo[1]) > 0.020) || Math.Abs(jieguo[2]) > 0.020);Console.WriteLine("共进行了{0}次运算", q);Console.WriteLine("旋转矩阵为");matrixOut(r);for (i = 0; i < jieguo.GetLength(0); i++){Console.Write("第{0}个外方位元素为：{1}", i + 1, jieguo[i]);}}//矩阵转置public static double[,] matrixTrans(double[,] X){double[,] A = X;double[,] C = new double[A.GetLength(1),A.GetLength(0)];for (int i = 0; i < A.GetLength(1); i++)for (int j = 0; j < A.GetLength(0); j++){C[i, j] = A[j, i];}return C;}//矩阵输出public static void matrixOut(double[,] X){double[,] C = X;for (int i = 0; i < C.GetLength(0); i++){for (int j = 0; j < C.GetLength(1); j++){Console.Write(" {0}", C[i, j]);}Console.Write("\n");}}//二维矩阵相乘public static double[,] matrixChe(double[,] X, double[,] Y){int i, j, n; double m;double[,] C = X; double[,] D = Y;double[,] E = new double[C.GetLength(0),C.GetLength(0)];for (i = 0; i < C.GetLength(0); i++){for (n = 0; n < C.GetLength(0); n++){m = 0;for (j = 0; j < C.GetLength(1); j++){m = m + C[i, j] * D[j, n];}E[i, n] = m;}}return E;}//计算行列式的值public static double MatrixValue(double[,] MatrixList, int Level){double[,] dMatrix = new double[Level, Level];for (int i = 0; i < Level; i++)for (int j = 0; j < Level; j++)dMatrix[i, j] = MatrixList[i, j];double c, x;int k = 1;for (int i = 0, j = 0; i < Level && j < Level; i++, j++){if (dMatrix[i, j] == 0){int m = i;for (; dMatrix[m, j] == 0; m++) ;if (m == Level)return 0;else{for (int n = j; n < Level; n++){c = dMatrix[i, n];dMatrix[i, n] = dMatrix[m, n];dMatrix[m, n] = c;}k *= (-1);}}for (int s = Level - 1; s > i; s--){x = dMatrix[s, j];for (int t = j; t < Level; t++)dMatrix[s, t] -= dMatrix[i, t] * (x / dMatrix[i, j]);}}double sn = 1;for (int i = 0; i < Level; i++){if (dMatrix[i, i] != 0)sn *= dMatrix[i, i];elsereturn 0;}return k * sn;}//计算逆阵public static double[,] ReverseMatrix(double[,] dMatrix, int Level){double dMatrixValue = MatrixValue(dMatrix, Level);if (dMatrixValue == 0) return null;double[,] dReverseMatrix = new double[Level, 2 * Level];double x, c;for (int i = 0; i < Level; i++){for (int j = 0; j < 2 * Level; j++){if (j < Level)dReverseMatrix[i, j] = dMatrix[i, j];elsedReverseMatrix[i, j] = 0;}dReverseMatrix[i, Level + i] = 1;}for (int i = 0, j = 0; i < Level && j < Level; i++, j++){if (dReverseMatrix[i, j] == 0){int m = i;for (; dMatrix[m, j] == 0; m++) ;if (m == Level)return null;else{for (int n = j; n < 2 * Level; n++)dReverseMatrix[i, n] += dReverseMatrix[m, n];}}x = dReverseMatrix[i, j];if (x != 1){for (int n = j; n < 2 * Level; n++)if (dReverseMatrix[i, n] != 0)dReverseMatrix[i, n] /= x;}for (int s = Level - 1; s > i; s--){x = dReverseMatrix[s, j];for (int t = j; t < 2 * Level; t++)dReverseMatrix[s, t] -= (dReverseMatrix[i, t] * x);}}for (int i = Level - 2; i >= 0; i--){for (int j = i + 1; j < Level; j++)if (dReverseMatrix[i, j] != 0){c = dReverseMatrix[i, j];for (int n = j; n < 2 * Level; n++)dReverseMatrix[i, n] -= (c * dReverseMatrix[j, n]);}}double[,] dReturn = new double[Level, Level];for (int i = 0; i < Level; i++)for (int j = 0; j < Level; j++)dReturn[i, j] = dReverseMatrix[i, j + Level];return dReturn;}}}。
简述单像空间后方交会的程序设计步骤
简述单像空间后方交会的程序设计步骤
单像空间后方交会是一种用于测量摄影点在三维空间中位置的方法。
以下是简述的程序设计步骤：
1.读取摄影测量数据：首先，从摄影测量设备（如相机）中读取图像和相关的内参数据，包括相机的焦距、像点大小等。
2.图像处理：对读取的图像进行预处理。
可能需要进行去畸变操作，校正图像的畸变。
3.特征提取：从图像中提取关键点或特征点。
这些特征点可以是角点、边缘、斑点等。
提取出的特征点用于后方交会计算。
4.求解相机位姿：使用特征点的像素坐标和已知内参数，通过解非线性方程组的方法，计算相机在三维空间中的位姿（即相机的位置和方向）。
5.求解三维点坐标：对于每个特征点，使用单像模型，将像素坐标投影到相机坐标系中。
然后，通过解线性方程组的方法，计算特征点在三维空间中的坐标。
6.误差检测与优化：计算测量误差，并进行误差检测。
可以使用一些优化算法，如最小二乘法，来优化相机位姿和三维点坐标。
7.输出测量结果：将计算得到的三维点坐标输出，可以是数字格式或者可视化结果。
以上是单像空间后方交会的基本程序设计步骤。
每个步骤可能会有不同的具体实现，根据具体的应用场景和需求进行设计和调整。
空间后方交会程序
一． 实验目的： 掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。
二． 仪器用具及已知数据文件： 计算机windows xp 系统，编程软件（VISUAL C++6.0），地面控制点在摄影测量坐标系中的坐标及其像点坐标文件shuju.txt 。
三． 实验内容：单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标根据共线方程反求影像的外方位元素。
数学模型：共线条件方程式： )(3)(3)(3)(1)(1)(1Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f x -+-+--+-+--= )(3)(3)(3)(2)(2)(2Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f y -+-+--+-+--= 求解过程： （1）获取已知数据。
从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影测量坐标。
（2）量测控制点的像点坐标并做系统改正。
（3）确定未知数的初始值。
在竖直摄影且地面控制点大致分布均匀的情况下，按如下方法确定初始值，即： n X X S ∑=0,n Y Y S ∑=0,n Z mf Z S ∑=0 φ =ω=κ=0 式中；m 为摄影比例尺分母；n 为控制点个数。
（4）用三个角元素的初始值，计算个方向余弦，组成旋转矩阵R 。
（5）逐点计算像点坐标的近似值。
利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x ）、（y ）。
（6）逐点计算误差方程式的系数和常数项，组成误差方程式。
（7）计算法方程的系数矩阵A A T 和常数项l A T ，组成法方程式。
（8）解法方程，求得外方位元素的改正数dXs ,S dY ，s dZ ，d φ,d ω,d κ。
（9）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。
单像空间后方交会和双像解析空间后方-前方交会的算法程序实现
单像空间后方交会和双像解析空间后方-前方交会的算法程序实现遥感科学与技术摘要：如果已知每张像片的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是已知的。
全站仪后方交会法步骤和高程测量步骤
全站仪后方交会法步骤和高程测量步骤全站仪是一种常用于测量地面高程和水平角度的仪器。
在工程测量中，经常会使用全站仪后方交会法进行高程测量。
下面将详细介绍全站仪后方交会法的步骤和高程测量的步骤。
1.设置仪器：首先，需要选择一个适合的测量点作为基准点，并将全站仪放置在基准点上。
将全站仪水平放置，并通过调整三个螺丝调整水平仪气泡位于中心位置。
然后，使用全站仪的目标板对准基准点。
2.测量目标点：使用全站仪的望远镜和交会杆，在目标点上设置目标板。
目标板上的标识点应与全站仪的十字线对齐。
准确平稳地在目标点上设置目标板。
3.观测目标点：通过调整全站仪的望远镜，使其对准目标板上的标识点。
在读数之前，要确保全站仪已经稳定下来。
然后，记录望远镜的水平角和垂直角的读数。
4.移动到下一个目标点：移动全站仪到下一个目标点，并重复步骤2和步骤3、在每次观测之间，全站仪应保持在基准点上，并使用目标板进行校准。
5.数据处理：利用观测到的水平角和垂直角的读数，可以计算出各个目标点之间的坐标和高程差。
这种计算可以使用后方交会法进行，根据目标点在水平方向和垂直方向上的角度差，以及目标点之间的距离差，推导出目标点的空间坐标。
高程测量步骤如下：1.设置起始点：选择一个起始点作为基准点。
全站仪被放置在基准点上，并确保仪器水平放置。
2.目标点设置：将目标板设置在需要测量高程的点上。
目标板上的标识点应与全站仪的十字线对齐。
3.观测目标点：调整全站仪的望远镜，使其对准目标板上的标识点。
在记录读数之前，要确保全站仪稳定下来。
然后，记录望远镜的垂直角的读数。
4.移动到下一个目标点：移动全站仪到下一个需要测量高程的点，并重复步骤2和步骤35.高程差计算：根据每个目标点的垂直角的读数，可以计算出不同目标点之间的高程差。
通过将起始点的高程与每个目标点的高程差相加，可以得到每个目标点的实际高程。
6.数据处理：将所有测量得到的目标点的实际高程整理并记录。
进行必要的校正和调整，以获得更准确的高程数据。
9-空间后方交会
a11 a12 Ai a21 a22
a13 a23
a14 a24
a15 a25
a16 a26
T
X dXS
dYS
dZS
d d d
li l x
ly
T
vi vx
vy
T
把所有像控点的误差方程列出后，构成总 误差方程，根据最小二乘间接评差原理可 列出法方程式：
五、空间后方交会实践
如何获取像片的六个外方位元素？
1)利用雷达；全球定位系统GPS；惯性导航系统.
2)空间后方交会：利用一定数量的地面控制点， 根据共线方程，反求像片的外方位元素。（已知 像片的内方位元素，至少三个地面点坐标并测出 相应的像点坐标）
计算要点：
1）计算的数学模型：共线方程按泰勒级数展开， 取一次项（线性化）。
2）在像片的四角选取四个或更多地面控制点，利 用最小二乘法平差计算。
计算步骤： 1）获取已知数据：比例尺1/m;H;内定向;控制点 地面坐标。
2）测量控制点的像点坐标。标刺，测量像框坐 标；像主点改正。 3）确定未知参数的初始值：竖直摄影时，角元 素初始值为零；线元素中，Zs0=H=mf; Xs0,Ys0取 控制点坐标的均值。
投 影 中 心 的 系 数
二、线性化-续
X Y Z
x xs 1 1 R R R y y s z zs
1 1
其中，R
R R R
1
1
1
把各偏导数代入整理得
f XX b2 Z Z ZZ x XX f sin XY f cos a 15 fsin ZZ ZZ x Yf a 16 Z fX f b1 YY f b2 XY y b 3 a 24 f b1 Z ZZ ZZ y XY f sin YY f cos a 25 fcos ZZ ZZ y X f a 26 Z x a 14 f Yf
空间后方交会C++程序代码
摄影测量后方交会程序（c/c++）输入数据截图：结果截图：程序源代码（其中的矩阵求逆在前面已经有了，链接）：#include <stdio.h>#include <stdlib.h>#include <math.h>const double PRECISION=1e-5;typedef double DOUBLE[5];int InputData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f);int InverseMatrix(double *matrix,const int &row);int main(int argc, char* argv[]){DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if(Resection(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if (Data!=NULL){delete []Data;}printf("解算完毕...\n");do{printf("计算结果保存于\"结果.txt\"文件中\n""请选择操作（输入P打开结果数据，R打开原始数据，其它退出程序）：");fflush(stdin); //刷新输入流char order=getchar();if ('P'==order || 'p'==order){system("结果.txt");}else if ('R'==order || 'r'==order){system("data.txt");}elsebreak;system("cls");}while(1);system("PAUSE");return 0;}/***********************************************函数名：InputData*函数介绍：从文件（data.txt）中读取数据,*文件格式如下：*点数 m(未知写作0)* 内方位元素(f x0 y0)*编号 x y X Y Z*下面是一个实例：4 0153.24 0 01 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31*参数：(in/out)Num(点数)，*(in/out)Data(存放数据),m,f,x0,y0*返回值：int ，0成功，1文件打开失败，2控制点个*数不足，3文件格式错误*作者：vcrs*完成时间：09-10-4**********************************************/int InputData(int &Num, DOUBLE *&Data,double &m,double &f) {double x0,y0;FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fscanf(fp_input,"%d%lf",&Num,&m);if (Num<4){return 2;}fscanf(fp_input,"%lf%lf%lf",&f,&x0,&y0);f/=1000;if (m<0 || f<0){return 3;}Data=new DOUBLE[Num];double *temp= new double[Num-1];double scale=0;int i;for (i=0;i<Num;i++){//读取数据，忽略编号if(fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&Data[i][0],&Data[i][1],&Data[i][2],&Data[i][3],&Data[i][4])!=5){return 3;}//单位换算成mData[i][0]/=1000.0;Data[i][1]/=1000.0;}//如果m未知则归算其值if (0==m){for (i=0;i<Num-1;i++){temp[i]=(Data[i][2]-Data[i+1][2])/(Data[i][0]-Data[i+1][0])+ (Data[i][3]-Data[i+1][3])/(Data[i][1]-Data[i+1][1]);scale+=temp[i]/2.0;}m=scale/(Num-1);}fclose(fp_input);delete []temp;return 0;}/***********************************************函数名：MatrixMul*函数介绍：求两个矩阵的积,*参数：Jz1(第一个矩阵),row(第一个矩阵行数),*Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个*矩阵列数),(out)JgJz(存放结果矩阵)*返回值：void*作者：vcrs*完成时间：09-10-4**********************************************/void MatrixMul(double *Jz1,const int &row,double *Jz2,const int &line,const int &com,double *JgJz){for (int i=0;i<row;i++){for (int j=0;j<line;j++){double temp=0;for (int k=0;k<com;k++){temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j));}*(JgJz+i*line+j)=temp;}}}/***********************************************函数名：OutPut*函数介绍：向结果.txt文件输出数据*参数：Q协因数阵，m精度，m0单位权中误差，6个外*方位元素，旋转矩阵*返回值：int，0成功，1失败*作者：vcrs*完成时间：09-10-4**********************************************/int OutPut(const double *&Q,const double *&m,const double &m0, const double &Xs,const double &Ys,const double &Zs,const double &Phi,const double &Omega,const double &Kappa,const double *R){FILE *fp_out;if (!(fp_out=fopen("结果.txt","w"))){return 1;}FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"\n空间后方交会程序(C\\C++)\n遥感信息工程学院\n班级：""00000\n学号：0000000\n姓名：vcrs\n\n");fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"已知数据：\n\n已知点数：");int num;double temp,x,y;fscanf(fp_input,"%d%lf",&num,&temp);fprintf(fp_out,"%d\n",num);fprintf(fp_out,"摄影比例尺（0表示其值位置）：");fprintf(fp_out,"%10.0lf\n",temp);fprintf(fp_out,"内方位元素(f x0 y0)：");fscanf(fp_input,"%lf%lf%lf",&temp,&x,&y);fprintf(fp_out,"%10lf\t%10lf\t%10lf\n",temp,x,y);for (int i=0;i<num;i++){double temp[5];fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&temp[0],&temp[1],&temp[2],&temp[3],&temp[4]);fprintf(fp_out,"%3d\t%10lf\t%10lf\t%10lf\t%10lf\t%10lf\n", i+1,temp[0],temp[1],temp[2],temp[3],temp[4]);}fclose(fp_input);fprintf(fp_out,"**************************************""**************************************""*********************************\n");fprintf(fp_out,"计算结果如下：\n\n外方位元素：\n"); fprintf(fp_out,"\tXs=%10lf\n",Xs);fprintf(fp_out,"\tYs=%10lf\n",Ys);fprintf(fp_out,"\tZs=%10lf\n",Zs);fprintf(fp_out,"\tPhi=%10lf\n",Phi);fprintf(fp_out,"\tOmega=%10lf\n",Omega);fprintf(fp_out,"\tKappa=%10lf\n\n",Kappa);fprintf(fp_out,"旋转矩阵：\n");for (i=0;i<3;i++){fprintf(fp_out,"\t");for (int j=0;j<3;j++){fprintf(fp_out,"%10lf\t",*(R+i*3+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n单位权中误差：%10lf\n\n",m0);fprintf(fp_out,"协因数阵：\n");for (i=0;i<6;i++){fprintf(fp_out,"\t");for (int j=0;j<6;j++){fprintf(fp_out,"%20lf\t",*(Q+i*6+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n外方位元素精度：");for (i=0;i<6;i++){fprintf(fp_out,"%10lf\t",m[i]);}fprintf(fp_out,"\n");fprintf(fp_out,"**************************************""**************************************""*********************************\n");fclose(fp_out);return 0;}/***********************************************函数名：Resection*函数介绍：计算*参数：Num(点数)，Data(数据)，m,,f(焦距),x0,y0*返回值：int，0成功，其它失败*作者：vcrs*完成时间：09-10-4**********************************************/int Resection(const int &Num,const DOUBLE *&Data,const double &m, const double &f){double Xs=0,Ys=0,Zs=0;int i,j;//设置初始值for (i=0;i<Num;i++){Xs+=Data[i][2];Ys+=Data[i][3];}Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);double R[3][3]={0.0};double *L=new double[2*Num];typedef double Double6[6];Double6 *A=new Double6[2*Num];double *AT=new double[2*Num*6];double *ATA=new double[6*6];double *ATL=new double[6];double *Xg=new double[6];//迭代计算do{//旋转矩阵R[0][0]=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa);R[0][1]=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa);R[0][2]=-sin(Phi)*cos(Omega);R[1][0]=cos(Omega)*sin(Kappa);R[1][1]=cos(Omega)*cos(Kappa);R[1][2]=-sin(Omega);R[2][0]=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa);R[2][1]=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa);R[2][2]=cos(Phi)*cos(Omega);for (i=0;i<Num;i++){double X=R[0][0]*(Data[i][2]-Xs)+R[1][0]*(Data[i][3]-Ys)+ R[2][0]*(Data[i][4]-Zs);double Y=R[0][1]*(Data[i][2]-Xs)+R[1][1]*(Data[i][3]-Ys)+ R[2][1]*(Data[i][4]-Zs);double Z=R[0][2]*(Data[i][2]-Xs)+R[1][2]*(Data[i][3]-Ys)+ R[2][2]*(Data[i][4]-Zs);double xxx,yyy;xxx=-f*X/Z;yyy=-f*Y/Z;//常数项L[2*i]=Data[i][0]-(-f*X/Z);L[2*i+1]=Data[i][1]-(-f*Y/Z);A[2*i][0]=(R[0][0]*f+R[0][2]*(xxx))/Z;A[2*i][1]=(R[1][0]*f+R[1][2]*(xxx))/Z;A[2*i][2]=(R[2][0]*f+R[2][2]*(xxx))/Z;A[2*i][3]=(yyy)*sin(Omega)-(((xxx)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))+f*cos(Kappa))*cos(Omega);A[2*i][4]=-f*sin(Kappa)-((xxx)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i][5]=(yyy);A[2*i+1][0]=(R[0][1]*f+R[0][2]*(yyy))/Z;A[2*i+1][1]=(R[1][1]*f+R[1][2]*(yyy))/Z;A[2*i+1][2]=(R[2][1]*f+R[2][2]*(yyy))/Z;A[2*i+1][3]=-(xxx)*sin(Omega)-(((yyy)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))-f*sin(Kappa))*cos(Omega);A[2*i+1][4]=-f*cos(Kappa)-((yyy)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i+1][5]=-(xxx);}//求矩阵A的转置矩阵ATfor (i=0;i<2*Num;i++){for (j=0;j<6;j++){*(AT+j*2*Num+i)=A[i][j];}}//求ATAMatrixMul(AT,6,&A[0][0],6,2*Num,ATA);if(InverseMatrix(ATA,6))return 1;MatrixMul(AT,6,L,1,2*Num,ATL);MatrixMul(ATA,6,ATL,1,6,Xg);Xs+=Xg[0];Ys+=Xg[1];Zs+=Xg[2];Phi+=Xg[3];Omega+=Xg[4];Kappa+=Xg[5];} while(fabs(Xg[0])>=PRECISION ||fabs(Xg[1])>=PRECISION || fabs(Xg[2])>=PRECISION ||fabs(Xg[3])>=PRECISION ||fabs(Xg[4])>=PRECISION || (Xg[5])>=PRECISION);//注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值，//由于变换很小忽略double *Q=ATA;double *V=new double[2*Num];MatrixMul(&A[0][0],2*Num,Xg,1,6,V);double VTV=0;for(i=0;i<2*Num;i++){V[i]-=L[i];VTV+=V[i]*V[i];}double m0=sqrt(VTV/(2*Num-6));double *mm=new double[6];for (i=0;i<6;i++){mm[i]=sqrt(*(Q+i*6+i))*m0;}OutPut(Q,mm,m0,Xs,Ys,Zs,Phi,Omega,Kappa,&R[0][0]);delete []L;delete []A;delete []AT;delete []ATA;delete []ATL;delete []Xg;delete []mm;delete []V;return 0;}void swap(double &a,double &b){double temp=a;a=b;b=temp;}/***********************************************函数名：InverseMatrix*函数介绍：求矩阵的逆（高斯-约当法）*输入参数：(in/out)matrix（矩阵首地址），*（in）row（矩阵阶数）*输出参数：matrix（原矩阵的逆矩阵）*返回值：int ，0成功，1失败*调用函数：swap(double&,double&)*作者：vcrs*完成时间：09-10-4**********************************************/int InverseMatrix(double *matrix,const int &row) {double *m=new double[row*row];double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;i<row;i++)for (j=0;j<row;j++){*pt=*ptemp;ptemp++;pt++;}}int k;int *is=new int[row],*js=new int[row];for (k=0;k<row;k++){double max=0;//全选主元//寻找最大元素for (i=k;i<row;i++){for (j=k;j<row;j++){if (fabs(*(m+i*row+j))>max){max=*(m+i*row+j);is[k]=i;js[k]=j;}}}if (0 == max){return 1;}//行交换if (is[k]!=k){for (i=0;i<row;i++){swap(*(m+k*row+i),*(m+is[k]*row+i));}}//列交换if (js[k]!=k){for (i=0;i<row;i++){swap(*(m+i*row+k),*(m+i*row+js[k]));}}*(m+k*row+k)=1/(*(m+k*row+k));for (j=0;j<row;j++){if (j!=k){*(m+k*row+j)*=*((m+k*row+k));}}for (i=0;i<row;i++){if (i!=k){for (j=0;j<row;j++){if(j!=k){*(m+i*row+j)-=*(m+i*row+k)**(m+k*row+j);}}}}for (i=0;i<row;i++){if(i!=k){*(m+i*row+k)*=-(*(m+k*row+k));}}}int r;//恢复行列for (r=row-1;r>=0;r--){if (js[r]!=r){for (j=0;j<row;j++){swap(*(m+r*row+j),*(m+js[r]*row+j));}}if (is[r]!=r){for (i=0;i<row;i++){swap(*(m+i*row+r),*(m+i*row+is[r]));}}}ptemp=matrix;pt=m;for (i=0;i<row;i++){for (j=0;j<row;j++){*ptemp=*pt;ptemp++;pt++;}}delete []is;delete []js;delete []m;return 0;}。
C语言空间后方交会源代码
#include<stdio.h>#include<math.h>#define n 4 //控制点个数#define PI 3.14159265struct coordinate{double x; //像点坐标double y;double Xt; //控制点坐标double Yt;double Zt;};// void inverse(double c[6][6]) //矩阵求逆// {// int i,j,h,k;// double p;// double q[6][12];// for(i=0;i<6;i++)//构造高斯矩阵// for(j=0;j<6;j++)// q[i][j]=c[i][j];// for(i=0;i<6;i++)// for(j=6;j<12;j++)// {// if(i+6==j)// q[i][j]=1;// else// q[i][j]=0;// }// for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据// for(i=k+1;i<6;i++)// {// if(q[i][h]==0)// continue;// p=q[k][h]/q[i][h];// // p=q[i][h]/q[k][h];// for(j=0;j<12;j++)// {// q[i][j]*=p;// q[i][j]-=q[k][j];// }// }// for(h=k=5;k>0;k--,h--) // 消去对角线以上的数据// for(i=k-1;i>=0;i--)// {// if(q[i][h]==0)// continue;// p=q[k][h]/q[i][h];// // p=q[i][h]/q[k][h];// for(j=11;j>0;j--)// {// q[i][j]*=p;// q[i][j]-=q[k][j];// }// }// for(i=0;i<6;i++)//将对角线上数据化为1// {// p=1.0/q[i][i];// for(j=0;j<12;j++)// q[i][j]*=p;// }// for(i=0;i<6;i++) //提取逆矩阵// for(j=0;j<n;j++)// {// c[i][j]=q[i][j+6];// }// }void ContraryMatrix(double *const pMatrix, double *const _pMatrix, const int &dim) {int ii,jj,kk;int flag=0;double *tMatrix = new double[2*dim*dim];for (int i=0; i<dim; i++){for (int j=0; j<dim; j++)tMatrix[i*dim*2+j] = pMatrix[i*dim+j];}for (i=0; i<dim; i++){for (int j=dim; j<dim*2; j++)tMatrix[i*dim*2+j] = 0.0;tMatrix[i*dim*2+dim+i] = 1.0;}//Initialization over!for (i=0; i<dim; i++)//Process Cols{double base = tMatrix[i*dim*2+i];if (fabs(base) < 1E-300){for (ii=i;ii<dim;ii++){if(tMatrix[ii*dim*2+i]!=0){flag=1;for (jj=0;jj<2*dim;jj++){tMatrix[i*dim*2+jj]+=tMatrix[ii*dim*2+jj];}}}base = tMatrix[i*dim*2+i];if (flag==0){printf( "求逆矩阵过程中被零除，无法求解!") ;}// exit(0);}for (int j=0; j<dim; j++)//row{if (j == i) continue;double times = tMatrix[j*dim*2+i]/base;for (int k=0; k<dim*2; k++)//col{tMatrix[j*dim*2+k] = tMatrix[j*dim*2+k] - times*tMatrix[i*dim*2+k];}}for (int k=0; k<dim*2; k++){tMatrix[i*dim*2+k] /= base;}}for (i=0; i<dim; i++){for (int j=0; j<dim; j++)_pMatrix[i*dim+j] = tMatrix[i*dim*2+j+dim];}delete[] tMatrix;}void main(){double f,xo,yo; //内方位元素int m; //摄影比例尺分母f=0.15324;xo=0;yo=0;m=50000;struct coordinate coor[n]={{-0.08615,-0.06899,36589.41,25273.32,2195.17},{-0.05340,0.08221,37631.08,31324.51, 728.69},{-0.01478,-0.07663,39100.97,24934.98,2386.50},{0.01046,0.06443,40426.54,30319.81,757.31}};int i,j; //循环变量double Xs,Ys,Zs; //外方位线元素double phi,omega,kappa; //外方位角元素double R[3][3]; //旋转矩阵Rdouble (x)[n],(y)[n]; //控制点像点坐标的近似值double A[8][6]; //误差方程中的矩阵Adouble ATA[6][6]; //矩阵A的转置矩阵与A的乘积double L[8][1];double ATL[6][1]; //矩阵A的转置矩阵与L的乘积double _ATA[6][6];double X[6][1]; //未知数矩阵double V[8][1]; //改正数矩阵double DXs,DYs,DZs,Dphi,Domega,Dkappa; //6个外方位元素的改正值//与精度计算有关的量double m0; //单位权中误差double Q[6][6]; //协方差矩阵double m1,m2,m3,m4,m5,m6; //未知数Xs,Ys,Zs,phi,omega,kappa的中误差矩阵Xs=0;Ys=0;for(i=0;i<n;i++){Xs+=coor[i].Xt;Ys+=coor[i].Yt;}//外方位元素的初始值Xs/=n;Ys/=n;Zs=f*m;phi=0;omega=0;kappa=0; //在竖直摄影情况下do{//计算旋转矩阵R的各个元素R[0][0]=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);R[0][1]=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);R[0][2]=-sin(phi)*cos(omega);R[1][0]=cos(omega)*sin(kappa);R[1][1]=cos(omega)*cos(kappa);R[1][2]=-sin(omega);R[2][0]=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);R[2][1]=-sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa);R[2][2]=cos(phi)*cos(omega);//将未知数的近似值代人共线方程式，计算(x),(y)for(i=0;i<n;i++){(x)[i]=-f*(R[0][0]*(coor[i].Xt-Xs)+R[1][0]*(coor[i].Yt-Ys)+R[2][0]*(coor[i].Zt-Zs))/(R[0][2]*(coor[ i].Xt-Xs)+R[1][2]*(coor[i].Yt-Ys)+R[2][2]*(coor[i].Zt-Zs));(y)[i]=-f*(R[0][1]*(coor[i].Xt-Xs)+R[1][1]*(coor[i].Yt-Ys)+R[2][1]*(coor[i].Zt-Zs))/(R[0][2]*(coor[i].Xt -Xs)+R[1][2]*(coor[i].Yt-Ys)+R[2][2]*(coor[i].Zt-Zs));}//计算矩阵A的各个元素for(i=0;i<n;i++){A[2*i][0]=(R[0][0]*f + R[0][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][1]=(R[1][0]*f + R[1][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][2]=(R[2][0]*f + R[2][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][3]=(coor[i].y-yo)*sin(omega)-((coor[i].x-xo)/f*((coor[i].x-xo)*cos(kappa)-(coor[i].y-yo)* sin(kappa))+f*cos(kappa))*cos(omega);A[2*i][4]=-f*sin(kappa)-(coor[i].x-xo)/f*((coor[i].x-xo)*sin(kappa)+(coor[i].y-yo)*cos(kappa));A[2*i][5]=coor[i].y-yo;A[2*i+1][0]=(R[0][1]*f + R[0][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][1]=(R[1][1]*f + R[1][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][2]=(R[2][1]*f + R[2][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][3]=-(coor[i].x-xo)*sin(omega)-((coor[i].y-yo)/f*((coor[i].x-xo)*cos(kappa)-(coor[i].y-y o)*sin(kappa))-f*sin(kappa))*cos(omega);A[2*i+1][4]=-f*cos(kappa)-(coor[i].y-yo)/f*((coor[i].x-xo)*sin(kappa)+(coor[i].y-yo)*cos(kappa ));A[2*i+1][5]=-(coor[i].x-xo);}for (i=0;i<6;i++){for (j=0;j<6;j++){printf("%f ",A[i][j]);if (j%5==0&&j!=0){printf("\n");}}}printf("'\n");//计算ATA的各个元素for(i=0;i<6;i++)for(j=0;j<6;j++)ATA[i][j]=A[0][i]*A[0][j]+A[1][i]*A[1][j]+A[2][i]*A[2][j]+A[3][i]*A[3][j]+A[4][i]*A[4][j]+A[5][i]*A[5][j]+A[6][i]*A[6][j]+A[7][i]*A[7][j];for (i=0;i<6;i++){for (j=0;j<6;j++){_ATA[i][j]=ATA[i][j];printf("%f ",ATA[i][j]);if (j%5==0&&j!=0){printf("\n");}}}printf("\n");//计算矩阵L的各个元素for(i=0;i<n;i++){L[2*i][0]=coor[i].x+f*(R[0][0]*(coor[i].Xt-Xs)+R[1][0]*(coor[i].Yt-Ys)+R[2][0]*(coor[i].Zt-Zs))/(R [0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));L[2*i+1][0]=coor[i].y+f*(R[0][1]*(coor[i].Xt-Xs)+R[1][1]*(coor[i].Yt-Ys)+R[2][1]*(coor[i].Zt-Zs))/(R[0 ][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));}//计算矩阵ATL的各个元素值for(i=0;i<6;i++)ATL[i][0]=A[0][i]*L[0][0]+A[1][i]*L[1][0]+A[2][i]*L[2][0]+A[3][i]*L[3][0]+A[4][i]*L[4][0]+A[5][i]*L[5][0]+A[6][i]*L[6][0]+A[7][i]*L[7][0];/* for (i=0;i<6;i++){for (j=0;j<1;j++){ATL[i][j]=0;for (int k=0;k<8;k++){ATL[i][j]+=A[k][i]*L[k][j];}}}*///调用函数计算矩阵ATA的逆矩阵// inverse(ATA);ContraryMatrix(*_ATA, *ATA, 6);for (i=0;i<6;i++){for (j=0;j<6;j++){printf("%f ",ATA[i][j]);if (j%5==0&&j!=0){printf("\n");}}}//计算矩阵X的各个元素值for(i=0;i<6;i++){X[i][0]=ATA[i][0]*ATL[0][0]+ATA[i][1]*ATL[1][0]+ATA[i][2]*ATL[2][0]+ATA[i][3]*ATL[3][0] +ATA[i][4]*ATL[4][0]+ATA[i][5]*ATL[5][0];printf("%f ",X[i][0]);}// for (i=0;i<6;i++)// {// for (j=0;j<1;j++)// {// X[i][j]=0;// for (int k=0;k<6;k++)// {// X[i][j]+=ATA[i][k]*ATL[k][j];//// }// printf("%f",X[i][j]);// }// }DXs=X[0][0];DYs=X[1][0];DZs=X[2][0];Dphi=X[3][0];Domega=X[4][0];Dkappa=X[5][0];Xs+=DXs;Ys+=DYs;Zs+=DZs;phi+=Dphi;omega+=Domega;kappa+=Dkappa;//计算矩阵V的各个元素for(i=0;i<n;i++){V[2*i][0]=A[2*i][0]*DXs+A[2*i][1]*DYs+A[2*i][2]*DZs+A[2*i][3]*Dphi+A[2*i][4]*Domega+A[ 2*i][5]*Dkappa-L[2*i][0];V[2*i+1][0]=A[2*i+1][0]*DXs+A[2*i+1][1]*DYs+A[2*i+1][2]*DZs+A[2*i+1][3]*Dphi+A[2*i+1][ 4]*Domega+A[2*i+1][5]*Dkappa-L[2*i+1][0];}/* //像点坐标改正后的值for(i=0;i<n;i++){coor[i].x+=V[2*i][0];coor[i].y+=V[2*i+1][0];}*/}//判断限差的条件while(!(fabs(Dphi) < 0.1/60/180*PI && fabs(Domega) < 0.1/60/180*PI && fabs(Dkappa) < 0.1/60/180*PI));//外方位元素计算完毕//有关精度的计算for(i=0;i<6;i++)Q[i][i]=ATA[i][i];m0=sqrt((V[0][0]*V[0][0]+V[1][0]*V[1][0]+V[2][0]*V[2][0]+V[3][0]*V[3][0]+V[4][0]*V[4][0]+V [5][0]*V[5][0]+V[6][0]*V[6][0]+V[7][0]*V[7][0])/(2*n-6));//计算各未知数的中误差m1=sqrt(Q[0][0])*m0;m2=sqrt(Q[1][1])*m0;m3=sqrt(Q[2][2])*m0;m4=sqrt(Q[3][3])*m0;m5=sqrt(Q[4][4])*m0;m6=sqrt(Q[5][5])*m0;printf("改正后的相对坐标及其对应的地面点坐标（x,y,Xt,Yt,Zt）:\n");for(i=0;i<n;i++){printf("%lf\t%lf\t%lf\t%lf\t%lf",(x)[i],(y)[i],coor[i].Xt,coor[i].Yt,coor[i].Zt);printf("\n");}printf("旋转矩阵R：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%lf\t",R[i][j]);printf("\n");}printf("外方位元素(Xs,Ys,Zs,phi,omega,kappa)及改正值：\n");printf("%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n",Xs,DXs,Ys,DYs,Zs,DZs, phi,Dphi,omega,Domega,kappa,Dkappa);printf("单位权中误差:\n");printf("%0.9f\n",m0);printf("外方位元素(Xs,Ys,Zs,phi,omega,kappa)的中误差:\n");printf("%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n",m1,m2,m3,m4,m5,m6);//计算完毕，输出结果,以文件方式保存printf("计算结果存储在文件中\n");FILE *fp;fp=fopen("计算结果.txt","w");fprintf(fp,"改正后的相对坐标及其对应的地面点坐标（x,y,Xt,Yt,Zt）:\n");for(i=0;i<n;i++){fprintf(fp,"%lf\t%lf\t%lf\t%lf\t%lf",coor[i].x,coor[i].y,coor[i].Xt,coor[i].Yt,coor[i].Zt);fprintf(fp,"\n");}fprintf(fp,"旋转矩阵R：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)fprintf(fp,"%lf\t",R[i][j]);fprintf(fp,"\n");}fprintf(fp,"外方位元素(Xs,Ys,Zs,phi,omega,kappa)及改正值：\n");fprintf(fp,"%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n",Xs,DXs,Ys,DYs,Zs,D Zs,phi,Dphi,omega,Domega,kappa,Dkappa);fprintf(fp,"单位权中误差:\n");fprintf(fp,"%0.9f\n",m0);fprintf(fp,"外方位元素(Xs,Ys,Zs,phi,omega,kappa)的中误差:\n");fprintf(fp,"%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n",m1,m2,m3,m4,m5,m6);fclose(fp);}。
单片空间后方交会matlab程序
单片空间后方交会matlab程序f=0.15;fai=0;omi=0;ka=0;Xd=[5083.205,5780.02,5210.879,5909.264];Yd=[5852.099,5906.365,4258.446,4314.283];Zd=[527.925,571.549,461.81,455.484];%以上为输入gcp1到gcp4的地面坐标xx=[-0.07393,-0.005252,-0.079122,-0.009887];yx=[0.078705,0.078184,-0.078879,-0.080089];%右片像点坐标Xs=(5083.205+5780.02+5210.879+5909.264)/4;Ys=(5852.099+5906.365+4258.446+4314.283)/4;Sx=sqrt((xx(2)-xx(1))^2+(yx(2)-yx(1))^2);Sd=sqrt((Xd(2)-Xd(1))^2+(Yd(2)-Yd(1))^2);m=Sd/SxZs=m*f+1/4*(Zd(1)+Zd(2)+Zd(3)+Zd(4));X=[1;1;1;1;1;1];while 1a1=cos(fai)*cos(ka)-sin(fai)*sin(omi)*sin(ka);a2=-cos(fai)*sin(ka)-sin(fai)*sin(omi)*cos(ka);a3=-sin(fai)*cos(omi);b1=cos(omi)*sin(ka);b2=cos(omi)*cos(ka);b3=-sin(omi);c1=sin(fai)*cos(ka)+cos(fai)*sin(omi)*sin(ka);c2=-sin(fai)*sin(ka)+cos(fai)*sin(omi)*cos(ka);c3=cos(fai)*cos(omi);R=[a1,b1,c1;a2,b2,c2;a3,b3,c3]for n=1:1:4s1=[Xd(n)-Xs,Yd(n)-Ys,Zd(n)-Zs]';X_(n)=[a1,b1,c1]*s1;Y_(n)=[a2,b2,c2]*s1;Z_(n)=[a3,b3,c3]*s1;x(n)=-f*X_(n)/Z_(n);y(n)=-f*Y_(n)/Z_(n);end%以上循环用于求解像空间坐标x(1),y(1)到x(4),y(4)for p=1:1:4a11(p)=1/Z_(p)*(a1*f+a3*x(p));a12(p)=1/Z_(p)*(b1*f+b3*x(p)); a13(p)=1/Z_(p)*(c1*f+c3*x(p));a21(p)=1/Z_(p)*(a2*f+a3*y(p));a22(p)=1/Z_(p)*(b2*f+b3*y(p)); a23(p)=1/Z_(p)*(c2*f+c3*y(p));a14(p)=y(p)*sin(omi)-(x(p)/f*(x(p)*cos(ka)-y(p)*sin(ka))+f*cos(ka))*cos(omi);a15(p)=-f*sin(ka)-x(p)/f*(x(p)*sin(ka)+y(p)*cos(ka)); %计算A 矩阵a16(p)=y(p);a24(p)=-x(p)*sin(omi)-(y(p)/f*(x(p)*cos(ka)-y(p)*sin(ka))-f*sin(ka))*cos(omi);a25(p)=-f*cos(ka)-y(p)/f*(x(p)*sin(ka)+y(p)*cos(ka));a26(p)=-x(p) ;%右片endA=[a11(1),a12(1),a13(1),a14(1),a15(1),a16(1);a21(1),a22(1),a2 3(1),a24(1),a25(1),a26(1);a11(2),a12(2),a13(2),a14(2),a15(2),a16(2 );a21(2),a22(2),a23(2),a24(2),a25(2),a26(2);a11(3),a12(3),a13(3),a1 4(3),a15(3),a16(3);a21(3),a22(3),a23(3),a24(3),a25(3),a26(3);a11(4 ),a12(4),a13(4),a14(4),a15(4),a16(4);a21(4),a22(4),a23(4),a24(4),a2 5(4),a26(4)]lx=xx-x;ly=yx-y;%右片L=[lx(1),ly(1),lx(2),ly(2),lx(3),ly(3),lx(4),ly(4)]'X=inv(A'*A)*A'*Ldx=X(1,1);dy=X(2,1);dz=X(3,1);df=X(4,1);do=X(5,1);dk=X(6,1 );Xs=Xs+dx;Ys=dy+Ys ;Zs=dz+Zs;fai=fai+df;omi=omi+do;ka=ka+dk;if all(abs(X)<0.0003)break;endendV=A*X-LXs=Xs+dxYs=dy+YsZs=dz+Zsfai=fai+dfomi=omi+doka=ka+dkVV=0;for n=1:8VV=VV+V(n,1)*V(n,1);endQii=inv(A'*A)m0=sqrt(VV/2)mi=m0*sqrt(Qii)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
空间后方交会程序————————————————————————————————作者：————————————————————————————————日期：ﻩ一． 实验目的:掌握摄影测量空间后方交会的原理,利用计算机编程语言实现空间后方交会外方位元素的解算。
二． 仪器用具及已知数据文件：计算机wind ｏws xp 系统,编程软件(VI ＳＵA Ｌ C ＋+6．0)，地面控制点在摄影测量坐标系中的坐标及其像点坐标文件shu ｊu.txt 。
三． 实验内容:单张影像的空间后方交会:利用已知地面控制点数据及相应像点坐标根据共线方程反求影像的外方位元素。
数学模型：共线条件方程式:)(3)(3)(3)(1)(1)(1Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f x -+-+--+-+--=)(3)(3)(3)(2)(2)(2Zs Z c Ys Y b Xs X a Zs Z c Ys Y b Xs X a f y -+-+--+-+--= 求解过程:（1)获取已知数据。
从航摄资料中查取平均航高与摄影机主距;获取控制点的地面测量坐标并转换为地面摄影测量坐标。
（2)量测控制点的像点坐标并做系统改正。
(３)确定未知数的初始值。
在竖直摄影且地面控制点大致分布均匀的情况下,按如下方法确定初始值,即:n X X S ∑=0,n Y Y S ∑=0,nZ mf Z S ∑=0φ =ω=κ=0式中;ｍ为摄影比例尺分母;ｎ为控制点个数。
（４)用三个角元素的初始值,计算个方向余弦，组成旋转矩阵R 。
(5)逐点计算像点坐标的近似值。
利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x ）、（y ）。
（6)逐点计算误差方程式的系数和常数项,组成误差方程式。
(7）计算法方程的系数矩阵A A T 和常数项l A T ,组成法方程式。
（8）解法方程，求得外方位元素的改正数dXs ,S dY ，s dZ ,d φ,ｄω，d κ。
（９）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。
(10)将求得的外方位元素改正数与规定的限差比较，若小于限差则迭代结束。
否则用新的近似值重复（4)~(９)，直到满足要求为止。
四．实验结果：程序的源代码如下所示:#iｎcｌude<stdiｏ.h>＃include<sｔdliｂ.h>#iｎclude＜ｍaｌloc．h>#inｃlude<ｍaｔｈ.h＞#ｉncｌude<coｎiｏ.h＞#ｄefｉｎe N 4voｉd turｎ(doｕble *Ａ,dｏｕble Ａ2[],int ｍ,int n）//计算矩阵的转置｛int i，ｊ;ﻩfｏｒ（i=0;i<m;ｉ+＋）for(j=0;j<n;j+＋)ﻩA2[ｊ*m＋i］=A[i＊ｎ＋j];｝vｏid mulAB(doｕble *Ａ,ｄoublｅ*Ｂ,ｄouｂle ＊C,inｔａm,int ａｎ,ｉnｔｂm，int bn) ／/计算两矩阵相乘{int i，j,l,u;iｆ(an！=bm）{ﻩpｒintf("erroｒ!ｃanｎot do thｅｍultiplicａｔiｏn．＼n");ﻩﻩreturn;ﻩ}ｆor(i=0;i<am;i＋＋)ﻩfor(j＝0;j<bn;j＋+）ﻩ{ﻩﻩu=ｉ＊bn+ｊ;C[ｕ］=0．０;ﻩfor(l=0；l<ａn；l+＋)ﻩﻩﻩC[u]+=A[i*ａn+l］*B［l＊bn+j];ﻩ}ﻩrｅturｎ;}douｂle *ｉnv(doｕblｅ*ａ,ｉｎt n) //计算矩阵的逆,本程序的难点｛//采用高斯-约旦－全选主元法iｎt *ｉs,*jｓ,ｉ,ｊ,ｋ,l，u,ｖ;douｂｌeｄ，p;ｉs＝(int*)mａlloc(n*sizｅof（ｉnt));js=(int＊)malloc（n*sｉzｅｏf(int））;ｆoｒ(k=０; ｋ<=ｎ-1; k++){ﻩd=0.０；fｏr(i＝ｋ;ｉ<ｎ;ｉ+＋)fｏｒ(j=ｋ;j<n;j++）{ l＝i*n+j;ﻩp=fabｓ(a［l]);ｉf (p>d)ﻩ{ﻩﻩﻩﻩd=p;ﻩﻩis［ｋ]=i;ﻩﻩjs[k]＝j；ﻩ}}ｉf (d+１.0=＝1.0）{ｆrｅe(is）;fｒｅe(js);printf("error not ｉnｖ\n");return NULL;｝if (ｉs[k]!＝k）for （ｊ＝0；ｊ<ｎ;j+＋){u=k*n+j;ﻩv=is［k]*n+j;p=ａ[u］;ﻩａ[u］＝a[ｖ];ﻩﻩa[v］=p;｝if(js[k]！=k)for （i=0;i<ｎ；i++)｛ｕ=i*n+ｋ;ﻩv=ｉ*n+js[k］;p＝a[u];ﻩa[ｕ]=a[v];ﻩa[v]=p；}l=k*n+ｋ;a[l]=1.0/a［ｌ］;foｒ(j=0；j<n;j++)iｆ(j!=k）｛ﻩu=ｋ＊n+j;ﻩﻩａ［ｕ］=a[u]*a[l]；ﻩ｝foｒ(i=0;i<n;i++)iｆ(i!=k)ｆor (ｊ=0;j＜n;ｊ+＋)ｉf (ｊ!=ｋ){ｕ＝i*n＋j;ａ［u]=a［u]-a[i*n+k］*a[k*n+ｊ];｝for(i＝0;ｉ<n;i＋＋)ｉf (i！=k){ u=ｉ＊n+k；ﻩa［u]=-ａ[u]＊a[l］;ﻩﻩ}｝ﻩfor (k=n－1;k>=0;k--){ if（jｓ[ｋ］!=k）ｆｏr (j=0;j<＝ｎ－1；j+＋){ u=k＊ｎ+j；ﻩﻩv=ｊｓ[ｋ]*n＋j；p=a［u];ﻩａ[u]=ａ［v］;ﻩﻩﻩa［ｖ]=p;}if(is[ｋ]!=ｋ)for (i=0;i＜n;ｉ++){u=i*n+k;ﻩv＝i*n+is[ｋ];p＝ａ［u];ﻩa［u］＝a［ｖ];ﻩﻩa[v］＝ｐ；｝｝ｆrｅe(ｉs）;ﻩfree(js);rｅtｕrn ａ;}voｉd ｐrintmａtrix(dｏuble M[］,ｉnt ｍ,int n) ／/矩阵的输出｛int i,j;for(i=０;i＜ｍ；ｉ+＋){for(j＝0;j＜n;j＋+)printf("%.５f＂，M[i*ｎ+j]）；ﻩpｒinｔf(＂\ｎ"）;}pｒiｎtf("\n"）;｝mａin() /／主函数,空间后方交会的计算{ﻩFＩLE*ｆp；//定义一个文件指针ｆpﻩint m，i，j=０；ﻩdｏubｌe f,t,w,k,Ｓ１=0．0,S2=0.0，S3=０.0,x[Ｎ],y[Ｎ］,x0[N］，y0［Ｎ],X［N］,Y[Ｎ],Z ［N］，Ｘｓ０，Ys0,Ｚs0;ﻩｄｏuｂle a[3],ｂ[3]，c[3],A［2*N＊6],AT[6*2＊Ｎ],ＡTA[6*6]，*ATＡ＿=NULL,l[２*N],ＡTl[6],V［6];if(（fp=fopen（"e：\\ｓｈuju.txt","r"))=＝NUＬＬ）/／使fp指向被打开的shujｕ.tｘt文件ﻩ{ //并判断文件是否打开成功ﻩｐriｎtｆ（＂\nerroｒon ｏpｅn sｈｕju．ｔｘt\n"）;geｔcｈ（）；eｘｉt(１);｝fｏｒ（i=0;i<Ｎ;i＋+）ﻩﻩfscanf(fp,"%lｆ％lf％lf%lf%lf",＆ｘ[i］,＆ｙ［i],&X［ｉ],＆Y[ｉ],&Ｚ[i])；//将文件中的数据赋给主函数定义的变量ﻩpｒｉnｔf("原始数据:＼n");pｒｉnｔf("ｘ\ｔ\ty\t\t\X\ｔ＼tY＼t＼tZ\ｔ\ｔ\n"); /／输出文件中的原始数据ｆor（i＝0；i<N;i++)ﻩﻩｐriｎtｆ("%lf %ｌf %lf %lf ％lf\n"，x[ｉ]，y[i],X[i],Y[i]，Z［i]);pｒintｆ(＂＼n请输入摄影机主距和摄影比例尺分母；f，m：");ｓcanf("％lf,%lf",&f,&ｍ); //输入f，mf=f/10０0.0;fｏr(i=0;i＜N;i++）ﻩ{ﻩx［ｉ]／=100０．0;ﻩﻩy［i]/=１00０.0;S1＋=X[i]；S２＋＝Y［i];ﻩS３+＝Z[i];}Xs０＝S１/N；ﻩYs0=S2/Ｎ; //计算外方位元素的初始值ﻩZs０=m*f+Ｓ３／N;t=0.0;w=０.0;ｋ=０.0；whiｌｅ(j<3)｛pｒｉnｔf(＂--－----－-－-----－－-－-－－-－--------－第%d次计算-----－-－---－－-－－-----－-－－---－-＼n"，j+1);ﻩa[0]=ｃｏs(t）＊cｏs(k）－siｎ（ｔ)＊sin(w)＊sｉｎ（k);a[1]=-cｏs（t)*ｓin(k)-ｓin（t)*sin(w)*cos(ｋ）;ﻩa［2]=-sin（t)*cｏs(ｗ);ｂ［０］=cos(ｗ)＊ｓin(k）；ﻩb[1]=cｏｓ(ｗ)＊ｃos(k); /／计算旋转矩阵b[2］=-sin(w);c[0]＝ｓin(t)*ｃos(k）+cos(t)*ｓin(w）*ｓin（k)；c［1]=-sin(ｔ)*sin(k)+cos（t)*ｓｉｎ(w）*coｓ（ｋ）;c[２]=ｃos(t)*cos(w);fｏr(ｉ=０;i<N；i++)ﻩ{ﻩx0[i]＝－f*（a［０]*（X[i]-Xｓ0)＋b[０］*(Ｙ[i]-Yｓ0)＋ｃ[0]*(Z[ｉ]-Zs０）)／（a[2]*（X[i]-Xs0)+ｂ[2]*(Ｙ［i]-Ｙs0）+c［2]*(Z[i]-Zs0)）; /／计算像点坐标近似值y0[i]＝-f*(a[1］*(X[i]-Ｘs0）+b[１]*(Ｙ[i］-Ys０)+c［1]*(Z［i]-Zs0))/(ａ[２]*(X［i］-Xｓ0）＋ｂ[2]*(Ｙ[i]-Ｙs0）＋c［2］*(Z[i]-Ｚs0));ﻩG[i］=ａ［２]*（X[i]-Xs０)＋b[2］*(Y［ｉ]-Ys０)+ｃ［2］*(Z［ｉ]－Zｓ0)；ﻩ｝for（i=0;i<N;i＋+)｛ﻩA[i*1２+0]=（a[0]*f+a[2］＊x［ｉ])/G［i]；A［ｉ*12＋1]=（b[0]＊ｆ+b[2]＊x［ｉ］)/G［i］;A［i*12+2］＝(c［0]*f＋ｃ[2］*x[i］)/G[i]；ﻩﻩA［i*1２+３]=y[ｉ]*sin(w)-(x[i]＊(x[i]*cｏs(k)-ｙ[i]*sｉn(ｋ）)/f＋f*cos(k）)＊cos(w）; ／/计算误差方程的系数阵以及lﻩA[i*1２+4]=-f＊siｎ(k)-ｘ[i］*（x[i］*sｉn(k)+ｙ[i］*cos(k））/f；ﻩＡ［i*1２+5]=ｙ[i];ﻩA[i*12+6］=(a[1]*ｆ+a[2]*y[i])/G［ｉ]；ﻩA［ｉ＊12+7]=（ｂ[1]*ｆ+b［2]*y[i]）／G［ｉ];ﻩﻩA［i*1２+8］=(c[1］*ｆ＋ｃ[2]＊y[ｉ］）/G[i]；ﻩＡ[i＊１２+９]＝-x[ｉ］*sｉn(ｗ)-(y[i］＊(x[i］＊cos（k）-y[i］*ｓｉｎ(k))/f-ｆ*sin(k））*cos(w);ﻩA[i＊12+10]＝-ｆ＊coｓ（k)-y[ｉ]*(x[i]*siｎ(k)＋ｙ[i］＊cos(k))/f;ﻩA[i*12+11]=-ｘ[i];ﻩﻩｌ[i＊2+0]=x［ｉ］－ｘ0[i];ﻩﻩl［ｉ*2+１]=ｙ[ｉ]－ｙ0[i］;ﻩ}// printf("ｏutput ｍaｔrix: Ａ\n"）；／/ﻩｐrintmatrｉx(Ａ,2*N,6）;//ﻩｐｒintf（"output matrix: ｌ\ｎ")；／/prｉｎtmatrix(l，2*N，1);turn(Ａ,AT,２*Ｎ,6);／/ prinｔf（"ｏutｐut matrix: AＴ\n")；／/ prｉnｔｍａtrｉx（A T,6,2*Ｎ)；mulAB(ＡＴ，A,ATA,６,2＊N,2*N,6); ／/间接平差法计算外方位元素,中间计算的矩阵可以根据需要//ﻩpｒiｎｔf(＂output matrｉx：ATA\n")；//选择性的输出,这里将其屏蔽,为了在打印输出结果的时候／/ﻩprintmatrix(AＴA,6，6);//节约资源ﻩATA_=ｉnv(ＡＴA,6);//ﻩprｉntf（"ｏuｔput mａtｒiｘ:A TA＿\n＂);//ﻩpｒintmatｒｉx（AＴＡ_,6,６)；mｕlAB(AT,l,ATl,６,２*N，2*Ｎ，1)；//ﻩprｉntf("ｏutpit matrinｘ: ＡTl\n");// printmａtriｘ(ATl，６,１)；mulAB（ATＡ_,A Tl,Ｖ,6，6,６,１)；// pｒinｔf（"ouｔpuｔmatriｘ: V\n");／/ﻩpriｎtｍatrix(Ｖ,６,1）；ﻩXs0+=V［0];Yｓ0+=Ｖ[1］;ﻩZs０+=V[2];ﻩt＋=Ｖ[3];ﻩw+=V[4］;k＋=V[5］;prｉｎtf（"第%ｄ次计算的外方位元素为：\n",++j）；ｐｒintf("Ｘｓ＝%.5lｆ,Ys=%.5lｆ，Zs=%.５lf，ｔ＝%.５lf,w=％．5ｌf,k=%.５ｌf\ｎ",Xs０,Ys0,Zs0,t,ｗ,k);}prｉntｆ(＂\n外方位元素为:\n");ｐrintｆ（"Xｓ＝％.5lf,Yｓ=％.5ｌf,Ｚs＝％．5ｌf,t=%.5lf,ｗ=％.5lf,k＝%．5lｆ\n",Xｓ0,Yｓ0,Zs０，ｔ,ｗ,ｋ)；ｆclｏｓe(ｆｐ);｝程序运行的结果为：五.实验总结：通过这次的实验我学到了很多的东西,通过编程加深了对摄影测量空间后方交会相关知识的理解。