单张影像空间后方交会程序

-----单张影像空间后方交会程序

实习时间 2013.5.20-2013.5.24 学生班级 10测绘（1）班 学生姓名 房新明 学生学号 ********** 所在院系 矿业工程学院 指导老师 张会战 邵亚琴

一．实习目的

1．深入理解单张影像空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。

2．利用Visual C++编写一个完整的单张影像空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并进行评定精度。

3．通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的动手能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。

4、掌握空间后方交会的定义和实现算法

(1) 定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。

(2) 算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。

5、了解空间后方交会的基本过程

(1) 获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。

(2) 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。

(3) 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。

(4) 计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。 (5) 逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。

(6) 逐点计算误差方程式的系数和常数项，组成误差方程式。 (7) 计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。

(8) 解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。

(9) 检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.1′，当3个改正数均小于0.1′时，迭代结束。否则用新的近似值重复(4)～（8）步骤的计算，直到满足要求为止。

二、实习要求

1、认真复习单张影像空间后方交会的原理及有关内容；

2、编写程序，计算左片号23和右片号24的外方位元素并进行精度评定； 3、验证数据部分外方位元素和旋转矩阵R已经给出，请验证自己编写的程序是否正确；

三、实习环境

1、硬件环境：windows操作系统

2、软件环境：VC++

四、实习原理

1、以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素 Xs，Ys ，zs ，∮ ，ω，κ共线条件方程如下： x-x0=-f*[a1(X-Xs)+b1(Y-Ys)+c1(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] y-y0=-f*[a2(X-Xs)+b2(Y-Ys)+c2(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] 其中：

x,y为像点的像平面坐标； Xs,Ys ,f为影像的外方位元素； Xs，Ys ,Zs为摄站点的物方空间坐标； X,Y,Z为物方点的物方空间坐标； 旋转矩阵为R ;

由于此共线条件方程是非线性方程，先对其进行线性化，像点观测值一般视为等权，即P=I；

矩阵形式：V=AX-L，P=I；

通过间接平差，为提高精度，增加多余观测方程，根据最小二乘平差原理，可计算出外方位元素的改正数。经过迭代计算，每次迭代用未知数的近似值与上次迭代计算的改正数之和作为新的近似值，重复计算，求出新的改正数，这样反复趋近，直到改正数小于某个限值为止。 2、精度评定

Mi=m0*√Qii 其中m0=±√【VV】/（2N-6）

3、公式推导

Fx(XS,YS,ZS,,,)Fx0()0Fx0XS0(XSXS)Fx0YS(YSYS0)Fx0ZS0(ZSZS)Fx0()0Fx000(0)Fx(XS,YS0,ZS,0,0,0)Fy(XS,YS,ZS,,,)Fy0()0Fy0XS0(XSXS)Fy0YS(YSYS0)Fy0ZS0(ZSZS)Fy0()0Fy000(0)Fy(XS,YS0,ZS,0,0,0)

同理

FyFy0XSdXSFy0YSdYSFy0ZSdZSFy0dFy0dFy0dFy0所以

Fx0XSdXSFx0YSdYSFx0ZSdZSFx0dFx0dFx0dFx00Fy0XS

dXSFy0YSdYSFy0ZSdZSFy0dFy0dFy0dFy00

1、基本原理

共线条件线性方程式为：

fxx2xydXsdZs(f)ddyd(xx计)0ffZZ2fdYsydZsxyd(fy)dxd(yy)0计ffZZ误差方程式为：

2fxxxy dXsdZs(f)ddyd(xx计)vxffZZ 2fdYsydZsxyd(fy)dxd(yy)vy计ffZZ

计算改正值： dXs,dYs,dZs,d,d,d计算改后的外方位元素：

XsXsdXs,YsYsdYs,ZsZsdZsd,d,d六

、

程

序

过

程

框

：

图

开始 获取已知数据：比例尺1/M 地面点坐标： 像点坐标：x(4),y(4) 内方位元素初始值：Xs0,Ys0,Zs0,Phi0,Ome0,Kap0 定义旋转矩阵R 计算x(i),y(i)的近视值Jinsi(i),Jinsi(i) 组成矩阵A(8,6)和L(8,1) 计算内方位元素的改正值： dXs, dYs, dZs, dp,do,dk 否 判断内方位元素改正值是否收敛 是 输出Xs, Ys, Zs, φ, ω, κ 结束

七．实习数据

1、模拟像片一对：左片23号 右片24号 2、相片比例尺：1/30000 3、航摄机主距：f=150mm 4、每张相片有四个控制点 5、各片像点坐标及地面坐标 片点像点坐标（mm） 号 号 1 x -91.596 -94.230 95.207 96.797 y X 地面坐标（m） Y 137500.0 142500.00 137500.00 142500.00 137500.00 142500.00 137500.00 142500.00 Z 11.00 36.00 42.00 56.00 90.00 31.00 7.00 5.00 -74.859 100000.0 81.446 100000.0 23 3 7 9 4 6 -75.512 106000.0 83.077 106000.0 -102.695 -79.618 103000.0 -99.904 86.890 88.904 81.754 103000.0 24 10 12 -77.540 109000.0 76.257 109000.0 八．验证数据

已知四对点的影像坐标和地面坐标及航摄仪内方位元素

f=153.42mm，X0=Y0=0

影像坐标 地面坐标 1 2 3 4 x(mm) -86.15 -53.40 -14.78 10.46 Y(mm) -68.99 82.21 -76.63 64.43 X（mm） Y（mm） Z（mm） 36589.41 25273.32 2195.17 37631.08 31324.51 728.69 39100.97 24934.98 2386.50 40426.54 30319.81 757.31 3、结算结果

Xs=39795.45 0.99771 0.06753 0.00399 Ys=27476.46 R= -0.06753 0.99772 -0.00221 Zs=7572.69 -0.00412 0.00184 0.99990

九．实现程序

输入文件形式如下：

C++源程序如下：

#include

#include #include #include #include using namespace std; const int n=6;

void inverse (double c[n][n]);

templatevoid transpose (T1*mat1,T2*mat2,int a,int b); templatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);

templatevoid input (T*mat,int a,int b);

templatevoid output(T*mat,char*s,int a,int b); int main() {

ofstream outFile; cout.precision(5);

double x0=0.0, y0=0.0; double fk=0.15324; //内方位元素 double m=39689; //估算比例尺

double B[4][5]={0.0},R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];

input (B,4,5); //从文件中读取控制点的影像坐标和地面坐标，存入数组B double Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0; double X,Y,Z,L[8][1],A[8][6]; //确定未知数的出始值 for(int i=0;i<4;i++) {Xs=Xs+B[i][2]; Ys=Ys+B[i][3]; Zs=Zs+B[i][4]; }

Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk; int f=0;

do//迭代计算 {f++;

//组成旋转矩阵

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R[0][2]=-sin(Q)*cos(W); R[1][0]=cos(W)*sin(K); R[1][1]=cos(W)*cos(K); R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K); R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K); R[2][2]=cos(Q)*cos(W); //计算系数阵和常数项

