空间后方交会编程实习报告
一实习目的
用程序设计语言（Visual C++或者C语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
二 实习内容
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。
三 实习数据
已知航摄仪的内方位元素：fk=153.24mm，x0=y0=0.0mm，摄影比例尺为1:50000；
4个地面控制点的地面坐标及其对应像点的像片坐标：
四实习原理
如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。
单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，t，w，k。
五 实习流程
（1）获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。
（2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
（3）确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：
Z0s=H=m*f+ΣZi/4；
X0s=Σxi/n;
Y0s=ΣYi/n;
t=ω=κ=0;
式中：m为摄影比例尺分母；
(4)计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。
（5）逐点计算像点坐标的近似值。利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值（x）、（y）。
（6）按下式逐点计算误差方程式的系数和常数项，组成误差方程。
(7)计算法方程的系数矩阵ATA与常数项ATL，组成法方程；
（8）解求外方位元素。根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。
（9）检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对t、ω、κ、Xs、Ys、Zs的改正数Δt，Δω，Δκ，ΔXs，ΔYs，ΔZs给予限差，当改正数小于限差时，迭代结束。否则用新的近似值重复（4）——（8）步骤计算，直到满足要求为止。
（10）空间后方交会的精度估计：
按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵（ATA）-1来解求，此时视像点坐标为等精度不相关观测值。因为ATA）-1中第i个主对角线上的元素Qii就是法方程式中第i个未知数的权倒数，若单位权中误差为m0，则第i个未知数的中误差为：
mi=
0
当参加空间后方交会的控制点有n个时，则单位权中误差可按下式计算：
六 主要代码与详解
void R(double t,double w,double k,double *a,double *b,double *c)
{
a[0]=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);
a[1]=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k);
a[2]=-sin(t)*cos(w);
b[0]=cos(w)*sin(k);
b[1]=cos(w)*cos(k);
b[2]=-sin(w);
c[0]=sin(t)*cos(k)+cos(t)*sin(w)*sin(k);
c[1]=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k);
c[2]=cos(t)*cos(w);
}
//子函数计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。
void main()
{
int i,m,num;
double t,w,k,Xs,Ys,Zs,f; //六个外方位元素与焦距f
double x[N]={-86.15,-53.40,-14.78,10.46},y[N]={-68.99,82.21,-76.63,64.43};
//4个像点坐标
double X[N]={36589.41,37631.08,39100.97,40426.54},Y[N]={25273.32,31324.51,24934.98,30319.81},Z[N]={2195.17,728.69,2386.50,757.31};//四个控制点的空间坐标
m=50000;f=153.24;// 影像比例尺1/m 与焦距f
// 以上主要是已知数据的定义
double V[6]={0},a[3],b[3],c[3],Xo[N],Yo[N],Zo[N],A[6*M],B[6*M],l[M],C[36],D[6],E[8];
for(i=0;i<N;i++)
{
x[i]=x[i]/1000.0;
y[i]=y[i]/1000.0;//像点坐标单位换算
}
// 以下为确定未知数的初始值。
Xs=(X[0]+X[1]+X[2]+X[3])/N;
Ys=(Y[0]+Y[1]+Y[2]+Y[3])/N;
Zs=m*f+(Z[0]+Z[1]+Z[2]+Z[3])/N;// Xs，Ys，Zs 为摄站点的空间坐标初始值
f=f/1000.0;
t=w=k=0.0; //角元素的近似值
for(num=1;num<10;num++) //num控制循环迭代次数
{
R(t,w,k,a,b,c);// 计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。
for(i=0;i<N;i++)//下面的循环用来逐点计算误差方程式的系数和常数项，组成误差方程。
{
//以下是用共线条件方程计算控制点像点坐标的近似值（x）、（y）。 Xo[i]=-f*(a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs)); Yo[i]=-f*(a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs));
Zo[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);//Zo便于后面计算
A[12*i+0]=(a[0]*f+a[2]*x[i])/Zo[i];
A[12*i+1]=(b[0]*f+b[2]*x[i])/Zo[i];
A[12*i+2]=(c[0]*f+c[2]*x[i])/Zo[i];
A[12*i+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);
A[12*i+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;
A[12*i+5]=y[i];
A[12*i+6]=(a[1]*f+a[2]*y[i])/Zo[i];
A[12*i+7]=(b[1]*f+b[2]*y[i])/Zo[i];
A[12*i+8]=(c[1]*f+c[2]*y[i])/Zo[i]; A[12*i+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);
A[12*i+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;
A[12*i+11]=-x[i];
l[2*i]=x[i]-Xo[i];
l[2*i+1]=y[i]-Yo[i];//计算l
}
zhuanzhi(A,B,8,6); //转置A
xiangchen(B,A,C,6,8,6);//求ATA
xiangchen(B,l,D,6,8,1); //求常数项ATL
qiuni(C,6);//求ATA的逆
xiangchen(C,D,V,6,6,1);//求改正数V
Xs+=V[0];
Ys+=V[1];
Zs+=V[2];
t+=V[3];
w+=V[4];
k+=V[5]; //结果改正
if((fabs(V[0])<=0.00001)&&(fabs(V[1])<=0.00001)&&(fabs(V[2])<=0.00001)&&(fabs(V[3])<=0.00001)&&(fabs(V[4])<=0.00001)&&(fabs(V[5])<=0.00001))
break; //限差控制。检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较。当改正数小于限差时，迭代结束。否则用新的近似值重复计算，直到满足要求为止 。
}
//以下是计算单位权中误差
double s=0,m0,v[8];
xiangchen(A,V,E,8,6,1); //计算AX
for(i=0;i<8;i++)
{
v[i]=E[i]-l[i];//计算AX-L
s+=v[i]*v[i];//计算[vv]
}
m0=sqrt(s/2);// 单位权中误差按下式计算：，这里n=4.
R(t,w,k,a,b,c);
//输出结果
printf("像主点的空间坐标为：\n");
printf("Xs=%.2f\n",Xs);
printf("Ys=%.2f\n",Ys);
printf("Zs=%.2f\n",Zs);
printf("旋转矩阵R为:\n");
for(i=0;i<3;i++)
printf(" %.5f",a[i]);
printf("\n");
for(i=0;i<3;i++)
printf(" %.5f",b[i]);
printf("\n");
for(i=0;i<3;i++)
printf(" %.5f",c[i]);
printf("\n");
printf("单位权中误差为:",m0);
cout<<m0<<endl;
}
七 实习结果
八 实习总结
通过本次实习我深刻地理解单片空间后方交会的原理，尤其是共线方程。学会了在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
这次实习我是先在课后花时间去掌握空间后方交会的原理与计算过程，从整体上去把握整个过程，然后再去编程实现这个过程。
实习的难点在于程序中的变量过多，关系式较多，需要好好去理清，冷静去处理。可以说细节决定了实习的成败，一个小小的等式错误，不仅结果错误，而且很难去检查出来。
实习中我还学会了利用断点，单步执行和变量查看来快速的定位错误和找到错误的原因。可以说这种VC调试手段是很有用的。
总之，这次实习加强了我综合运用所学知识解决实际问题的能力。
第二篇：单像空间后方交会实习报告
摄影测量学
单像空间后方交会实习报告
一、 实习目的
1. 掌握空间后方交会的定义和实现算法
(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。
(2)算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。
2. 了解摄影测量平差的基本过程
(1) 获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标X，Y，Z。
(2) 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
(3) 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。
(4) 计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。
(5) 逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。
(6) 逐点计算误差方程式的系数和常数项，组成误差方程式。
(7) 计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。
(8) 解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。
(9) 检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.000001弧度，当3个改正数均小于0.000001弧度时，迭代结束。否则用新的近似值重复(4)～（8）步骤的计算，直到满足要求为止。
3. 通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。
4. 实习过程：4.1学习单张像片空间后方交会的基本理论，掌握其基本思想。如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。而单像空间后方交会就是用于测定像片的外方位元素的，它的基本思想是：以单幅影像为基础，从影像所覆盖的地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,p,w,k.由于空间后方交会所采用的数学模型共线方程是非线性函数，为了便于外方位元素的解求，首先将其线性化。4.2在纸上绘出空间后方交会的计算机程序框图。为了能够在宏观上指导我们编写程序，我们需要在草稿纸上绘出程序框图。
二、 源代码
using System;
usingSystem.Collections.Generic;
usingSystem.Linq;
usingSystem.Text;
namespace 单像空间的后方交会
{
public class calculate
{
private double j, k, l, Xs, Ys, Zs;//六个外方位元素
private double f = 28.1539;//主距
struct point //像点和地面点坐标
{
public double x, y, X, Y, Z;
}
private point[] p = new point[4];//存取控制点的坐标
private double[] R = new double[9];//旋转矩阵
private double[] a = new double[8];//近似值坐标
private double[] L = new double[8];//误差方程常数项
private double[,] A = new double[8, 6];//误差方程系数项
private int count = 0;
public void Y(double[] q)
{
for (int n = 0; n < 4; n++)
{
int m = n * 5;
Console.WriteLine("请输入第｛0｝控制点的坐标", n + 1);
q[m] = Convert.ToDouble(Console.ReadLine());
q[m + 1] = Convert.ToDouble(Console.ReadLine());
q[m + 2] = Convert.ToDouble(Console.ReadLine());
q[m + 3] = Convert.ToDouble(Console.ReadLine());
q[m + 4] = Convert.ToDouble(Console.ReadLine());
p[n].x = q[m];
p[n].y = q[m + 1];
p[n].X = q[m + 2];
p[n].Y = q[m + 3];
p[n].Z = q[m + 4];
}
double ave = 0, sum = 0;// 求比例尺分母，用来求外方位元素
for (int n = 0; n < 3; n++)
{
for (int m = n + 1; m < 4; m++)
{
sum += Math.Sqrt(Math.Pow(p[n].X - p[m].X, 2) + Math.Pow(p[n].Y - p[m].Y, 2)) / Math.Sqrt(Math.Pow(p[n].x - p[m].x, 2) + Math.Pow(p[n].y - p[m].y, 2));
}
}
ave = sum / 6;
double j = 0.054882; k = 0.057034; l = -0.036175;//六个外方位元素的近似值
doubleXs = 500215.49;
doubleYs = 4185301.89;
doubleZs = 1475.56;
}
private double sin(double m)
{
returnMath.Sin(m);
}
private double cos(double m)
{
returnMath.Cos(m);
}
public void calX()//计算旋转矩阵
{
R[0] = cos(j) * cos(k) - sin(j) * sin(l) * sin(k);
R[1] = -cos(j) * sin(k) - sin(j) * sin(l) * cos(k);
R[2] = -sin(j) * cos(l);
R[3] = cos(l) * sin(k);
R[4] = cos(l) * cos(k);
R[5] = -sin(l);
R[6] = sin(j) * cos(k) + cos(j) * sin(l) * sin(k);
R[7] = -sin(j) * sin(k) + cos(j) * sin(l) * cos(k);
R[8] = cos(j) * cos(l);
}
public void calJSZ()//计算像点坐标近似值
{
for (int n = 0; n < 4; n++)
{
a[2 * n] = -f * (R[0] * (p[n].X - Xs) + R[3] * (p[n].Y - Ys) + R[6] * (p[n].Z - Zs)) / (R[2] * (p[n].X - Xs) + R[5] * (p[n].Y - Ys) + R[8] * (p[n].Z - Zs));
a[2 * n + 1] = -f * (R[1] * (p[n].X - Xs) + R[4] * (p[n].Y - Ys) + R[7] * (p[n].Z - Zs)) / (R[2] * (p[n].X - Xs) + R[5] * (p[n].Y - Ys) + R[8] * (p[n].Z - Zs));
}
}
public void calXSJZandCSX()//计算系数矩阵和常数项
{
for (int n = 0; n < 4; n++)//计算常数项
{
L[2 * n] = p[n].x - a[2 * n];
L[2 * n + 1] = p[n].y - a[2 * n + 1];
}
for (int n = 0; n < 4; n++)//计算系数矩阵
{
double z = R[2] * (p[n].X - Xs) + R[5] * (p[n].Y - Ys) + R[8] * (p[n].Z - Zs);
A[2 * n, 0] = (R[0] * f + R[2] * p[n].x) / z;
A[2 * n, 1] = (R[3] * f + R[5] * p[n].x) / z;
A[2 * n, 2] = (R[6] * f + R[8] * p[n].x) / z;
A[2 * n, 3] = p[n].y * sin(l) - (p[n].x * (p[n].x * cos(k) - p[n].y * sin(k)) / f + f * cos(k)) * cos(l);
A[2 * n, 4] = -f * sin(k) - p[n].x / f * (p[n].x * sin(k) + p[n].y * cos(k));
A[2 * n, 5] = p[n].y;
A[2 * n + 1, 0] = (R[1] * f + R[2] * p[n].y) / z;
A[2 * n + 1, 1] = (R[4] * f + R[5] * p[n].y) / z;
A[2 * n + 1, 2] = (R[7] * f + R[8] * p[n].y) / z;
A[2 * n + 1, 3] = p[n].x * sin(l) - (p[n].y * (p[n].x * cos(k) - p[n].y * sin(k)) / f - f * sin(k)) * cos(l);
A[2 * n + 1, 4] = -f * cos(k) - p[n].y / f * (p[n].x * sin(k) + p[n].y * cos(k));
A[2 * n, 5] = -p[n].x;
}
}
public double calJSGZS()//计算改正数
{
double[,] AT = new double[6, 8];//A转置
double[,] temp = new double[6, 6];//A转置与A的乘积
double[] X = new double[6];//改正数
double[] ATL = new double[6];//A转置与L相乘的积
int n, m, s;
for (n = 0; n < 8; n++)//求A的转置矩阵AT
for (m = 0; m < 6; m++)
{
AT[m, n] = AT[n, m];
}
for (n = 0; n < 6; n++)//求A转置与A的乘积
for (m = 0; m < 6; m++)
{
temp[n, m] = 0;
for (s = 0; s < 8; s++)
{
temp[n, m] += AT[n, s] * A[s, m];
}
}
MA(temp);//求逆
for (n = 0; n < 6; n++)//求A的转置与L的乘积
{
for (m = 0; m < 8; m++)
{
ATL[n] += AT[n, m] * L[m];
}
}
for (n = 0; n < 6; n++)//求改正数X
{
for (m = 0; m < 6; m++)
{
X[n] += temp[n, m] * ATL[m];
}
}
Xs += X[0];
Ys += X[1];
Zs += X[2];
j += X[3];
k += X[4];
l += X[5];
returnXs;
returnYs;
returnZs;
return j;
return k;
return l;
}
public void Main(string[] args)//计算流程n
{
while (count >= 4)
{
calX();
calJSZ();
calXSJZandCSX();
Console.WriteLine("第{0}次迭代的结果", count);
calJSGZS();
count++;
}
}
private void MA(double[,] c)
{
int i, j, h, m;
constint n = 6;
double l;
double[,] q = new double[n, 12];
for (i = 0; i < n; i++)//求高斯矩阵
{
for (j = 0; j < 12; j++)
{
q[i, j] = c[i, j];
}
}
for (i = 0; i < n; i++)//
{
for (j = 0; j < 12; j++)//构造单位阵
{
if (i + 6 == j)
q[i, j] = 1;
else
q[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];
}
}
}
}
}
三、 计算结果
起算数据文件为：
文件存储的计算结果为：
