空间后方交会

空间后方交会实习报告

遥感信息工程学院沈翔第一部分:空间后方交会的基本原理

空间后方交会的基本原理：已知影像范围内一定数量的控制点的地面坐标(或地面摄影坐标)及相对应的影像坐标,求的该影像的外方位元素的方法。主要分为：共线条件方程法、角锥体法、直接线性变换法。本程序采用共线条件方程法。

第二部分:空间后方交会的计算机程序框图

其主要公式如下：

共线方程：

xfa1(XXs)b1(YYs)c1(ZZs)x0a3(XXs)b3(YYs)c3(ZZs)a2(XXs)b2(YYs)c2(ZZs)yfy0a3(XXs)b3(YYs)c3(ZZs)

经线性化后，得到的误差方程式为：

xxxxxxXsYsZsXsYsZs yyyyyyvx(y)yXsYsZsXsYsZsvx(x)x

组成法方程：

ATPAXATPL

解方程得：

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

//函数模型

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

for(i=0;i

{

for(j=0;j

cout

}

cout

for(i=0;i

{

cout

for(j=0;j

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分代表的弧度

//开始循环

{

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

GetR(phi,omega,kappa,R);

cout

for(i=0;i

{

cout

for(j=0;j

cout

}

return 0;

}

m0=sqrt(v/(2*N-6));

for(j=0;j

{

mi[j]=m0*sqrt(*(ata_+7*j));

}

//结果输出

cout

= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

cout

GetR(phi,omega,kappa,R);

cout

for(i=0;i

{

cout

for(j=0;j

cout

}

return 1; PRS021沈翔 [1**********]6 }

2004.11.1

空间后方交会实习报告

遥感信息工程学院沈翔第一部分:空间后方交会的基本原理

第二部分:空间后方交会的计算机程序框图

其主要公式如下：

共线方程：

xfa1(XXs)b1(YYs)c1(ZZs)x0a3(XXs)b3(YYs)c3(ZZs)a2(XXs)b2(YYs)c2(ZZs)yfy0a3(XXs)b3(YYs)c3(ZZs)

经线性化后，得到的误差方程式为：

组成法方程：

ATPAXATPL

解方程得：

X(ATA)1ATL

第三部分：试验结果

附录:源程序代码

影像外方位元素并确定其精度

其已知变量如下：

内方位元素：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

//函数模型

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

for(i=0;i

{

for(j=0;j

cout

}

cout

for(i=0;i

{

cout

for(j=0;j

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分代表的弧度

//开始循环

{

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

GetR(phi,omega,kappa,R);

cout

for(i=0;i

{

cout

for(j=0;j

cout

}

return 0;

}

m0=sqrt(v/(2*N-6));

for(j=0;j

{

mi[j]=m0*sqrt(*(ata_+7*j));

}

//结果输出

cout

= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

cout

GetR(phi,omega,kappa,R);

cout

for(i=0;i

{

cout

for(j=0;j

cout

}

return 1; PRS021沈翔 [1**********]6 }

2004.11.1

相关文章