单片空间后方交会程序设计

合集下载

空间后方交会程序设计

单片空间后方交会程序设计专业：测绘工程班级：测绘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。

第五讲单片空间后方交会

x12 − f (1 + 2 ) f xy − 1 1 f
2 x2 − f (1 + 2 ) f
−
x1 y1 f
y12 − f (1 + 2 ) f − x2 y2 f
x y − 2 2 f
2 x3 − f (1 + 2 ) f
2 y2 − f (1 + 2 ) f
−
x3 y3 f
xy − 3 3 f
Y B
A
C X
利用航摄像片上三个以上像点坐标和对应像点坐标和对应地面点坐标，计算像片外方位元素的工作，称为单张像片的空间后方交会。进行空间后方交会运算，常用的一个基本公式是前面提到的共线方程。式中的未知数，是六个外方位元素。由于一个已知点可列出两个方程式，如有三个不在一条直线上的已知点，就可列出六个独立的方程式，解求六个外方位元素。由于共线条件方程的严密关系式是非线性函数，不便于计算机迭代计算。为此，要由严密公式推导出一次项近似公式，即变为线性函数。
（5）用所取未知数的初始值和控制点的地面坐标，代入共线方程式，逐）用所取未知数的初始值和控制点的地面坐标，代入共线方程式，点计算像点坐标的近似值 ( x), ( y ) 并计算 lx , l y a ( X − X S ) + b1 (Y − YS ) + c1 ( Z − Z S ) x=−f 1 a3 ( X − X S ) + b3 (Y − YS ) + c3 ( Z − Z S ) a ( X − X S ) + b2 (Y − YS ) + c2 ( Z − Z S ) y=−f 2 a3 ( X − X S ) + b3 (Y − YS ) + c3 ( Z − Z S ) （6）组成误差方程式。）组成误差方程式。 7）计算法方程式的系数矩阵与常数项，组成法方程式。（7）计算法方程式的系数矩阵与常数项，组成法方程式。（8）解算法方程，迭代求得未知数的改正数。）解算法方程，迭代求得未知数的改正数。

单像空间后方交会程序《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．实验前认真复习单片空间后方交会的有关内容。

六、实验报告
原理、实验步骤及程序框图
上机调试程序并打印结果。

单向空间后方交会报告

