摄影测量学单像空间后方交会程序设计作业
摄影测量后方交会
单张相片后方交会目录●作业任务 (3)●解算原理 (3)●具体过程 (4)●算法描述及程序流程 (4)●计算结果 (7)●结果分析 (8)●心得体会及建议 (8)●参考文献 (9)一，作业任务已知摄影机主距f=153.24mm，四对点的像点坐标与相应地面坐标列入下表：表1-1计算近似垂直摄影情况下后方交会解。
二，解算原理【关键词1】中心投影构像方程在摄影测量学中，最重要的方程就是中心投影构像方程（图2-1）。
这个方程将地面点在地面摄影测量坐标系中的坐标（物方坐标）和地面点对应像点的像平面坐标联系起来。
在解析摄影测量与数字摄影测量中是极其有用的。
在以后将要学习到的双像摄影测量光束法、解析测图仪原理及数字影像纠正等都要用到该式。
图2-1在上述公式中：x和y分别为以像主点为原点的像点坐标，相应地面点坐标为X,Y,Z，相片主距f以及外方位元素Xs，Ys，Zs，ψ，ω，κ。
而在此次作业中，就是已知四个地面控制点的坐标以及其对应的像点坐标，通过间接平差原理来求解此张航片的外方位元素。
【关键词2】间接平差在一个平差问题中，当所选的独立参数X的个数等于必要观测值t时，可将每个观测值表达成这t个参数的函数，组成观测方程，然后依据最小二乘原理求解，这种以观测方程为函数模型的平差方法，就是间接平差方法间接平差的函数模型为：随机模型为：平差准则为：VtPV=min【关键词3】单像空间后方交会利用至少三个已知地面控制点的坐标A（Xa，Ya，Za）、B（Xb，Yb，Zb）、Z（Xc，Yc，Zc），与其影像上对应的三个像点的影像坐标a（xa，ya）、b（xb，yb）、c（xc，yc），根据共线方程，反求该像点的外方位元素Xs，Ys，Zs，ψ，ω，κ。
这种解算方法是以单张像片为基础，亦称单像空间后方交会。
在此次作业中，就是已知四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标。
由此可以列出8个误差方程，存在两个多余观测数，则n=2。
摄影测量学单像空间后方交会编程实习报告(精品资料).doc
【最新整理，下载后即可编辑】摄影测量学单像空间后方交会编程实习报告班级：130x姓名：xx学号：2013302590xxx指导老师：李欣一、实习目的通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。
了解摄影测量平差的基本过程，掌握空间后方交会的定义和实现算法。
二、实习内容根据学习的单像空间后方交会的知识，用程序设计语言（C++或C语言）编写一个完整的单像空间后方交会程序，通过对提供的数据进行计算，输出像片的外方位元素并评定精度。
三、实习数据已知航摄仪的内方位元素：fk =153.24mm，x=y=0，摄影比例尺为1:15000；4个地面控制点的地面坐标及其对应像点的像片坐标：四、实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。
因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。
可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。
单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。
五、 实习流程1. 获取已知数据。
从摄影资料中查取影像比例尺1/m ，平均摄影距离（航空摄影的航高、内方位元素x 0，y 0，f ；获取控制点的空间坐标X t ，Y t ，Z t 。
2. 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
单片空间后方交会程序设计
单片空间后方交会程序设计
1 目的
用程序设计语言（VC或者VB）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。
本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
2. 内容
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。
3. 数据准备
已知航摄仪的内方位元素：f
k =153.24mm，x
=y
=0.0mm，摄影比例尺为1:50000；
4个地面控制点的地面坐标及其对应像点的像片坐标：
4. 操作步骤
上机调试程序并打印结果。
摄影测量-空间前交、后交【精选文档】
空间后交—前交程序设计(实验报告)姓名：班级：学号：时间:空间后交-前交程序设计一、实验目的用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序⑴提交实习报告：程序框图、程序源代码、计算结果、体会⑵计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度二、实验数据f=150。
000mm，x0=0，y0=0三、实验思路1。
利用空间后方交会求左右像片的外方位元素(1）.获取m（于像片中选取两点，于地面摄影测量坐标系中选取同点,分别计算距离，距离比值即为m)，x，y，f,X，Y,Z（2).确定未知数初始值Xs,Ys，Zs,q，w,k（3）.计算旋转矩阵R(4）.逐点计算像点坐标的近似值(x），（y）(5）。
组成误差方程式(6)。
组成法方程式（7).解求外方位元素（8）。
检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤(3）-（7）2。
利用求出的外方位元素进行空间前交，求出待定点地面坐标(1）.用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1，R2(2）。
根据左、右像片的外方位元素，计算摄影基线分量Bx，By，Bz(3)。
计算像点的像空间辅助坐标（X1,Y1，Z1）和（X2，Y2,Z2）（4）.计算点投影系数N1和N2(5）。
计算未知点的地面摄影测量坐标四、实验过程⑴程序框图函数AandL%求间接平差时需要的系数％%％已知％a=像点坐标x，b=像点坐标y,f内方位元素主距%φ=q，ψ=w，κ=k%像空间坐标系X,Y，Z%地面摄影测量坐标系Xs,Ys,Zsfunction ［A1，L1，A2，L2］=AandL（a，b，f,q,w,k,X，Y，Z,Xs,Ys,Zs) %％％％％%％％%%％选择矩阵元素a1=cos（q)＊cos(k)—sin(q）＊sin(w）*sin(k）;a2=-cos(q)＊sin(k)—sin(q）＊sin(w）＊cos（k)；a3=-sin(q）*cos(w);b1=cos(w）*sin(k)；b2=cos（w）*cos（k)；b3=—sin(w);c1=sin（q)*cos(k)+cos（q)＊sin(w）*sin（k)；c2=—sin（q）＊sin(k)+cos（q)*sin(w）*cos（k)；c3=cos(q)＊cos(w);％%％%％％％共线方程的分子分母X_=a1*(X—Xs）+b1*（Y-Ys)+c1＊（Z-Zs);Y_=a2*(X-Xs)+b2＊（Y—Ys)+c2＊（Z-Zs);Z_=a3＊(X—Xs)+b3＊（Y—Ys）+c3*(Z-Zs）;％％％％%％%近似值x=-f＊X_/Z_；y=-f＊Y_/Z_;％％％％%％％A组成L组成a11=1/Z_＊（a1＊f+a3＊x)；a12=1/Z_＊(b1＊f+b3＊x）；a13=1/Z_*(c1*f+c3＊x);a21=1/Z_*(a2＊f+a3＊y）;a22=1/Z_*（b2＊f+b3*y);a23=1/Z_＊(c2＊f+c3＊y）;a14=y*sin（w)-（x/f＊(x*cos(k)—y*sin（k））+f＊cos（k）)*cos(w）;a15=-f*sin(k)—x/f＊（x*sin（k）+y*cos(k));a16=y；a24=—x*sin(w）-（y/f＊（x*cos(k)-y*sin(k)）—f＊sin(k））*cos（w）；a25=-f*cos（k）-y/f*(x*sin（k）+y*cos(k));a26=-x;lx=a—x；ly=b-y;％%%%%%％％%组成一个矩阵，并返回A1=[a11,a12，a13,a14,a15，a16];A2=［a21,a22,a23，a24,a25，a26];L1=lx;L2=ly;函数deg2dms％%％%％%％％角度转度分秒function y=deg2dms(x）a=floor(x)；b=floor(（x-a)＊60);c=(x-a—b/60）＊3600；y=a+（b/100）+(c/10000);函数dms2deg%%％％%度分秒转度function y=dms2deg（x）a=floor(x）；b=floor((x-a)＊100）；c=（x-a—b/100)*10000；y=a+b/60+c/3600;函数ok%％％％％%%％%%%％％％目的是为了保证各取的值的有效值%％xy为n*1，a为1＊nfunction result=ok(xy，a)format short gi=size（xy，1)；for n=1:io=xy（n)—floor(xy(n,1));o=round（o＊（10^a（n))）/(10^a(n）);xy(n，1)=floor(xy（n,1）)+o;endformat long gresult=xy;函数rad2dmsxy%%％%求度分秒表现形式的三个外方位元素，三个角度function xydms=rad2dmsxy(xy)［a,b,c，d，e,f]=testvar(xy）;d=deg2dms（rad2deg(d））；e=deg2dms（rad2deg(e))；f=deg2dms(rad2deg(f）);xydms=[a,b，c,d,e，f］';函数spacehoujiao％％%％%%%空间后交％%% f%％输入p（2*n,1）%％像点坐标x,y,X,Y，Z,均为（n,1）function ［xy,m，R］=spacehoujiao(p，x,y，f,X，Y,Z）format long；%%％%％权的矢量化,这是等精度时的，如果非,将函数参数改为PP=diag（p)；％％求nj=size（X，2）；%％初始化Xs=0；Ys=0;Zs=0;for n=1：jXs=Xs+X（n）；Ys=Ys+Y(n）;Zs=Zs+Z(n）;endSx=sqrt((x（2）-x(1）)^2+(y（2)—y(1)）^2);％%％%两像点之间距离Sd=sqrt(（X(2）-X（1))^2+(Y（2）-Y（1))^2);%％%％两地面控制点之间距离m=Sd/Sx; ％％％％图像比例系数Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j；m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6）；l=rand（2＊j，1);％％％％for n=1:j[a（2＊n—1,:)，l(2*n—1,1）,a(2＊n,：），l(2*n,1）］=AandL(x(n)，y（n),f,q，w，k，X(n)，Y（n），Z(n），Xs,Ys,Zs);enddet=inv（a’＊P*a）*transpose（a)*P*l;％%%%%%％％％循环体while 1％％％％％%%％%％％％％%％％[dXs,dYs,dZs，dq,dw,dk]=testvar(det）;detXs=abs（dXs）；detYs=abs（dYs)；detZs=abs(dZs)；detq=abs(dq)；detw=abs（dw)；detk=abs（dk)；%%％％%%％%%if （（detXs<0。
摄影测量学空间后方交会实验报告
摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘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．实验结果四．实验总结此次实验让我深入了解单像空间后方交会的计算过程，加强了对空间后方交会基本公式和误差方程式，法线方程式的记忆。
摄影测量实验报告(空间后方交会—前方交会)
空间后方交会—空间前方交会程序编程实验一．实验目的要求掌握运用空间后方交会-空间前方交会求解地面点的空间位置.学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。
然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程.二．仪器用具计算机、编程软件(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;}}}。
单像空间后方交会实验报告(c++版)
单像空间后方交会姓名：学号：时间：Echo did this for you .2013/4/25目录一、作业任务 .............................................................................................................. - 4 -二、计算原理 .............................................................................................................. - 4 -三、算法流程 .............................................................................................................. - 8 -四、源程序 .................................................................................................................. - 9 -五、计算结果 .............................................................................................................. - 9 -六、结果分析 .............................................................................................................. - 9 -七、心得与体会 .......................................................................................................... - 9 -八、附页 ...................................................................................................................... - 9 -1.c++程序........................................................................................................... - 9 -2.C++程序截图..................................................................................................- 16 -3.matlb程序.....................................................................................................- 17 -一、 作业任务 已知条件：摄影机主距f=153.24mm ，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。
五上、数字摄影测量学单片空间后方交会
总误差方程
法方程
V Ax L
x (AT A) 1 (AT L)
X s Ys V1 A1 l1 Z V2 A2 l2 s V , A , L , x , Vn An ln T T li xi ( xi ) yi ( yi ) , Vi v xi v yi a11 a12 a13 a14 a15 a16 Ai a21 a22 a23 a24 a25 a26
已知点必须多余点， 数据处理方法采用 最小二乘法！
这是所有测量的一个统一的基本原则！ 摄影测量也不例外。
二、误差方程与法方程
已知值 x0 , y0 , f ，m, X, Y, Z 观测值 x , y 相应改正数 vx,vy 未知数 Xs, Ys, Zs, , , 泰勒级数展开
四、空间后方交会的精度
求解各未知数的精度可以通过法方程系数矩阵 求逆的方法，解出相应的权倒数 Qii
mi m0 Qii 按下式计算第i未知数的中误差：
式中，m0为单位权中误差，计算公式 为： m [VV ] 0 2n 6 ,其中n为控制点的点数。
空间后方交会用到的已知点越多，空间后方交会 的精度越高，此外空点的分布也空间后方交会计算 的精度。空间后方交会使用的控制点应当避免位于 一个圆柱面上，否则，会出现解不唯一的情况。
偏导数 1
x f X Z 2 ( Z X) X s Z X s X s f 2 ( a1Z a3 X ) Z 1 X (a1 f f a3 ) Z Z 1 (a1 f a3 x) Z
偏导数 2
x f X Z 2 ( Z X) Z
摄影测量作业3-空间后方交会计算
CFileDialog dlgOpenFile(TRUE, _T("txt"), NULL, OFN_FILEMUSTEXIST, _T("(文本文件)|*.txt|(所有文件）|*.*)||"));
if (dlgOpenFile.DoModal() == IDCANCEL) return;//如果选择取消按钮，则退出
原理、算法流程、源程序、计算结果、结果分析、心得体会等。
三．实验所用到的数学公式及程序计算步骤。
单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标 根据共线方程反 求影像的外方位元素。 数学模型：共线条件方程式：
3
求解过程： （1）获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测
CMatrix X,_A,_AA,N_AA; _A = ~A;//A 的转置 _AA = _A*A; N_AA = _AA.Inv();//_AA 的逆矩阵 X = N_AA*_A*L; return X; }
CMatrix CKongJianHouFangJiaoHuiDlg::GetA(CMatrix xyXYZ, double f, CMatrix XX)//计算系数矩 阵A {
CMatrix CKongJianHouFangJiaoHuiDlg::GetL(CMatrix xyXYZ, double f, CMatrix XX)//计算 L 矩阵 {
int iRow = xyXYZ.Row(); CMatrix L(2 * iRow, 1); double XS = XX(0, 0); double YS = XX(0, 1); double ZS = XX(0, 2);
A(2*i, 3) = y*sin(w) - (x*(x*cos(k) - y*sin(k)) / f + f*cos(k))*cos(w); A(2*i, 4) = -f*sin(k) - x*(x*sin(k) + y*cos(k)) / f; A(2*i, 5) = y; A(2*i+1, 0) = (a2*f + a3*y) / _Z; A(2 * i + 1, 1) = (b2*f + b3*y) / _Z; A(2 * i + 1, 2) = (c2*f + c3*y) / _Z; A(2 * i + 1, 3) = -x*sin(w) - (y*(x*cos(k) - y*sin(k)) / f - f*sin(k))*cos(w); A(2 * i + 1, 4) = -f*cos(k) - y/ f*(x*sin(k) + y*cos(k)); A(2 * i + 1, 5) = -x; } return A; }
简述单像空间后方交会的程序设计步骤
简述单像空间后方交会的程序设计步骤
单像空间后方交会是一种用于测量摄影点在三维空间中位置的方法。
以下是简述的程序设计步骤：
1.读取摄影测量数据：首先，从摄影测量设备（如相机）中读取图像和相关的内参数据，包括相机的焦距、像点大小等。
2.图像处理：对读取的图像进行预处理。
可能需要进行去畸变操作，校正图像的畸变。
3.特征提取：从图像中提取关键点或特征点。
这些特征点可以是角点、边缘、斑点等。
提取出的特征点用于后方交会计算。
4.求解相机位姿：使用特征点的像素坐标和已知内参数，通过解非线性方程组的方法，计算相机在三维空间中的位姿（即相机的位置和方向）。
5.求解三维点坐标：对于每个特征点，使用单像模型，将像素坐标投影到相机坐标系中。
然后，通过解线性方程组的方法，计算特征点在三维空间中的坐标。
6.误差检测与优化：计算测量误差，并进行误差检测。
可以使用一些优化算法，如最小二乘法，来优化相机位姿和三维点坐标。
7.输出测量结果：将计算得到的三维点坐标输出，可以是数字格式或者可视化结果。
以上是单像空间后方交会的基本程序设计步骤。
每个步骤可能会有不同的具体实现，根据具体的应用场景和需求进行设计和调整。
摄影测量学单像空间后方交会编程实习报告
摄影测量学单像空间后方交会编程实习报告本次实习中，我使用编程语言进行了单像空间后方交会的实现，并取得了一定的成果。
首先，我了解了单像空间后方交会的基本原理。
根据光线在透镜上的成像规律，可以推导出物体点在像平面上的坐标与图像点在像平面上的坐标之间的关系式。
通过已知的摄像机内外方位元素和图像点坐标，可以反求得物体点的坐标。
在程序编写过程中，我采用了Python编程语言。
首先，我定义了一个类，用于存储摄像机的内外方位元素和图像点坐标。
然后，我编写了一个函数，用于计算物体点的坐标。
该函数根据已知的内外方位元素和图像点坐标，使用逆向投影的方式反求物体点的坐标。
最后，我编写了一个主函数，通过读取输入文件中的数据，调用计算函数，并将结果保存到输出文件中。
在实现的过程中，我遇到了一些问题。
首先，由于摄像机的内外方位元素需要提前获取，因此我通过测量方法获得了实际的内外方位元素。
然而，测量的过程中存在一定的误差，因此在计算物体点坐标时可能存在一定的误差。
其次，图像坐标与物体点坐标之间的关系式中存在一些参数，如焦距、主点坐标等，这些参数也需要提前获取。
在程序中，我将这些参数作为输入参数，通过外部文件进行输入。
在实习的过程中，我充分运用了自己所学的摄影测量学知识，并将其与编程技能相结合。
在实现过程中，我遇到了一些难题，但通过查阅资料和与老师的讨论，最终得以解决。
通过编程实习，我深入理解了单像空间后方交会的原理，并通过实际操作提高了自己的计算能力。
总的来说，本次实习使我对摄影测量学有了更深入的认识，也提升了我的计算和编程能力。
通过此次实习，我对摄影测量学的兴趣更加浓厚，也更加期待在今后的学习和研究中能够进一步深入探索。
(完整word版)单像空间后方交会程序报告
单像空间后方交会程序报告指导老师：刘老师班级：测绘 101姓名：尚锋学号：19号1、应用程序的主进口部分的代码：using System;using System.Collections.Generic;using System.Linq;using System.Windows.Forms;namespace 单像空间后方交会{static class Program{///<summary>///应用程序的主进口点。
///</summary>[ STAThread]static void Main(){Application .EnableVisualStyles();Application .SetCompatibleTextRenderingDefault( false );Application .Run( new Form1());}}}2、方法解算类（通用）部分的代码：using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace 单像空间后方交会{class Tongyong{struct image_point // 一个像点结构，包括像点坐标和地面点坐标{public double x;public double y;public double X;public double Y;public double Z; }private private private private private private private private private private private private private publicdouble f;// 主距double u;//u 为外方向元素，下边 5个同样double w;double k;double Xs;double Ys;double Zs;image_point [] p = new image_point [4];// 四个控制点double [] R = new double [9]; // 旋转矩阵double [] a = new double [8];// 像点坐标近似值double [,] A =new double [8, 6];// 偏差方程式系数double [] L = new double [8];// 偏差方程式常数项int count = 0;// 统计代次数Tongyong( double g, double [] q)// 结构函数，初始化各变量 , 单位 m{f = g;for ( int i = 0; i < 4; i++){int j = i * 5;p[i].x = q[j];p[i].y = q[j + 1];p[i].X = q[j + 2];p[i].Y = q[j + 3];p[i].Z = q[j + 4];}double ave = 0, sum = 0;for ( int i = 0; i < 3; i++)// 求比率尺分母{for ( int j = i + 1; j < 4; j++){sum += Math.Pow(p[i].Y - p[j].Y, 2)) / + Math.Pow(p[i].y - p[j].y, 2));} Math.Sqrt(Math.Pow(p[i].X - p[j].X, 2) + Math.Sqrt( Math.Pow(p[i].x - p[j].x, 2)}ave = sum / 6;u = 0;// 给定外方向元素的初始值w = 0;k = 0;Xs = (p[0].X + p[1].X + p[2].X + p[3].X) / 4; , 角度均设置为0//Xs 为四个控制点 X的均匀值，Ys近似Ys = (p[0].Y + p[1].Y + p[2].Y + p[3].Y) / 4;Zs = (p[0].Z + p[1].Z + p[2].Z + p[3].Z) / 4 + ave * f;}private double sin( double m) // 正弦，为简化而写 , 下同{return Math.Sin(m);}private double cos( double m){return Math.Cos(m);}private void calcos()// 计算旋转矩阵{R[0] = cos(u) * cos(k) - sin(u) * sin(w) * sin(k);R[1] = -cos(u) * sin(k) - sin(u) * sin(w) * cos(k);R[2] = -sin(u) * cos(w);R[3] = cos(w) * sin(k);R[4] = cos(w) * cos(k);R[5] = -sin(w);R[6] = sin(u) * cos(k) + cos(u) * sin(w) * sin(k);R[7] = cos(u) * sin(w) * cos(k) - sin(u) * sin(k);R[8] = cos(u) * cos(w);}private void calabout() // 像点坐标的近似值{int i;for (i = 0; i < 4; i++){a[2 * i] = -f * (R[0] * (p[i].X - Xs) + R[3] * (p[i].Y-Ys) + R[6] * (p[i].Z - Zs)) / (R[2] * (p[i].X - Xs) + R[5] * (p[i].Y-Ys) + R[8] * (p[i].Z - Zs));a[2* i + 1] = -f * (R[1] * (p[i].X - Xs) + R[4] * (p[i].Y -Ys) + R[7] * (p[i].Z - Zs)) / (R[2] * (p[i].X - Xs) + R[5] * (p[i].Y-Ys) + R[8] * (p[i].Z - Zs));}}private void calxx() // 偏差方程式的系数和常数项{int for i;(i = 0; i < 4; i++) // 系数{double z = R[2] * (p[i].X - Xs) + R[5] * (p[i].Y - Ys) + R[8] * (p[i].Z - Zs);int n = i * 2;A[n, 0] = (R[0] * f + R[2] * p[i].x) / z;A[n, 1] = (R[3] * f + R[5] * p[i].x) / z;A[n, 2] = (R[6] * f + R[8] * p[i].x) / z;A[n, 3] = p[i].y * sin(w) - f * cos(w) * cos(k) - p[i].x/f * (p[i].x * cos(w) * cos(k) - p[i].y * cos(w) * sin(k)); A[n,4] = -f * sin(k) - p[i].x / f * (p[i].x * sin(k) +p[i].y * cos(k));A[n, 5] = p[i].y;A[n + 1, 0] = (R[1] * f + R[2] * p[i].y) / z;A[n + 1, 1] = (R[4] * f + R[5] * p[i].y) / z;A[n + 1, 2] = (R[7] * f + R[8] * p[i].y) / z;A[n + 1, 3] = -p[i].x * sin(w) + f * cos(w) * sin(k) -p[i].x / f * (p[i].x * cos(w) * cos(k) - p[i].y * sin(k) * cos(w));A[n + 1, 4] = -f * cos(k) - p[i].y / f * (p[i].x * sin(k) + p[i].y * cos(k));A[n + 1, 5] = -p[i].x;}for (i = 0; i < 4; i++)// 常数项{L[2 * i] = p[i].x - a[2 * i];L[2 * i + 1] = p[i].y - a[2 * i + 1];}}private double calAdd(){double [,] temp =new double [6, 6];//A 的转置与 A相乘的积double [,] ANew = new double [6, 8];//A 的转置double [] t =new double [6];//A 的转置与 L相乘的积double [] X = new double [6];// 更正数int i, j, n;for (i = 0; i < 8; i++)// 求A的转置 ANew{for (j = 0; j < 6; j++){ANew[j, i] = A[i, j];}}for (i = 0; i < 6; i++)// 求A的转置与 A相乘的积 temp {for (j = 0; j < 6; j++){temp[i, j] = 0;for (n = 0; n < 8; n++){temp[i, j] += ANew[i, n] * A[n, j];}}}MATINV(temp);//temp for (i = 0; i < 6; i++) 的逆，保留在自己矩阵中// 求A的转置与 L的乘积 t{t[i] = 0;for (j = 0; j < 8; j++){t[i] += ANew[i, j] * L[j];}}for (i = 0; i < 6; i++)// 求更正数 X{X[i] = 0;for (j = 0; j < 6; j++){X[i] += temp[i, j] * t[j];}}Xs += X[0];// 外方向元素初始值加上更正数Ys += X[1];Zs += X[2];u += X[3];w += X[4];k += X[5];return maxone(X);// 返回判断条件 , 最大的更正数的值}public void makeSure() // 计算流程控制函数{calcos();calabout();calxx();double VALUE = calAdd();count++;while (VALUE > 0.00001)// 迭代至最大更正数为止{calcos();calabout();calxx();VALUE = calAdd();count++;}}private void MATINV(double [,] c) // 求6阶矩阵的逆{int i, j, h, m;const int n = 6;double l;double [,] q = new double [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 = 0, m = 0; m < n - 1; m++, h++)// 消去对角线以下的数据{for (i = m + 1; i < n; i++){if(q[i, h] == 0d)continue ;l = q[m, h] / q[i, h];for (j = 0; j < 12; j++){q[i, j] *= l;q[i, j] -= q[m, j];}}}for (h = n - 1, m = n - 1; m > 0; m--, h--)//消去对角线以上的数据{for (i = m - 1; i >= 0; i--){if(q[i, h] == 0d)continue ;l = q[m, h] / q[i, h];for (j = 0; j < 12; j++){q[i, j] *= l;q[i, j] -= q[m, j];}}}for (i = 0; i < n; i++)// 将对角线上数据化为 1{l = 1.0 / q[i, i];for (j = 0; j < 12; j++){q[i, j] *= l;}}for (i = 0; i < n; i++)// 提取逆矩阵{for (j = 0; j < n; j++){c[i, j] = q[i, j + 6];}}}private double maxone(double [] Arr)// 返回六个元素中的最大值{double [] ARR =new double [6];for ( int i = 0; i < 6; i++){ARR[i] = Arr[i];if (ARR[i] < 0d)// 取正ARR[i] = -ARR[i];}Array .Sort(ARR, 0, 6);return ARR[5];}public int COUNT{get{return count;}}public double U {get{return u;}}public double W {get{return w;}}public double K {get{return k;}}public double XS {get{return Xs;}}public double YS {get{return Ys;}}public double ZS {get{return Zs;}}}3、窗体一部分的代码：using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;using System.Diagnostics;using System.Data.OleDb;using System.IO;namespace 单像空间后方交会{public partial class Form1 : Form {private private double [] data= double f = 0;new double [20];// 主距// 保留表中的数据public Form1(){InitializeComponent();}private void Form1_Load( object sender, EventArgs e) {//TODO: 这行代码将数据加载到表“ database1DataSet1.data ”中。
摄影测量程序汇总(后方交会+前方交会+单模型光束法平差)
程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug 下。
1.单像空间后方交会C语言程序：#include <stdio.h>#include <math.h>#include <iostream>double *readdata();void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa);void transpose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,double * mat2,double * result,int a,int b,int c);void inverse(double *a,int n)/*正定矩阵求逆*/{int i,j,k;for(k=0;k<n;k++){for(i=0;i<n;i++){if(i!=k)*(a+i*n+k)=-*(a+i*n+k)/(*(a+k*n+k));}*(a+k*n+k)=1/(*(a+k*n+k));for(i=0;i<n;i++){if(i!=k){for(j=0;j<n;j++){if(j!=k)*(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k);}}}for(j=0;j<n;j++){if(j!=k)*(a+k*n+j)*=*(a+k*n+k);}}}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 multi(double *mat1,double *mat2,double * result,int a,int b,int c) { int i,j,k;for(i=0;i<a;i++){for(j=0;j<c;j++){result[i*c+j]=0;for(k=0;k<b;k++)result[i*c+j]+=mat1[i*b+k]*mat2[k*c+j];}}return;}double *readdata(){FILE *fp;int i,j;int number;char datacatolog[100];//scanf("%s",datacatolog);if ((fp=fopen("控制点坐标.txt","r"))==NULL){printf("读取数据出错!\n");return false;}fscanf(fp,"%d",&number);double *cordata=new double[number*5];for (i=0;i<number;i++){for (j=0;j<5;j++){fscanf(fp,"%lf",cordata+i*5+j);}}printf("控制点坐标数据读取成功！\n");return cordata;}void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa){FILE *fp;char *file1="结算数据.txt";fp=fopen(file1,"w");fprintf(fp,"---------原始坐标数据为--------：\n");for (int i=0;i<hang;i++){for (int j=0;j<5;j++){fprintf(fp,"%7.4lf ",data[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------误差方程系数阵为：--------：\n");for (int i=0;i<hang*2;i++){for (int j=0;j<6;j++){fprintf(fp,"%7.4lf ",xishuarray[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------法方程系数阵为：--------：\n");for (int i=0;i<6;i++){for (int j=0;j<6;j++){fprintf(fp,"%7.5lf ",faxishu[i*5+j]);}fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------误差方程常数项为:--------：\n");for (int i=0;i<hang*2;i++){fprintf(fp,"%lf ",l[i]);fprintf(fp,"\n");}fprintf(fp,"--------------------------------\n");fprintf(fp,"---------迭代次数为:--------：\n");fprintf(fp,"%d\n",i);fprintf(fp,"--------------------------------\n");fprintf(fp,"-----------外方位元素为：---------\n");fprintf(fp," Xs= %lf, Ys=%lf, Zs=%lf\n",xs,ys,zs);fprintf(fp," fai= %lf, oumiga=%lf, kapa=%lf\n",fai,oumiga,kapa);fprintf(fp,"--------------------------------\n");fclose(fp);return;}void main(){int i,j;int ii,jj;int diedainumber=0;double x0=0.0,y0=0.0,f=0.0;double m=50000; //估算比例尺double fai=0,oumiga=0,kapa=0,Xs=0,Ys=0,Zs=0;double R[3][3]={0.0};double X=0.0,Y=0.0,Z=0.0,L[8][1]={0.0},A[8][6]={0.0};double correct[6][1]={0.0},AT[6][8]={0.0},ATA[6][6]={0.0},ATL[6][1]={0.0};int row; //row用于存放坐标行数double *controlpoint;controlpoint=readdata();row=sizeof(controlpoint);for (i=0;i<row;i++){for (j=0;j<5;j++){printf("%3.3lf ",controlpoint[i*5+j]);}printf("\n");}/*----------内方位元素----------*/printf("请输入像片的内方位元素(mm)：\n");printf("x0=");x0/=1000.0;scanf("%lf",&x0); //double类型数据要用%lfprintf("y0=");y0/=1000.0;scanf("%lf",&y0);printf("f=");scanf("%lf",&f);f=f/1000.0;/*------------------------------*///-------确定未知数初始值------for(int i=0;i<row;i++){Xs=Xs+controlpoint[i*5+2];Ys=Ys+controlpoint[i*5+3];Zs=Zs+controlpoint[i*5+4];}Xs/=row;Ys/=row;Zs=Zs/row+m*f;//-----------------------------do{diedainumber++;//---------组成旋转矩阵--------R[0][0]=cos(fai)*cos(kapa)-sin(fai)*sin(oumiga)*sin(kapa);R[0][1]=-cos(fai)*sin(kapa)-sin(fai)*sin(oumiga)*cos(kapa);R[0][2]=-sin(fai)*cos(oumiga);R[1][0]=cos(oumiga)*sin(kapa);R[1][1]=cos(oumiga)*cos(kapa);R[1][2]=-sin(oumiga);R[2][0]=sin(fai)*cos(kapa)+cos(fai)*sin(oumiga)*sin(kapa);R[2][1]=-sin(fai)*sin(kapa)+cos(fai)*sin(oumiga)*cos(kapa);R[2][2]=cos(fai)*cos(oumiga);//-----------------------------//计算系数阵和常数项for(int i=0,k=0,j=0;i<=3;i++,k++,j++){X=R[0][0]*(controlpoint[i*5+2]-Xs)+R[1][0]*(controlpoint[i*5+3]-Ys)+R[2][0]*(co ntrolpoint[i*5+4]-Zs);Y=R[0][1]*(controlpoint[i*5+2]-Xs)+R[1][1]*(controlpoint[i*5+3]-Ys)+R[2][1]*(con trolpoint[i*5+4]-Zs);Z=R[0][2]*(controlpoint[i*5+2]-Xs)+R[1][2]*(controlpoint[i*5+3]-Ys)+R[2][2]*(controlpoint[i*5+4]-Zs);L[j][0]=controlpoint[i*5+0]-(x0-f*X/Z);L[j+1][0]=controlpoint[i*5+1]-(y0-f*Y/Z);j++;A[k][0]=(R[0][0]*f+R[0][2]*(controlpoint[i*5+0]-x0))/Z;A[k][1]=(R[1][0]*f+R[1][2]*(controlpoint[i*5+0]-x0))/Z;A[k][2]=(R[2][0]*f+R[2][2]*(controlpoint[i*5+0]-x0))/Z;A[k][3]=(controlpoint[i*5+1]-y0)*sin(oumiga)-((controlpoint[i*5+0]-x0)*((control point[i*5+0]-x0)*cos(kapa)-(controlpoint[i*5+1]-y0)*sin(kapa))/f+f*cos(kapa))*cos(o umiga);A[k][4]=-f*sin(kapa)-(controlpoint[i*5+0]-x0)*((controlpoint[i*5+0]-x0)*sin(kapa) +(controlpoint[i*5+1]-y0)*cos(kapa))/f;A[k][5]=controlpoint[i*5+1]-y0;A[k+1][0]=(R[0][1]*f+R[0][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][1]=(R[1][1]*f+R[1][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][2]=(R[2][1]*f+R[2][2]*(controlpoint[i*5+1]-y0))/Z;A[k+1][3]=-(controlpoint[i*5+0]-x0)*sin(oumiga)-((controlpoint[i*5+1]-y0)*((cont rolpoint[i*5+0]-x0)*cos(kapa)-(controlpoint[i*5+1]-y0)*sin(kapa))/f-f*sin(kapa))*cos( oumiga);A[k+1][4]=-f*cos(kapa)-(controlpoint[i*5+1]-y0)*((controlpoint[i*5+0]-x0)*sin(ka pa)+(controlpoint[i*5+1]-y0)*cos(kapa))/f;A[k+1][5]=-(controlpoint[i*5+0]-x0);k++;}transpose(A[0],AT[0],8,6);multi(AT[0],A[0],ATA[0],6,8,6);inverse(ATA[0],6);multi(AT[0],L[0],ATL[0],6,8,1);multi(ATA[0],ATL[0],correct[0],6,6,1);Xs=Xs+correct[0][0];Ys=Ys+correct[1][0];Zs=Zs+correct[2][0];fai=fai+correct[3][0];oumiga=oumiga+correct[4][0];kapa=kapa+correct[5][0];}while(correct[3][0]>=6.0/206265.0||correct[4][0]>=6.0/206265.0||correct[5][0] >=6.0/206265.0);printf("迭代次数为：%d\n",diedainumber);printf("---------误差方程系数为：--------\n");for (i=0;i<8;i++){for (j=0;j<6;j++){printf("%4.4lf ",A[i][j]);}printf("\n");}printf("--------------------------------\n");printf("求解得到的外方位元素为：\n");printf(" Xs= %lf\n",Xs);printf(" Ys= %lf\n",Ys);printf(" Zs= %lf\n",Zs);printf(" fai= %lf\n",fai);printf(" oumiga= %lf\n",oumiga);printf(" kapa= %lf\n",kapa);savedata(row,controlpoint,A[0],ATA[0],L[0],diedainumber,Xs,Ys,Zs,fai,oumiga,kap a);printf("-----------------解算结束！--------------\n");system("pause");}解算结果：已知左右像片外方位元素，给出像点坐标：C语言代码：#include <stdio.h>#include <iostream>#include <math.h>double *readdata();void savedata(int hang,double *data);double *readdata(){FILE *fp;int i,j,k;int number;char datacatolog[100];char leftdata[300];//scanf("%s",datacatolog);if ((fp=fopen("像点坐标数据.txt","r"))==NULL){printf("读取数据出错!\n");system("pause");exit(0);}fscanf(fp,"%d",&number);double *c=new double[number*4];for (k=0;k<2;k++){fread(&leftdata,14,1,fp);for (i=0;i<number;i++){for (j=0;j<2;j++){fscanf(fp,"%lf",c+k*2+i*4+j);}}}fclose(fp);return c;}void savedata(int hang,double *data){FILE *fp;char *file1="地面点坐标数据.txt";fp=fopen(file1,"w");fprintf(fp,"---------像点对应地面点坐标为--------：\n");fprintf(fp,"\n");for (int i=0;i<hang;i++){fprintf(fp,"第%d点：",i+1);for (int j=0;j<3;j++){fprintf(fp,"%7.4lf ",data[i*3+j]);}fprintf(fp,"\n\n");}fprintf(fp,"-----------------------------------------");fclose(fp);return;}void main(){double *imagepoint;int row;int i,j;imagepoint=readdata();row=sizeof(imagepoint);//--------------------------------------------double f=24;f/=1000;doublefai1=-0.0061,oumiga1=0.0327,kapa1=0.1711,Ys1=397367.171,Xs1=3445820.098,Zs1=1486.212;doublefai2=0.0063,oumiga2=0.0178,kapa2=0.1489,Ys2=397367.234,Xs2=3445959.266,Zs2=1490.096;// printf("请输入左像片的外方位元素：\n");//printf("Xs1= ");//scanf("%lf",&Xs1);//printf("Ys1= ");//scanf("%lf",&Ys1);//printf("Zs1= ");//scanf("%lf",&Zs1);//printf("fai1= ");//scanf("%lf",&fai1);//printf("oumiga1= ");//scanf("%lf",&oumiga1);//printf("kapa1= ");//scanf("%lf",&kapa1);//printf("请输入右像片的外方位元素：\n");//printf("Xs2= ");//scanf("%lf",&Xs2);//printf("Ys2= ");//scanf("%lf",&Ys2);//printf("Zs2= ");//scanf("%lf",&Zs2);//printf("fai2= ");//scanf("%lf",&fai2);//printf("oumiga2= ");//scanf("%lf",&oumiga2);//printf("kapa2= ");//scanf("%lf",&kapa2);double Bx=Xs2-Xs1,By=Ys2-Ys1,Bz=Zs2-Zs1;double N1=0,N2=0;double X1=0,Y1=0,Z1=0,X2=0,Y2=0,Z2=0;double R1[3][3]={0.0};double R2[3][3]={0.0};double GEOdata[4][3]={0.0};for (i=0;i<row;i++){//---------组成左影像旋转矩阵--------R1[0][0]=cos(fai1)*cos(kapa1)-sin(fai1)*sin(oumiga1)*sin(kapa1);R1[0][1]=-cos(fai1)*sin(kapa1)-sin(fai1)*sin(oumiga1)*cos(kapa1);R1[0][2]=-sin(fai1)*cos(oumiga1);R1[1][0]=cos(oumiga1)*sin(kapa1);R1[1][1]=cos(oumiga1)*cos(kapa1);R1[1][2]=-sin(oumiga1);R1[2][0]=sin(fai1)*cos(kapa1)+cos(fai1)*sin(oumiga1)*sin(kapa1);R1[2][1]=-sin(fai1)*sin(kapa1)+cos(fai1)*sin(oumiga1)*cos(kapa1);R1[2][2]=cos(fai1)*cos(oumiga1);//-----------------------------------//---------组成右影像旋转矩阵--------R2[0][0]=cos(fai2)*cos(kapa2)-sin(fai2)*sin(oumiga2)*sin(kapa2);R2[0][1]=-cos(fai2)*sin(kapa2)-sin(fai2)*sin(oumiga2)*cos(kapa2);R2[0][2]=-sin(fai2)*cos(oumiga2);R2[1][0]=cos(oumiga2)*sin(kapa2);R2[1][1]=cos(oumiga2)*cos(kapa2);R2[1][2]=-sin(oumiga2);R2[2][0]=sin(fai2)*cos(kapa2)+cos(fai2)*sin(oumiga2)*sin(kapa2);R2[2][1]=-sin(fai2)*sin(kapa2)+cos(fai2)*sin(oumiga2)*cos(kapa2);R2[2][2]=cos(fai2)*cos(oumiga2);//-------------像空辅系坐标-------------X1=R1[0][0]*imagepoint[i*4+0]+R1[0][1]*imagepoint[i*4+1]-R1[0][2]*f;Y1=R1[1][0]*imagepoint[i*4+0]+R1[1][1]*imagepoint[i*4+1]-R1[1][2]*f;Z1=R1[2][0]*imagepoint[i*4+0]+R1[2][1]*imagepoint[i*4+1]-R1[2][2]*f;X2=R2[0][0]*imagepoint[i*4+2]+R2[0][1]*imagepoint[i*4+3]-R2[0][2]*f;Y2=R2[1][0]*imagepoint[i*4+2]+R2[1][1]*imagepoint[i*4+3]-R2[1][2]*f;Z2=R2[2][0]*imagepoint[i*4+2]+R2[2][1]*imagepoint[i*4+3]-R2[2][2]*f;//--------------------------------------//------------点投影系数-------------N1=(Bx*Z2-Bz*X2)/(X1*Z2-Z1*X2);N2=(Bx*Z1-Bz*X1)/(X1*Z2-Z1*X2);//-----------------------------------//------------计算地面点坐标------------GEOdata[i][0]=Xs1+N1*X1;GEOdata[i][1]=Ys1+By+N2*Y2;GEOdata[i][2]=Zs1+N1*Z1;//--------------------------------------}//--------------------------------------------for (i=0;i<4;i++){printf("第%d个地面点坐标：",i+1);for (j=0;j<3;j++){printf("%lf ",GEOdata[i][j]);}printf("\n\n");}savedata(row,GEOdata[0]);system("pause");}测试结果：3.单模型光束法严密平差缺少已知数据进行验证，因此如果有已知数据请代入已知数据进行验证。
单像空间后方交会
摄影测量学实习报告遥感07011班吴倩200732590254一、实习目的1.掌握空间后方交会的定义和实现算法(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。
(2)算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。
实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。
2.了解摄影测量平差的基本过程(1)获取已知数据。
从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。
(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
(3)确定未知数的初始值。
单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。
或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。
(4)计算旋转矩阵R。
利用角元素近似值计算方向余弦值，组成R阵。
(5)逐点计算像点坐标的近似值。
利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。
(6)逐点计算误差方程式的系数和常数项，组成误差方程式。
(7)计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。
(8)解求外方位元素。
根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。
(9)检查计算是否收敛。
将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.1′，当3个改正数均小于0.1′时，迭代结束。
实习一 单张影像空间后方交会程序设计
实习一单张影像空间后方交会程序设计一、实习目的
学生通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，为该门课程的学习打下好的基础。
二、实习内容
根据实习要求及试验数据实现单张影像空间后方交会的计算，并将最后计算得到的旋转矩阵及外方位元素输出。
图１空间后方交会程序流程图
三、预备知识
1、熟悉一门编程语言；
2、分析已有数据的特点；
3、弄清矩阵的存储与操作方法；
4、矩阵的运算模块（函数）可直接在相关资料上查找并直接使用。
四、实习步骤
已有的控制点像点坐标及物方坐标可通过输入界面进行输入，已可定义成文本文件由程序读取文本文件；定义变量和数组来存储数据；最终计算结果的输出（可输出为文本文件）。
实习结束后提交每人提交一份实习报告。
五、思考题
1、矩阵在计算机中是怎样存储的？怎样实现对矩阵中每个元素的操作（读取与存储）？
2、单张影像空间后方交会计算时在什么条件下迭代结束？
六、实验数据
已知航摄仪的内方位元素：f=153.24mm，x0=y0=0.0mm，摄影比例尺为1:50000；4个地面控制点的地面坐标及其对应像点的像片坐标如下表。
试验一 单片空间后方交会程序设计
三、存储数据的变量及数组
针对空间后方交会计算中数据的特点，需定义相 关的变量和数组来存储数据。 存储已知数据的数组：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 否 所有点完否？
摄影测量实习报告-单片空间后方交会8页word文档
摄影测量实习报告实习内容:单片空间后方交会编程实习者:李友兵学号:0810050121指导老师:张金平老师实习时间:2019.05.30——2019.06.03一、实习目的与任务此次摄影测量实习主要是要自主编程实现单像空间后方交会，通过已知的内方位元素和控制点像点坐标和地面坐标求解六个外方位元素，在此过程中深入理解单像空间后方交会的原理和对编程的熟悉和理解（我用的是C语言编程），和对时间的合理运用，对知识的综合运用，培养理论的实际运用能力，任务是在一个星期内自主完成。
二、单片空间后方交会理论基础单像空间后方交会：是通过以像点平面坐标为观测值，以控制点坐标为已知值，利用共线条件方程和最小二乘原理，运用间接平差方法，通过迭代求解6个外方位元素。
程序设计的思路与流程三、程序设计的思路与流程（编程框架和步骤）1、根据内方位元素和已知的数据将控制点的框标坐标转换为像平面坐标系坐标：x=x′-x0，y=y′-y02、确定未知数的初始值：角元素初始值κ0=ω0=Φ0=0；线元素初始值Zs0=H=mf，Xs0=(X1+X2+X3+X4)/4,Ys0=（Y1+Y2+Y3+Y4）/43、利用角元素初始值计算方向余弦值组成旋转矩阵Ra1=cosΦ*cosκ-sinΦ*sinω*sinκ,a2=-cosΦ*sinκ-sinΦ*sin ω*cosκ，a3=-sinΦ*cosωb1=cosω*sinκ,b2=cosω*cosκ,b3=-sinωc1=sinΦ*cosκ+cosΦ*sinω*sinκ,c2=-sinΦ*sinκ+cosΦ*sin ω*cosκ,c3=cosΦ*cosω4、逐点计算控制点的像点坐标的近似值，共线方程为:x=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys) +c3*(Z-Zs))y=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys) +c3*(Z-Zs))(x)=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Y s)+c3*(Z-Zs))(y)=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Y s)+c3*(Z-Zs))5、组成误差方程：Vx=a11*dXs+a12*dYs+a13*dZs+a14*dΦ+a15*dω+a16*dκ+(x)-xVy=a21*dXs+a22*dYs+a23*dZs+a24*dΦ+a25*dω+a26*dκ+(y)-y系数求解(是共线方程分别外方位元素求导，是共线方程线性化的系数)：变量代换A= a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)B= a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)C= a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)a11=(a1*f+a3*x)/C,a12=(b1*f+b3*x)/C,a13=(c1*f+c3*f)/C,a14=y *sinω-((x*(x*cosκ-y*sinκ))/f+f*cosκ)cosω,a15=-f*sin κ-(x*(x*sinκ+y*cosκ)/f),a16=ya21=(a2*f+a3*x)/C,a22=(b2*f+b3*x)/C,a23=（c2*f+c3*x）/C,a24=-xsinω-((x*(x*cosκ-y*sinκ))/f-f*sinκ)cosω,a25=-f*cosκ-(y*(x*sinκ+y*cosκ)/f),a26 =-x 误差方程的常系数（是像点坐标观测值与计算的近似值的差值）：lx=x-(x) ,ly=y-(y)6、组成法方程，解求外方位元素改正数X=(A T A)-1A T L(A为误差方程的系数矩阵，L为误差方程的常系数矩阵，通过步骤5求得，此处先求A T A再求矩阵的逆矩阵，解得的改正数加上相应的近似值得到外方位元素新的近似值)7、检查计算是否收敛：将求得的外方位元素改正数与规定的限差比较，大于限差继续迭代，小于限差则终止。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
{System;System.Collections.Generic;System.Linq;System.Text;namespace 单像空间后方交会{class Program{static void Main( string [] args)for (j = 0; j < 5; j++)if (j < 3)"请输入第 {0} 个点的第 {1} 个地面坐标： ", i + 1, j+ 1); double .Parse( Console .ReadLine());"请输入第 {0} 个点的第 {1} 个像点坐标： ", i + 1, j - 2);double .Parse( Console .ReadLine());Console .WriteLine();// 归算像点坐标(i = 0; i < 4; i++)for (j = 3; j < 5; j++)if (j == 3)zuobiao[i, j] = zuobiao[i, j] - x0;else zuobiao[i, j] = zuobiao[i, j] - y0;// 计算和确定初值double zs0 = m * f, xs0 = 0, ys0 = 0; for (i = 0; i < 4; i++)elseusing using using using x0 = y0 = int x0, y0, i, j; double f, m;Console .Write( " 请输入像片比例尺： ");double .Parse( Console .ReadLine());Console .Write( " 请输入像片的内方位元素 x0:" ); // 均以毫米为单位 int .Parse( Console .ReadLine());Console .Write( " 请输入像片的内方位元素 y0:" );int .Parse( Console .ReadLine());Console .Write( " 请输入摄影机主距 f:" );double .Parse( Console .ReadLine());Console .WriteLine();// 输入坐标数据 double [,] zuobiao = new double [4, 5];(i = 0; i < 4; i++)forConsole .Write( zuobiao[i, j] = Console .Write( zuobiao[i, j] = forxs0 = 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;double [,] r = double [,] jinsi = double [] chazhi = double [] jieguo ={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));int q = 0;new double [3, 3]; new double [4, 2];new double [8];new double [6];r[0, 0] = Math.Sin(kappa0);r[0, 1] = -Math.Cos(kappa0); r[0, 2] = - r[1, 0] = r[1, 1] = r[1, 2] = - r[2,0] = Math.Sin(kappa0);r[2, 1] = - double [,] zhong = matrixChe(dReturn,matrixTrans(xishu)); do// 计算旋转矩阵 r Math.Cos(phi0) *Math.Cos(kappa0) - Math.Sin(phi0) *Math.Sin(omega0) * Math.Cos(phi0) * Math.Sin(kappa0) - Math.Sin(phi0) * Math.Sin(omega0) * Math.Sin(phi0) * Math.Cos(omega0) * Math.Cos(omega0) * Math.Sin(omega0);Math.Sin(phi0) *Math.Cos(omega0);Math.Sin(kappa0);Math.Cos(kappa0);Math.Cos(kappa0) + Math.Cos(phi0) * Math.Sin(omega0) * Math.Sin(phi0) * Math.Sin(kappa0) + Math.Cos(phi0) *Math.Sin(omega0) *Math.Cos(phi0) * //计算x,y 的近似值for (i = 0; i < 4; i++)Math.Cos(omega0);}for ( int j = 0; j < C.GetLength(1); j++){}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[2]) > 0.020); Math.Abs(jieguo[1]) > 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; ( int i = 0; i < C.GetLength(0);i++)for Console .Write( " {0}" , C[i, j]);for ( int s = Level - 1; s > i; s--){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, 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;k = 1;( int i = 0, j = 0; i < Level && j < Level;i++, j++)int for if (dMatrix[i, j] ==0)int m = i;for (; dMatrix[m, j] == 0;m++) ; if (m == Level)return 0;elsefor ( int n = j; n < Level;n++){c = dMatrix[i, n];dMatrix[i, n] = dMatrix[m,n]; dMatrix[m, n] = c;}k *= (-1);x = dMatrix[s, j];for ( int t = j; t < Level;t++)int Level)} for ( int s = Level - 1; s > i; s--){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];else return 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]; else dReverseMatrix[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; x = dReverseMatrix[s,j];} }for ( int t = j; t < 2 * Level; t++) dReverseMatrix[s, t] -= (dReverseMatrix[i, t] * x);( 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;for。