空间后方交会实习报告
遥感信息工程学院 沈翔 第一部分:空间后方交会的基本原理
空间后方交会的基本原理: 已知影像范围内一定数量的控制点的地面坐标(或地面摄影坐标)及相对应的影像坐标,求的该影像的外方位元素的方法。主要分为:共线条件方程法、角锥体法、直接线性变换法。本程序采用共线条件方程法。
第二部分:空间后方交会的计算机程序框图
其主要公式如下:
共线方程:
xfa1(XXs)b1(YYs)c1(ZZs)x0a3(XXs)b3(YYs)c3(ZZs)a2(XXs)b2(YYs)c2(ZZs)yfy0a3(XXs)b3(YYs)c3(ZZs)
经线性化后,得到的误差方程式为:
xxxxxxXsYsZsXsYsZs yyyyyyvx(y)yXsYsZsXsYsZsvx(x)x
组成法方程:
ATPAXATPL
解方程得:
X(ATA)1ATL
第三部分:试验结果
附录:源程序代码
/*本程序为空间后方交会计算程序(在windows 2003及Visual C++6.0平台中编译通过) 原理:由控制点(至少三个)的地面坐标及相应的影像坐标利用共线方程及最小二乘法来求得
影像外方位元素并确定其精度
其已知变量如下:
内方位元素:x0,y0,f
影象比例尺:m
控制点影象坐标:photo[N][2]
控制点地面坐标:earth[N][3]
其主要中间变量如下:
计算的像点坐标: x[N],y[N]
A阵元素: a[2][6]
L阵元素: l[2]
待求变量如下:
外方位元素:xs,ys,zs,phi,omega,kappa
R阵元素:R[3][3]
单位权中误差:m0
精度:mi[6]*/
const int N=4;
#include
#include
//函数模型
void GetR(double phi,double omega,double kappa,double R[3][3]);//求R阵
void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double
x0,double y0,double f,double xs,double ys,double zs);//计算控制点像点
坐标近似值
void GetA(double a[2][6],double R[3][3],double earth[N][3],double photo[N][2],double
x0,double y0,double f,double phi,double omega,double kappa,double
zs,int i);//求A阵
void GetL(double l[2],double photo[N][2],double x[N],double y[N],int i);//求L阵
void MatrixT(double a[2][6],double at[6][2],int m,int n);//矩阵转秩
double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,int bcol);//矩阵相乘 double MatrixDet(double x[],int n); //求矩阵行列式
double *MatrixInver(double *x,int n);//求逆矩阵
////////////////////////////////////////////////////////////////////////////////////////////////////
void GetR(double phi,double omega,double kappa,double R[3][3])//求R阵
{
R[0][0]=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);
R[0][1]=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);
R[0][2]=-sin(phi)*cos(omega);
R[1][0]=cos(omega)*sin(kappa);
R[1][1]=cos(omega)*cos(kappa);
R[1][2]=-sin(omega);
R[2][0]=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);
R[2][1]=-sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa);
R[2][2]=cos(phi)*cos(omega);
}
void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double
x0,double y0,double f,double xs,double ys,double zs)//计算控制点像点
坐标近似值
{
int i;
for(i=0;i
{
x[i]=x0-f*(R[0][0]*(earth[i][0]-xs)+R[1][0]*(earth[i][1]-ys)+R[2][0]*(earth[i][2]-zs))/(R[0][ 2]*(earth[i][0]-xs)+R[1][2]*(earth[i][1]-ys)+R[2][2]*(earth[i][2]-zs));
y[i]=y0-f*(R[0][1]*(earth[i][0]-xs)+R[1][1]*(earth[i][1]-ys)+R[2][1]*(earth[i][2]-zs))/(R[0][
2]*(earth[i][0]-xs)+R[1][2]*(earth[i][1]-ys)+R[2][2]*(earth[i][2]-zs));
}
}
void GetA(double a[2][6],double R[3][3],double earth[N][3],double photo[N][2],double
x0,double y0,double f,double phi,double omega,double kappa,double
zs,int i)//求A阵
{
a[0][0]=(R[0][0]*f+R[0][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][1]=(R[1][0]*f+R[1][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][2]=(R[2][0]*f+R[2][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][3]=(photo[i][1]-y0)*sin(omega)-((photo[i][0]-x0)*((photo[i][0]-x0)*cos(kappa)-(ph
oto[i][1]-y0)*sin(kappa))/f+f*cos(kappa))*cos(omega);
a[0][4]=-f*sin(kappa)-(photo[i][0]-x0)*((photo[i][0]-x0)*sin(kappa)+(photo[i][1]-y0)*c
os(kappa))/f;
a[0][5]=photo[i][1]-y0;
a[1][0]=(R[0][1]*f+R[0][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][1]=(R[1][1]*f+R[1][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][2]=(R[2][1]*f+R[2][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][3]=-(photo[i][0]-x0)*sin(omega)-((photo[i][1]-y0)*((photo[i][0]-x0)*cos(kappa)-(p
hoto[i][1]-y0)*sin(kappa))/f-f*sin(kappa))*cos(omega);
a[1][4]=-f*cos(kappa)-(photo[i][1]-y0)*((photo[i][0]-x0)*sin(kappa)+(photo[i][1]-y0)*c
os(kappa))/f;
a[1][5]=-(photo[i][0]-x0);
}
void GetL(double l[2],double photo[N][2],double x[N],double y[N],int i)//求L阵
{
l[0]=photo[i][0]-x[i];
l[1]=photo[i][1]-y[i];
}
void MatrixT(double a[2][6],double at[6][2],int m,int n)//矩阵转秩
{
int i,j;
for(i=0;i
for(j=0;j
{
at[j][i]=a[i][j];
}
}
double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,int bcol)//矩阵相乘 {
int i,j,p;
double *c;
c=new double[arow*bcol];
for(i=0;i
for(j=0;j
{
*(c+i*bcol+j)=0;
for(p=0;p
*(c+i*bcol+j)+=*(a+i*acol+p)* *(b+p*bcol+j);
}
return c;
}
double MatrixDet(double x[],int n) //求矩阵行列式
{
int i,j,p,ss;
double det=0.0;
if(n= =2)
{
det=x[0]*x[3]-x[1]*x[2];
}
else
{
double *z;
z=new double[(n-1)*(n-1)];
for(p=0;p
{
for(i=1;i
for(j=0;j
{
if(j
z[(i-1)*(n-1)+j]=x[i*n+j];
if(j>p)
z[(i-1)*(n-1)+j-1]=x[i*n+j];
}
if(p%2==0)
ss=1;
else
ss=-1;
det+=x[p]*ss*MatrixDet(z,n-1);
}
delete []z;
}
return det;
}
double *MatrixInver(double *x,int n)//求逆矩阵
{
int i,j,p,q,ss;
double det,det_1;
double *xx=NULL;
xx=new double[n*n];
det=MatrixDet(x,n);
for(i=0;i
for(j=0;j
{
double *x_1;//余子式
x_1=new double[(n-1)*(n-1)];
for(p=0;p
for(q=0;q
{
if((p
if((pj)) x_1[p*(n-1)+q-1]=x[p*n+q];
if((p>i)&&(q
if((p>i)&&(q>j)) x_1[(p-1)*(n-1)+q-1]=x[p*n+q];
det_1=MatrixDet(x_1,n-1);
if((i+j)%2==0) ss=1;
else ss=-1;
xx[j*n+i]=ss*det_1/det;
}
}
return xx;
}
int main()
{
//定义变量
int i,j=0,n=0;//i,j:循环次数//n:收敛次数
double photo[N][2]={-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,
0.01046,0.06443};//像点坐标
double earth[N][3]={{36589.41,25273.32,2195.17},{37631.08,31324.51,728.69}
,{39100.97,24934.98,2386.50},{40426.54,30319.81,757.31}};//地面点坐标
double x0=0,y0=0,f=0.15324;//像片内方位元素
double m=50000;//比例尺
double xs,ys,zs,phi,omega,kappa;//像片外方位元素
double delta;//角元素限值
double R[3][3];//R阵元素
double x[N],y[N];//计算的像点坐标
double a[2][6],at[6][2];//A阵及A阵的转帙阵
double l[2];//L阵
double *ata=new double[36],*ata_=new double[36],*atl=new double[6];//指向中间结果的
指针
double *X;//改正数的指针
double v;//v的值
double m0,mi[6];//单位权中误差及未知数的精度
//输出初始值
cout
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cout
cout
cout
for(i=0;i
{
for(j=0;j
cout
cout
}
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
//确定未知数的初始值
for(i=0,xs=0,ys=0;i
{
xs+=earth[i][0];
ys+=earth[i][1];
}
xs=xs/N;
ys=ys/N;
zs=m*f;
phi=0.0;
omega=0.0;
kappa=0.0;
delta=0.1*60/206265;//0.1分代表的弧度
//开始循环
do
{
v=0;
GetR(phi,omega,kappa,R);//求R阵
Earth2Photo(earth,R,x,y,x0,y0,f,xs,ys,zs);//计算控制点坐标的近似值
for(i=0;i
{
GetA(a,R,earth,photo,x0,y0,f,phi,omega,kappa,zs,i);//求得A阵
GetL(l,photo,x,y,i);//求L阵
for(j=0;j
{
v+=l[j]*l[j];
}
MatrixT(a,at,2,6);//a矩阵转帙为at矩阵
if(i==0)
{
ata=MatrixMutiple(at[0],6,2,a[0],2,6);//at阵与a阵相乘得ata阵 atl=MatrixMutiple(at[0],6,2,l,2,1);
}
else
{
for(j=0;j
*(ata+j)+=*(MatrixMutiple(at[0],6,2,a[0],2,6)+j);
for(j=0;j
*(atl+j)+=*(MatrixMutiple(at[0],6,2,l,2,1)+j);
}
}
ata_=MatrixInver(ata,6);//ata阵求逆得ata_阵
X=MatrixMutiple(ata_,6,6,atl,6,1);//ata_at阵与l阵相乘得X阵
xs+=*(X+0);//近似值与改正数相加得新的近似值
ys+=*(X+1);
zs+=*(X+2);
phi+=*(X+3);
omega+=*(X+4);
kappa+=*(X+5);
n++;
}while((*(X+3)>=delta||*(X+4)>=delta||*(X+5)>=delta)&&(n5)
{
cout
x0 cout
cout
cout
cout
GetR(phi,omega,kappa,R);
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
return 0;
}
m0=sqrt(v/(2*N-6));
for(j=0;j
{
mi[j]=m0*sqrt(*(ata_+7*j));
}
//结果输出
cout
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cout
cout
cout
cout
cout
cout
cout
cout
cout
cout
GetR(phi,omega,kappa,R);
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
return 1; PRS021沈翔 [1**********]6 }
2004.11.1
空间后方交会实习报告
遥感信息工程学院 沈翔 第一部分:空间后方交会的基本原理
空间后方交会的基本原理: 已知影像范围内一定数量的控制点的地面坐标(或地面摄影坐标)及相对应的影像坐标,求的该影像的外方位元素的方法。主要分为:共线条件方程法、角锥体法、直接线性变换法。本程序采用共线条件方程法。
第二部分:空间后方交会的计算机程序框图
其主要公式如下:
共线方程:
xfa1(XXs)b1(YYs)c1(ZZs)x0a3(XXs)b3(YYs)c3(ZZs)a2(XXs)b2(YYs)c2(ZZs)yfy0a3(XXs)b3(YYs)c3(ZZs)
经线性化后,得到的误差方程式为:
xxxxxxXsYsZsXsYsZs yyyyyyvx(y)yXsYsZsXsYsZsvx(x)x
组成法方程:
ATPAXATPL
解方程得:
X(ATA)1ATL
第三部分:试验结果
附录:源程序代码
/*本程序为空间后方交会计算程序(在windows 2003及Visual C++6.0平台中编译通过) 原理:由控制点(至少三个)的地面坐标及相应的影像坐标利用共线方程及最小二乘法来求得
影像外方位元素并确定其精度
其已知变量如下:
内方位元素:x0,y0,f
影象比例尺:m
控制点影象坐标:photo[N][2]
控制点地面坐标:earth[N][3]
其主要中间变量如下:
计算的像点坐标: x[N],y[N]
A阵元素: a[2][6]
L阵元素: l[2]
待求变量如下:
外方位元素:xs,ys,zs,phi,omega,kappa
R阵元素:R[3][3]
单位权中误差:m0
精度:mi[6]*/
const int N=4;
#include
#include
//函数模型
void GetR(double phi,double omega,double kappa,double R[3][3]);//求R阵
void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double
x0,double y0,double f,double xs,double ys,double zs);//计算控制点像点
坐标近似值
void GetA(double a[2][6],double R[3][3],double earth[N][3],double photo[N][2],double
x0,double y0,double f,double phi,double omega,double kappa,double
zs,int i);//求A阵
void GetL(double l[2],double photo[N][2],double x[N],double y[N],int i);//求L阵
void MatrixT(double a[2][6],double at[6][2],int m,int n);//矩阵转秩
double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,int bcol);//矩阵相乘 double MatrixDet(double x[],int n); //求矩阵行列式
double *MatrixInver(double *x,int n);//求逆矩阵
////////////////////////////////////////////////////////////////////////////////////////////////////
void GetR(double phi,double omega,double kappa,double R[3][3])//求R阵
{
R[0][0]=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);
R[0][1]=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);
R[0][2]=-sin(phi)*cos(omega);
R[1][0]=cos(omega)*sin(kappa);
R[1][1]=cos(omega)*cos(kappa);
R[1][2]=-sin(omega);
R[2][0]=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);
R[2][1]=-sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa);
R[2][2]=cos(phi)*cos(omega);
}
void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double
x0,double y0,double f,double xs,double ys,double zs)//计算控制点像点
坐标近似值
{
int i;
for(i=0;i
{
x[i]=x0-f*(R[0][0]*(earth[i][0]-xs)+R[1][0]*(earth[i][1]-ys)+R[2][0]*(earth[i][2]-zs))/(R[0][ 2]*(earth[i][0]-xs)+R[1][2]*(earth[i][1]-ys)+R[2][2]*(earth[i][2]-zs));
y[i]=y0-f*(R[0][1]*(earth[i][0]-xs)+R[1][1]*(earth[i][1]-ys)+R[2][1]*(earth[i][2]-zs))/(R[0][
2]*(earth[i][0]-xs)+R[1][2]*(earth[i][1]-ys)+R[2][2]*(earth[i][2]-zs));
}
}
void GetA(double a[2][6],double R[3][3],double earth[N][3],double photo[N][2],double
x0,double y0,double f,double phi,double omega,double kappa,double
zs,int i)//求A阵
{
a[0][0]=(R[0][0]*f+R[0][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][1]=(R[1][0]*f+R[1][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][2]=(R[2][0]*f+R[2][2]*(photo[i][0]-x0))/(earth[i][2]-zs);
a[0][3]=(photo[i][1]-y0)*sin(omega)-((photo[i][0]-x0)*((photo[i][0]-x0)*cos(kappa)-(ph
oto[i][1]-y0)*sin(kappa))/f+f*cos(kappa))*cos(omega);
a[0][4]=-f*sin(kappa)-(photo[i][0]-x0)*((photo[i][0]-x0)*sin(kappa)+(photo[i][1]-y0)*c
os(kappa))/f;
a[0][5]=photo[i][1]-y0;
a[1][0]=(R[0][1]*f+R[0][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][1]=(R[1][1]*f+R[1][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][2]=(R[2][1]*f+R[2][2]*(photo[i][1]-y0))/(earth[i][2]-zs);
a[1][3]=-(photo[i][0]-x0)*sin(omega)-((photo[i][1]-y0)*((photo[i][0]-x0)*cos(kappa)-(p
hoto[i][1]-y0)*sin(kappa))/f-f*sin(kappa))*cos(omega);
a[1][4]=-f*cos(kappa)-(photo[i][1]-y0)*((photo[i][0]-x0)*sin(kappa)+(photo[i][1]-y0)*c
os(kappa))/f;
a[1][5]=-(photo[i][0]-x0);
}
void GetL(double l[2],double photo[N][2],double x[N],double y[N],int i)//求L阵
{
l[0]=photo[i][0]-x[i];
l[1]=photo[i][1]-y[i];
}
void MatrixT(double a[2][6],double at[6][2],int m,int n)//矩阵转秩
{
int i,j;
for(i=0;i
for(j=0;j
{
at[j][i]=a[i][j];
}
}
double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,int bcol)//矩阵相乘 {
int i,j,p;
double *c;
c=new double[arow*bcol];
for(i=0;i
for(j=0;j
{
*(c+i*bcol+j)=0;
for(p=0;p
*(c+i*bcol+j)+=*(a+i*acol+p)* *(b+p*bcol+j);
}
return c;
}
double MatrixDet(double x[],int n) //求矩阵行列式
{
int i,j,p,ss;
double det=0.0;
if(n= =2)
{
det=x[0]*x[3]-x[1]*x[2];
}
else
{
double *z;
z=new double[(n-1)*(n-1)];
for(p=0;p
{
for(i=1;i
for(j=0;j
{
if(j
z[(i-1)*(n-1)+j]=x[i*n+j];
if(j>p)
z[(i-1)*(n-1)+j-1]=x[i*n+j];
}
if(p%2==0)
ss=1;
else
ss=-1;
det+=x[p]*ss*MatrixDet(z,n-1);
}
delete []z;
}
return det;
}
double *MatrixInver(double *x,int n)//求逆矩阵
{
int i,j,p,q,ss;
double det,det_1;
double *xx=NULL;
xx=new double[n*n];
det=MatrixDet(x,n);
for(i=0;i
for(j=0;j
{
double *x_1;//余子式
x_1=new double[(n-1)*(n-1)];
for(p=0;p
for(q=0;q
{
if((p
if((pj)) x_1[p*(n-1)+q-1]=x[p*n+q];
if((p>i)&&(q
if((p>i)&&(q>j)) x_1[(p-1)*(n-1)+q-1]=x[p*n+q];
det_1=MatrixDet(x_1,n-1);
if((i+j)%2==0) ss=1;
else ss=-1;
xx[j*n+i]=ss*det_1/det;
}
}
return xx;
}
int main()
{
//定义变量
int i,j=0,n=0;//i,j:循环次数//n:收敛次数
double photo[N][2]={-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,
0.01046,0.06443};//像点坐标
double earth[N][3]={{36589.41,25273.32,2195.17},{37631.08,31324.51,728.69}
,{39100.97,24934.98,2386.50},{40426.54,30319.81,757.31}};//地面点坐标
double x0=0,y0=0,f=0.15324;//像片内方位元素
double m=50000;//比例尺
double xs,ys,zs,phi,omega,kappa;//像片外方位元素
double delta;//角元素限值
double R[3][3];//R阵元素
double x[N],y[N];//计算的像点坐标
double a[2][6],at[6][2];//A阵及A阵的转帙阵
double l[2];//L阵
double *ata=new double[36],*ata_=new double[36],*atl=new double[6];//指向中间结果的
指针
double *X;//改正数的指针
double v;//v的值
double m0,mi[6];//单位权中误差及未知数的精度
//输出初始值
cout
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cout
cout
cout
for(i=0;i
{
for(j=0;j
cout
cout
}
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
//确定未知数的初始值
for(i=0,xs=0,ys=0;i
{
xs+=earth[i][0];
ys+=earth[i][1];
}
xs=xs/N;
ys=ys/N;
zs=m*f;
phi=0.0;
omega=0.0;
kappa=0.0;
delta=0.1*60/206265;//0.1分代表的弧度
//开始循环
do
{
v=0;
GetR(phi,omega,kappa,R);//求R阵
Earth2Photo(earth,R,x,y,x0,y0,f,xs,ys,zs);//计算控制点坐标的近似值
for(i=0;i
{
GetA(a,R,earth,photo,x0,y0,f,phi,omega,kappa,zs,i);//求得A阵
GetL(l,photo,x,y,i);//求L阵
for(j=0;j
{
v+=l[j]*l[j];
}
MatrixT(a,at,2,6);//a矩阵转帙为at矩阵
if(i==0)
{
ata=MatrixMutiple(at[0],6,2,a[0],2,6);//at阵与a阵相乘得ata阵 atl=MatrixMutiple(at[0],6,2,l,2,1);
}
else
{
for(j=0;j
*(ata+j)+=*(MatrixMutiple(at[0],6,2,a[0],2,6)+j);
for(j=0;j
*(atl+j)+=*(MatrixMutiple(at[0],6,2,l,2,1)+j);
}
}
ata_=MatrixInver(ata,6);//ata阵求逆得ata_阵
X=MatrixMutiple(ata_,6,6,atl,6,1);//ata_at阵与l阵相乘得X阵
xs+=*(X+0);//近似值与改正数相加得新的近似值
ys+=*(X+1);
zs+=*(X+2);
phi+=*(X+3);
omega+=*(X+4);
kappa+=*(X+5);
n++;
}while((*(X+3)>=delta||*(X+4)>=delta||*(X+5)>=delta)&&(n5)
{
cout
x0 cout
cout
cout
cout
GetR(phi,omega,kappa,R);
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
return 0;
}
m0=sqrt(v/(2*N-6));
for(j=0;j
{
mi[j]=m0*sqrt(*(ata_+7*j));
}
//结果输出
cout
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cout
cout
cout
cout
cout
cout
cout
cout
cout
cout
GetR(phi,omega,kappa,R);
cout
for(i=0;i
{
cout
for(j=0;j
cout
cout
}
return 1; PRS021沈翔 [1**********]6 }
2004.11.1