单向空间后方交会程序设计报告课题名称计算单向空间后方交会姓名胡小伟学号 2009301610207学院、班级测绘学院0906班方向A2011年 10月18日※※※※※※※※※※※※※※※ 单向空间后方交会程序设计报告（1）编程工具：Matlab R2010a（2）程序代码：已知值及初始值控制点：x(mm) y(mm)-11.00 -13.8010.48 -6.4410.78 8.3122.712 8.32X(m) Y(m) Z(m)500017.19 4185066.00 885.40500452.19 4185197.50 915.82500470.16 4185494.50 928.76500302.16 4185489.00 949.71焦距：28.1359（mm）外方位元素初值：X(m) Y(m) Z(m) phi(rad) omega(rad) kappa(rad) 500215.49 4185301.89 1475.56 0.054882 0.057034 -0.036175%单向空间后方交会clc;%主距mf=28.1359/1000;%像点坐标mxy = [-11.00 -13.8010.48 -6.4410.78 8.3122.712 8.32]/1000;%地面坐标mXYZ = [500017.19 4185066.00 885.40;500452.19 4185197.50 915.82;500470.16 4185494.50 928.76;500302.16 4185489.00 949.71];%比例尺m=(sqrt((XYZ(1,1)-XYZ(2,1))^2+(XYZ(1,2)-XYZ(2,2))^2))/...(sqrt((xy(1,1)-xy(2,1))^2+(xy(1,2)-xy(2,2))^2))%初值XS0=500215.49YS0=4185301.89ZS0=1475.56H=m*fFai0=0.054882Oum0=0.057034Kpa0=-0.036175a1=cos(Fai0)*cos(Kpa0)-sin(Fai0)*sin(Oum0)*sin(Kpa0);a2=-cos(Fai0)*sin(Kpa0)-sin(Fai0)*sin(Oum0)*cos(Kpa0);a3=-sin(Fai0)*cos(Oum0);b1=cos(Oum0)*sin(Kpa0);b2=cos(Oum0)*cos(Kpa0);b3=-sin(Oum0);c1=sin(Fai0)*cos(Kpa0)+cos(Fai0)*sin(Oum0)*sin(Kpa0);c2=-sin(Fai0)*sin(Kpa0)+cos(Fai0)*cos(Oum0)*cos(Kpa0);c3=cos(Fai0)*cos(Oum0);for i=1:4f=28.1359/1000;x(i)=-f*(a1*(XYZ(i,1)-XS0)+b1*(XYZ(i,2)-YS0)+c1*(XYZ(i,3)-ZS0))/...(a3*(XYZ(i,1)-XS0)+b3*(XYZ(i,2)-YS0)+c3*(XYZ(i,3)-ZS0));y(i)=-f*(a2*(XYZ(i,1)-XS0)+b2*(XYZ(i,2)-YS0)+c2*(XYZ(i,3)-ZS0))/...(a3*(XYZ(i,1)-XS0)+b3*(XYZ(i,2)-YS0)+c3*(XYZ(i,3)-ZS0));lx(i)=xy(i,1)-x(i);ly(i)=xy(i,2)-y(i);endxyL=[lx(1);lx(1); lx(2);lx(2);lx(3);lx(3); lx(4);lx(4)]A1=[-f/H 0 -xy(1,1)/H -f*(1+(xy(1,1)^2)/(f^2)) -(xy(1,1))*(xy(1,2)/f) xy(1,2)0 -f/H -(xy(1,2))/H -(xy(1,1))*(xy(1,2)/f) -f*(1+(xy(1,2)^2)/(f^2)) -xy(1,2)]; A2=[-f/H 0 -xy(2,1)/H -f*(1+(xy(2,1)^2)/(f^2)) -(xy(2,1))*(xy(2,2)/f) xy(2,2)0 -f/H -(xy(2,2))/H -(xy(2,1))*(xy(2,2)/f) -f*(1+(xy(2,2)^2)/(f^2)) -xy(2,2)]; A3=[-f/H 0 -xy(3,1)/H -f*(1+(xy(3,1)^2)/(f^2)) -(xy(3,1))*(xy(3,2)/f) xy(3,2)0 -f/H -(xy(3,2))/H -(xy(3,1))*(xy(3,2)/f) -f*(1+(xy(3,2)^2)/(f^2)) -xy(3,2)]; A4=[-f/H 0 -xy(4,1)/H -f*(1+(xy(4,1)^2)/(f^2)) -(xy(4,1))*(xy(4,2)/f) xy(4,2)0 -f/H -(xy(4,2))/H -(xy(4,1))*(xy(4,2)/f) -f*(1+(xy(4,2)^2)/(f^2)) -xy(4,2)]; A=[A1;A2;A3;A4];X=(A'*A)\A'*LXS1=(XS0+X(1))YS1=(YS0+X(2))ZS1=(ZS0+X(3))Fai1=Fai0+X(4)Oum1=Oum0+X(5)Kpa1=Kpa0+X(6)XS1=(XS0+X(1))YS1=(YS0+X(2))ZS1=(ZS0+X(3))jishu=1;while (abs((Fai1-Fai0)>0.000003)||(abs(Oum1-Oum0)>0.000003)...||(abs(Kpa1-Kpa0))>0.000003)jishu=jishu+1XS0=XS1YS0=YS1ZS0=ZS1Fai0=Fai1Oum0=Oum1Kpa0=Kpa1a1=cos(Fai0)*cos(Kpa0)-sin(Fai0)*sin(Oum0)*sin(Kpa0);a2=-cos(Fai0)*sin(Kpa0)-sin(Fai0)*sin(Oum0)*cos(Kpa0);a3=-sin(Fai0)*cos(Oum0);b1=cos(Oum0)*sin(Kpa0);b2=cos(Oum0)*cos(Kpa0);b3=-sin(Oum0);c1=sin(Fai0)*cos(Kpa0)+cos(Fai0)*sin(Oum0)*sin(Kpa0);c2=-sin(Fai0)*sin(Kpa0)+cos(Fai0)*cos(Oum0)*cos(Kpa0);c3=cos(Fai0)*cos(Oum0);for i=1:4f=28.1359/1000;x(i)=-f*(a1*(XYZ(i,1)-XS0)+b1*(XYZ(i,2)-YS0)+c1*(XYZ(i,3)-ZS0))/...(a3*(XYZ(i,1)-XS0)+b3*(XYZ(i,2)-YS0)+c3*(XYZ(i,3)-ZS0));y(i)=-f*(a2*(XYZ(i,1)-XS0)+b2*(XYZ(i,2)-YS0)+c2*(XYZ(i,3)-ZS0))/...(a3*(XYZ(i,1)-XS0)+b3*(XYZ(i,2)-YS0)+c3*(XYZ(i,3)-ZS0));lx(i)=xy(i,1)-x(i);ly(i)=xy(i,2)-y(i);endxyL=[lx(1);lx(1); lx(2);lx(2);lx(3);lx(3); lx(4);lx(4)]A1=[-f/H 0 -xy(1,1)/H -f*(1+(xy(1,1)^2)/(f^2)) -(xy(1,1))*(xy(1,2)/f) xy(1,2)0 -f/H -(xy(1,2))/H -(xy(1,1))*(xy(1,2)/f) -f*(1+(xy(1,2)^2)/(f^2)) -xy(1,2)]; A2=[-f/H 0 -xy(2,1)/H -f*(1+(xy(2,1)^2)/(f^2)) -(xy(2,1))*(xy(2,2)/f) xy(2,2)0 -f/H -(xy(2,2))/H -(xy(2,1))*(xy(2,2)/f) -f*(1+(xy(2,2)^2)/(f^2)) -xy(2,2)]; A3=[-f/H 0 -xy(3,1)/H -f*(1+(xy(3,1)^2)/(f^2)) -(xy(3,1))*(xy(3,2)/f) xy(3,2)0 -f/H -(xy(3,2))/H -(xy(3,1))*(xy(3,2)/f) -f*(1+(xy(3,2)^2)/(f^2)) -xy(3,2)]; A4=[-f/H 0 -xy(4,1)/H -f*(1+(xy(4,1)^2)/(f^2)) -(xy(4,1))*(xy(4,2)/f) xy(4,2)0 -f/H -(xy(4,2))/H -(xy(4,1))*(xy(4,2)/f) -f*(1+(xy(4,2)^2)/(f^2)) -xy(4,2)]; A=[A1;A2;A3;A4];X=(A'*A)\A'*LFai1=Fai0+X(4)Oum1=Oum0+X(5)Kpa1=Kpa0+X(6)XS1=(XS0+X(1))YS1=(YS0+X(2))ZS1=(ZS0+X(3))end（3）运行结果：改正数：X =1.0e-003 *0.880693463596646-0.173862879363425-0.160454437823512-0.0014419260245700.0003069730686500.000107948847312外方位元素：Fai =0.054994013287714Oum =0.056571816845093Kpa =-0.036172713831137XS =5.002154276846778e+005YS =4.185302159274295e+006ZS =1.475497145016308e+003（4）心得体会通过本次编程，我对单向空间后方交会的算法有了直观、深入的了解，并对Matlab的使用有了进一步的认识。

摄影测量学空间后方交会实验报告

摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘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.评定精度：根据实验数据，评定求解得到的像片外方位元素的精度。

四、实验意义及能力培养单像空间后方交会实验有助于学生掌握摄影测量基本原理，了解单像空间后方交会的计算过程，提高动手实践能力。

通过实验，学生可以深入理解线性代数和微分学在摄影测量中的应用，为后续学习提供更扎实的基础。

此外，实验还可以培养学生的解决问题的能力和综合运用所学知识的能力，为未来从事相关领域工作打下坚实基础。

综上所述，单像空间后方交会是一种重要的摄影测量方法，通过程序设计实现对影像的定位和测量。

摄影测量学单像空间后方交会程序设计作业

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、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

单片空间后方交会程序设计
1 目的
用程序设计语言（VC或者VB）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。