for(int i=0,k=0,j=0;i<=3;i++,k++,j++) { X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs); Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs); Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs); L[j][0]=B[i][0]-(x0-fk*X/Z); L[j+1][0]=B[i][1]-(y0-fk*Y/Z);

j++;

A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z; A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z; A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;

A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk+fk*cos(K))*cos(W);

A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk; A[k][5]=B[i][1]-y0; A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z; A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z; A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;

A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][1]-y0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin(K))*cos(W);

A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk; A[k+1][5]=-(B[i][0]-x0); k++; }

transpose(A,AT,6,8); multi(AT,A,ATA,6,8,6); inverse(ATA);

multi(AT,L,ATL,6,8,1); multi(ATA,ATL,XG,6,6,1);

Xs=Xs+XG[0][0]; Ys=Ys+XG[1][0]; Zs=Zs+XG[2][0]; Q=Q+XG[3][0]; W=W+XG[4][0]; K=K+XG[5][0];

}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265.0);

cout<<\"迭代次数为：\"<double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6]; multi(A,XG,AXG,8,6,1);

for( i=0;i<8;i++) //计算改正数 V[i][0]=AXG[i][0]-L[i][0]; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV[0][0]/2; for(i=0;i<6;i++) for(int j=0;j<6;j++) D[i][j]=m0*ATA[i][j];

//屏幕输出误差方程系数阵、常数项、改正数 output(A,\"误差方程系数阵A为：\ output(L,\"常数项L为：\ output(XG,\"改正数为：\

outFile.open(\"aim.txt\ //打开并添加aim.txt文件 outFile.precision(10);

//以文件的形式输出像片外方位元素、旋转矩阵、方差阵 outFile<<\"一、像片的外方位元素为：\"<outFile<<outFile<\"旁向倾角为：\"<outFile<outFile<outFile<outFile.close(); return 0; }

templatevoid transpose(T1*mat1,T2*mat2,int a,int b) { int i,j;

for(i=0;itemplatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) { int i,j,k;

for(i=0;ireturn;

}

template void input (T*mat,int a,int b) { ifstream inFile;

inFile.open(\"控制点坐标.txt\"); while(!inFile.eof()) {for (int i=0;i>mat[i][j]; }

inFile.close(); return; }

templatevoid output(T*mat,char*s,int a,int b) { cout<cout<return; }

void inverse(double c[n][n]) { int i,j,h,k; double p;

double q[n][12];

for(i=0;ifor(j=n;j<12;j++) {if(i+6==j) q[i][j]=1; else

q[i][j]=0;}

for(h=k=0;kp=q[k][h]/q[i][h]; for(j=0;j<12;j++) { q[i][j]*=p;

q[i][j]-=q[k][j]; } }

for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据 for(i=k-1;i>=0;i--) {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;ifor(i=0;i程序的结果输出如下：（包括文本输出结果和荧屏输出中间数据）

十．

十一．实习心得

通过一周的实习，使我感觉得到了自己在摄影测量知识本身的欠缺和和实现语言方面知识的欠缺。

当任务布置下来时，就不知所措。听别人说，MATLAB和VB++简单容易实现这个程序。先从图书馆借了本MATLAB来学学，几天只学了个皮毛，对编

出这个程序来说远远不够，我们学的是C语言目，前来说编出这个程序的可能性不太大，也只能从网上找个别人已经实现过的程序。

先对上学期所学的该知识点进行了系统的复习，从整体上把握了单相空间后方交会原理及解算过程。也理解了将控制点的地面点坐标视为真值，把相应的像点坐标视为观测值，加入相应的改正数，最后列出每个点的误差方程，然后用矩阵表示，最后构成总误差方程，从而求出外方位元素的改正数。

在编程里也学习了好多以前没见过的函数，mult是矩阵相乘的函数、invers_矩阵求逆函数、transpose矩阵转换函数。这对以后的学习是大有裨益的。

此外，通过此次实习，更使我确定在这个假期了学习VB和MATLAB的信心。