空间后方交会程序.docx

实验目的：
掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。

仪器用具及已知数据文件：
计算机windowsxp系统，编程软件（VISUALC++6.0），地面控制点在摄影测量坐标系中的坐标及其像点坐标文件shuju.txt。

实验内容：
单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标根据共线方程反求影像的外方位元素。
数学模型：共线条件方程式：

求解过程：
（1）获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影测量坐标。
（2）量测控制点的像点坐标并做系统改正。
（3）确定未知数的初始值。在竖直摄影且地面控制点大致分布均匀的情况下，按如下方法确定初始值，即：
,,
φ=ω=κ=0
式中；m为摄影比例尺分母；n为控制点个数。
（4）用三个角元素的初始值，计算个方向余弦，组成旋转矩阵R。
（5）逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x）、（y）。
（6）逐点计算误差方程式的系数和常数项，组成误差方程式。
（7）计算法方程的系数矩阵和常数项，组成法方程式。
（8）解法方程，求得外方位元素的改正数,，，dφ,dω,dκ。
（9）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。
（10）将求得的外方位元素改正数与规定的限差比较，若小于限差则迭代结束。否则用新的近似值重复（4）～（9），直到满足要求为止。

四．实验结果：
程序的源代码如下所示：
#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include<math.h>
#include<conio.h>
#defineN4
voidturn(double*A,doubleA2[],intm,intn)//计算矩阵的转置
{
inti,j;
for(i=0;i<m;i++)
for(j=0;j<n;j++)
A2[j*m+i]=A[i*n+j];
}
voidmulAB(double*A,double*B,double*C,intam,intan,intbm,intbn)//计算两矩阵相乘
{
inti,j,l,u;
if(an!=bm)
{
printf("error!cannotdothemultiplication.\n");
return;
}
for(i=0;i<am;i++)
for(j=0;j<bn;j++)
{
u=i*bn+j;
C[u]=0.0;
for(l=0;l<an;l++)
C[u]+=A[i*an+l]*B[l*bn+j];
}
return;
}
double*inv(double*a,intn)//计算矩阵的逆，本程序的难点
{//采用高斯-约旦-全选主元法
int*is,*js,i,j,k,l,u,v;
doubled,p;
is=(int*)malloc(n*sizeof(int));
js=(int*)malloc(n*sizeof(int));
for(k=0;k<=n-1;k++)
{
d=0.0;
for(i=k;i<n;i++)
for(j=k;j<n;j++)
{l=i*n+j;
p=fabs(a[l]);
if(p>d)
{
d=p;
is[k]=i;
js[k]=j;
}
}
if(d+1.0==1.0)
{free(is);
free(js);
printf("errornotinv\n");
returnNULL;
}
if(is[k]!=k)
for(j=0;j<n;j++)
{u=k*n+j;
v=is[k]*n+j;
p=a[u];
a[u]=a[v];
a[v]=p;
}
if(js[k]!=k)
for(i=0;i<n;i++)
{u=i*n+k;
v=i*n+js[k];
p=a[u];
a[u]=a[v];
a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for(j=0;j<n;j++)
if(j!=k)
{
u=k*n+j;
a[u]=a[u]*a[l];
}
for(i=0;i<n;i++)
if(i!=k)
for(j=0;j<n;j++)
if(j!=k)
{u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for(i=0;i<n;i++)
if(i!=k)
{u=i*n+k;
a[u]=-a[u]*a[l];
}
}
f