while((XX[3]>6/206265 || XX[4]>6/206265 || XX[5]>6/206265)&&n<100)
{
/*----------------旋转矩阵R-----------------------*/
a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);
double X[N],Y[N],Z[N];
double x1[N],y1[N];
double m;
double L[2*N];
double XX[6];
double A[2*N][6];
double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];
{if(q[i][h]==0)
continue;
p=q[k][h]/q[i][h];
for(j=0;j<12;j++)
{q[i][j]*=p;
q[i][j]-=q[k][j];}}
for(i=0;i<n;i++)//将对角线上数据化为1
{ p=1.0/q[i][i];
for(j=0;j<12;j++)
c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);
c[2]=cos(q)*cos(w);
/*-----------------像点坐标计算值------------------*/
for(i=0;i<N;i++)
{
X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs);
for(i=0;i<N;i++)
{
A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];
A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i];
A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i];
A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin(k))/f+f*cos(k))
*cos(w);
A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;
A[2*i][5]=y[i]-y0;
L[2*i]=x[i]-x1[i];
A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];
int i,n=0;
double sum=0,m0;
/*---------------输入点地面坐标---------------------*/
X[0]=5083.205;X[1]=5780.02;X[2]=5210.879;X[3]=5909.264;
Y[0]=5852.099;Y[1]=5906.365;Y[2]=4258.446;Y[3]=4314.283;
*cos(w);
A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;
A[1+2*i][5]=-x[i]+x0;
L[1+2*i]=y[i]-y1[i];
}
/*-------------------解法方程--------------------*/
用矩阵表示为：
具体的流程图框：
是
否否
是
否
是
四、实验源程序：
#include<iostream>
#include<cmath>
using namespace std;
const int N=4;
const int n=6;
/*---------------矩阵相乘---------------------*/
Xs+=XX[0];
Ys+=XX[1];
Zs+=XX[2];
q+=XX[3];
w+=XX[4];
k+=XX[5];
n++;
}
/*----------------旋转矩阵R-----------------------*/
cout<<"迭代次数为:"<<n<<endl;
printf("\n像片外方位元素的解\n");
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;
else
q[i][j]=0;}
for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据
二、实验原理：
如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。
那么由共线条件方程可知：
其中 是只含三个独立参数ψ，ω，κ的九个方向余弦。
在方程中共有6个未知参数 所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，如下：
cout<<"航向顷角q:"<<q<<endl;
cout<<"旁向倾角w:"<<w<<endl;
cout<<"像片旋角k:"<<k<<endl;
cout<<"Xs "<<Xs<<endl;
cout<<"Ys "<<Ys<<endl;
cout<<"Zs "<<Zs<<endl;
cout<<endl;
cout<<"所有坐标已经内置,所以不显示了"<<endl; }
cout<<endl;
/*-----------------设定外方位元素初始值--------------*/
x0=0;y0=0;f=152;m=10000;
Xs=0;Ys=0;Zs=f*m/10000;
q=0;w=0;k=0;
XX[3]=1;
/*------------------迭代计算--------------------------*/
result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];
}
return;
}
/*---------------矩阵求逆---------------------*/
void inverse(double c[n][n])
{ int i,j,h,k;
double p;
double q[n][12];
用新的符号表示各偏导数后为
其中（x）、（y）是函数近似值， , 是外方位元素近似值的改正数，它们的系数为函数的偏导数。为了便于推导，令：
=
=
=
那么有
对于系数，其严密算法（以 , 为例）如下：
对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为：
为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数 , ，可列出每个点的误差方程式：
void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2)
{
int i,j,k;
for(i=0;i<i_1;i++)
for(j=0;j<j_2;j++)
{
result[i*j_2+j]=0.0;
for(k=0;k<j_12;k++)
江西理工大学应用科学学院
摄
影
测
量Байду номын сангаас
实
验
报
告
班级：测绘101班
姓名：邬志刚
学号：08090110125
一、实验目的：
掌握空间后方交会和前方交会的原理，根据所给控制点的地面物方坐标以及相应的像点在像平面坐标系中的坐标和同名像点在左右像片上的坐标，利用计算机编程语言实现空间后方交会、前方交会的过程，完成所给立体像对中两张像片各自的6个外方位元素的解求和所给立体像对上5对同名点对应的地面物方点的坐标计算。
Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs);
Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);