光束法区域网空中三角测量作业报告.docx

光束法区域网空中三角测量作业报告.docx

《光束法区域网空中三角测量作业报告.docx》由会员分享，可在线阅读，更多相关《光束法区域网空中三角测量作业报告.docx（23页珍藏版）》请在冰豆网上搜索。

光束法区域网空中三角测量作业报告

光束法区域网空中三角测量求像片外方位元素和加密点的地面坐标

一、原理：

首先给出的原始数据像素坐标转换成像平面直角坐标；然后利用后方交会原理并采用控制点的向平面直角坐标和地面坐标迭代求解每张像片的外方位元素；再通过前方交会求加密点的地面坐标；最后进行精度评定。

二、计算流程：

1、求出四个控制点及所有待求点的像点坐标（x1,y1）与（x2,y2）；

2、空间后方交会计算两像片的外方位元素：

根据编好的计算机程序读取原始文件中的数据控制点的地面坐标和相应的像点像素坐标（第一步已经转换成像平面直角坐标），对两张像片各自进行空间后方交会，计算各自的6个外方位元素Xs1,Ys1,Zs1,ψ1,w1,k1和Xs2,Ys2,Zs2,ψ2,w2,k2；

3、空间前方交会计算待定点的地面坐标：

用各片的外方位角元素计算左右片的方向余弦值，组成旋转矩阵R1与R2，逐点计算像点的像空间爱你辅助坐标（u1,v1,w1）及（u2,v2,w2）；根据外方位元素计算基线分量（Bu,Bv,Bw）;计算投影系数N1,N2；计算待定点在各自的像空间辅助坐标系中的坐标（U1,V1,W1）及（U2,V2,W2）；最后计算待定点的地面摄影测量坐标。

4、进行精度评定：

利用前方交会求出地面控制点的坐标，再与地面实测坐标求差，将两者差值视为真差值，由这些真差值计算出点位坐标精度。

三、基本公式：

1、像素坐标（I,J）转换成像平面直角坐标（x,y）：

x=（J-2136）*5.19/1000/1000;

y=（1424-I）*5.19/1000/1000;

2、后方交会：

共线条件方程：

外方位元素系数矩阵：

Zp=a3*（L_X-Xs）+b3*（L_Y-Ys）+c3*（L_Z-Zs）;

a11=1/Zp*（a1*f+a3*x）;

a12=1/Zp*（b1*f+b3*x）;

a13=1/Zp*（c1*f+c3*x）;

a14=y*Sin（w）-（x/f*（x*Cos（k）-y*Sin（k））+f*Cos（k））*Cos（w）;

a15=-f*Sin（k）-x/f*（x*Sin（k）+y*Cos（k））;

a16=y;

a21=1/Zp*（a2*f+a3*y）;

a22=1/Zp*（b2*f+b3*y）;

a23=1/Zp*（c2*f+c3*y）;

a24=-x*Sin（w）-（x/f*（x*Cos（k）-y*Sin（k））-f*Sin（k））*Cos（w）;

a25=-f*Cos（k）-y/f*（x*Sin（k）+y*Cos（k））;

a26=-x;

3、前方交会：

计算旋转矩阵：

a=Cos（ψ）*Cos（k）-Sin（w）*Sin（ψ）*Sin（k）;

a2=Cos（ψ）*Sin（k）-Sin（w）*Sin（ψ）*Cos（k）;

a3=-Sin（ψ）*Cos（w）;

b1=Cos（w）*Sin（k）;

b2=Cos（w）*Math.Cos（k）;

b3=Sin（w）;

c1==Sin（ψ）*.Cos（k）+Cos（ψ）*Sin（w）*Sin（k）;

c2=-Sin（ψ）*Sin（k）+Cos（ψ）*Sin（w）*Cos（k）;

c3=Cos（ψ）*Cos（w）;

基线分量：

Bu=Xs2-Xs1;

Bv=Ys2-Ys1;

Bw=Zs2-Zs1;

投影系数：

计算地面点坐标：

5、精度评定：

四、实例：

使用数码相机拍摄了两张航测相片Img261和Img262，量测了6对点的影像坐标，其中前4个点为地面控制点（GCP）,后两个点为加密点。

如下表所示。

相片参数：

大小为4272（列）*2848（行），像主点位于影像中心，像素大小为5.19micron，相机的焦距为24mm。

.

Img261的外方位元素初始值：

X=397510.760，Y=3445853.978,Z=1455.153，角元素为0；

Img262的外方位元素初始值：

X=397513.320，Y=3445979.811,Z=1453.685，角元素为0；

请采用光束法区域网空中三角测量方法，解求每张相片的外方位元素和加密点的地面坐标。

点名

Img261

Img262

地面坐标

x（pixel）

y（pixel）

x（pixel）

y（pixel）

X（m）

Y（m）

H（m）

1

3156.625

95.125

3064.875

903.625

397534.502

3446082.739

655.925

12

3206.857

983.316

3139.375

1792.375

397570.854

3445925.794

659.372

7

3705.375

462.625

3625.875

1251.625

397645.909

3446035.598

653.428

16

1660.823

1122.982

1598.875

1986.375

397303.238

3445852.908

677.949

3

3134.125

320.125

3047.875

1130.625

？

？

？

？

？

？

2

2188.625

624.875

2105.625

1491.125

？

？

？

？

？

？

原始数据文件：

生成结果文件：

程序源代码：

usingSystem.Collections.Generic;

usingSystem.Linq;

usingSystem.Text;

usingSystem.IO;

namespace摄影测量光束法编程作业

{

classProgram

{

staticvoidMain（string[]args）

{

//原始观测数据：

控制点和加密点的像素坐标（x1,y1,x2,y2），控制点的地面坐标（X,Y,Z）

FileStreamfs=newFileStream（@"控制点观测数据.txt",FileMode.Open,FileAccess.Read）;

FileStreamfs2=newFileStream（@"控制点观测数据.txt",FileMode.Open,FileAccess.Read）;

StreamReaderg_line=newStreamReader（fs）;

StreamReaderr=newStreamReader（fs2）;

intline=0;

while（g_line.ReadLine（）!

=null）

{

line++;

}

g_line.Close（）;

double[]L_X=newdouble[line];

double[]L_Y=newdouble[line];

double[]L_Z=newdouble[line];

double[]I1=newdouble[line];

double[]J1=newdouble[line];

double[]I2=newdouble[line];

double[]J2=newdouble[line];

double[]x1=newdouble[line];

double[]y1=newdouble[line];

double[]x2=newdouble[line];

double[]y2=newdouble[line];

intn=0;

for（n=0;n

{

stringasd=r.ReadLine（）;

string[]kon;

kon=asd.Split（newChar[]{',',',',','}）;

I1[n]=double.Parse（kon[0]）;

J1[n]=double.Parse（kon[1]）;

I2[n]=double.Parse（kon[2]）;

J2[n]=double.Parse（kon[3]）;

L_X[n]=double.Parse（kon[4]）;

L_Y[n]=double.Parse（kon[5]）;

L_Z[n]=double.Parse（kon[6]）;

}

//框标像素坐标转换成框标坐标（x1[line],y1[line],x2[line],y2[line]）

intj=0;

for（j=0;j

{

x1[j]=（J1[j]-2136）*5.19/1000/1000;

y1[j]=（1424-I1[j]）*5.19/1000/1000;

x2[j]=（J2[j]-2136）*5.19/1000/1000;

y2[j]=（1424-I2[j]）*5.19/1000/1000;

}

//相机焦距和相片外方位元素的初始值

doublef=0.024;

doubleXs1=397510.760;

doubleYs1=3445853.978;

doubleZs1=1455.153;

doublew1=0.0;

doubleψ1=0.0;

doublek1=0.0;

//计算旋转矩阵

MatrixR1=newMatrix（3,3）;

double[]dX1=newdouble[6];

doublexzhi=1.0/3600/180*Math.PI;

intd1=0;

do

{

d1++;

doublea1=R1[0,0]=Math.Cos（ψ1）*Math.Cos（k1）-Math.Sin（w1）*Math.Sin（ψ1）*Math.Sin（k1）;

doublea2=R1[0,1]=-Math.Cos（ψ1）*Math.Sin（k1）-Math.Sin（w1）*Math.Sin（ψ1）*Math.Cos（k1）;

doublea3=R1[0,2]=-Math.Sin（ψ1）*Math.Cos（w1）;

doubleb1=R1[1,0]=Math.Cos（w1）*Math.Sin（k1）;

doubleb2=R1[1,1]=Math.Cos（w1）*Math.Cos（k1）;

doubleb3=R1[1,2]=-Math.Sin（w1）;

doublec1=R1[2,0]=Math.Sin（ψ1）*Math.Cos（k1）+Math.Cos（ψ1）*Math.Sin（w1）*Math.Sin（k1）;

doublec2=R1[2,1]=-Math.Sin（ψ1）*Math.Sin（k1）+Math.Cos（ψ1）*Math.Sin（w1）*Math.Cos（k1）;

doublec3=R1[2,2]=Math.Cos（ψ1）*Math.Cos（w1）;

//计算像点坐标近似值，求出常数项lx1,ly1;

double[]jx1=newdouble[line];

double[]jy1=newdouble[line];

double[]lx1=newdouble[line];

double[]ly1=newdouble[line];

inti=0;

for（i=0;i

{

jx1[i]=-f*（a1*（L_X[i]-Xs1）+b1*（L_Y[i]-Ys1）+c1*（L_Z[i]-Zs1））/（a3*（L_X[i]-Xs1）+b3*（L_Y[i]-Ys1）+c3*（L_Z[i]-Zs1））;

jy1[i]=-f*（a2*（L_X[i]-Xs1）+b2*（L_Y[i]-Ys1）+c2*（L_Z[i]-Zs1））/（a3*（L_X[i]-Xs1）+b3*（L_Y[i]-Ys1）+c3*（L_Z[i]-Zs1））;

lx1[i]=x1[i]-jx1[i];

ly1[i]=y1[i]-jy1[i];

}

//像片的常数项为l1

double[]l1=newdouble[line*2];

for（intm2=0;m2

{

intm3=m2/2;

l1[m2]=lx1[m3];

l1[m2+1]=ly1[m3];

m2++;

}

//求出外方位元素系数矩阵HA1,HA2及其中的元素a11,a12,a13,a14,a15........

double[,]HA1=newdouble[line*2,6];

double[]a11=newdouble[line];

double[]a12=newdouble[line];

double[]a13=newdouble[line];

double[]a14=newdouble[line];

double[]a15=newdouble[line];

double[]a16=newdouble[line];

double[]a21=newdouble[line];

double[]a22=newdouble[line];

double[]a23=newdouble[line];

double[]a24=newdouble[line];

double[]a25=newdouble[line];

double[]a26=newdouble[line];

intn3=0;

for（n3=0;n3

{

intn1=n3/2;

doubleZp1=a3*（L_X[n1]-Xs1）+b3*（L_Y[n1]-Ys1）+c3*（L_Z[n1]-Zs1）;

a11[n1]=HA1[n3,0]=1/Zp1*（a1*f+a3*x1[n1]）;

a12[n1]=HA1[n3,1]=1/Zp1*（b1*f+b3*x1[n1]）;

a13[n1]=HA1[n3,2]=1/Zp1*（c1*f+c3*x1[n1]）;

a14[n1]=HA1[n3,3]=y1[n1]*Math.Sin（w1）-（x1[n1]/f*（x1[n1]*Math.Cos（k1）-y1[n1]*Math.Sin（k1））+f*Math.Cos（k1））*Math.Cos（w1）;

a15[n1]=HA1[n3,4]=-f*Math.Sin（k1）-x1[n1]/f*（x1[n1]*Math.Sin（k1）+y1[n1]*Math.Cos（k1））;

a16[n1]=HA1[n3,5]=y1[n1];

a21[n1]=HA1[n3+1,0]=1/Zp1*（a2*f+a3*y1[n1]）;

a22[n1]=HA1[n3+1,1]=1/Zp1*（b2*f+b3*y1[n1]）;

a23[n1]=HA1[n3+1,2]=1/Zp1*（c2*f+c3*y1[n1]）;

a24[n1]=HA1[n3+1,3]=-x1[n1]*Math.Sin（w1）-（x1[n1]/f*（x1[n1]*Math.Cos（k1）-y1[n1]*Math.Sin（k1））-f*Math.Sin（k1））*Math.Cos（w1）;

a25[n1]=HA1[n3+1,4]=-f*Math.Cos（k1）-y1[n1]/f*（x1[n1]*Math.Sin（k1）+y1[n1]*Math.Cos（k1））;

a26[n1]=HA1[n3+1,5]=-x1[n1];

n3++;

}

//求外方位元素的改正数

Matrixha1=newMatrix（HA1）;

MatrixL1=newMatrix（line*2,1,l1）;

MatrixHA1T=ha1.Transpose（）;

MatrixHA1THA1=HA1T.Multiply（ha1）;

HA1THA1.InvertGaussJordan（）;

Matrixz1=HA1THA1.Multiply（HA1T）;

Matrixz2=z1.Multiply（L1）;

dX1=z2;

Xs1=Xs1+dX1[0];

Ys1=Ys1+dX1[1];

Zs1=Zs1+dX1[2];

ψ1=ψ1+dX1[3];

w1=w1+dX1[4];

k1=k1+dX1[5];

}

while（Math.Abs（dX1[3]）>xzhi||Math.Abs（dX1[4]）>xzhi||Math.Abs（dX1[5]）>xzhi）;

//相机焦距和相片外方位元素的初始值

doubleXs2=397513.3200;

doubleYs2=3445979.811;

doubleZs2=1453.685;

doublew2=0.0;

doubleψ2=0.0;

doublek2=0.0;

//计算两张像片的外方位元素Xs1,Ys1,Zs1,ψ1,w1,k1，Xs2,Ys2,Zs2,ψ2,w2,k2

MatrixR2=newMatrix（3,3）;

double[]dX2=newdouble[6];

intd2=0;

do

{

d2++;

doubleaa1=R2[0,0]=Math.Cos（ψ2）*Math.Cos（k2）-Math.Sin（w2）*Math.Sin（ψ2）*Math.Sin（k2）;

doubleaa2=R2[0,1]=-Math.Cos（ψ2）*Math.Sin（k2）-Math.Sin（w2）*Math.Sin（ψ2）*Math.Cos（k2）;

doubleaa3=R2[0,2]=-Math.Sin（ψ2）*Math.Cos（w2）;

doublebb1=R2[1,0]=Math.Cos（w2）*Math.Sin（k2）;

doublebb2=R2[1,1]=Math.Cos（w2）*Math.Cos（k2）;

doublebb3=R2[1,2]=-Math.Sin（w2）;

doublecc1=R2[2,0]=Math.Sin（ψ2）*Math.Cos（k2）+Math.Cos（ψ2）*Math.Sin（w2）*Math.Sin（k2）;

doublecc2=R2[2,1]=-Math.Sin（ψ2）*Math.Sin（k2）+Math.Cos（ψ2）*Math.Sin（w2）*Math.Cos（k2）;

doublecc3=R2[2,2]=Math.Cos（ψ2）*Math.Cos（w2）;

//计算旋转矩阵

double[]jx2=newdouble[line];

double[]jy2=newdouble[line];

double[]lx2=newdouble[line];

double[]ly2=newdouble[line];

inti=0;

for（i=0;i

{

jx2[i]=-f*（aa1*（L_X[i]-Xs2）+bb1*（L_Y[i]-Ys2）+cc1*（L_Z[i]-Zs2））/（aa3*（L_X[i]-Xs2）+bb3*（L_Y[i]-Ys2）+cc3*（L_Z[i]-Zs2））;

jy2[i]=-f*（aa2*（L_X[i]-Xs2）+bb2*（L_Y[i]-Ys2）+cc2*（L_Z[i]-Zs2））/（aa3*（L_X[i]-Xs2）+bb3*（L_Y[i]-Ys2）+cc3*（L_Z[i]-Zs2））;

lx2[i]=x2[i]-jx2[i];

ly2[i]=y2[i]-jy2[i];

}

double[]l2=newdouble[line*2];

for（intm2=0;m2

{

intm3=m2/2;

l2[m2]=lx2[m3];

l2[m2+1]=ly2[m3];

m2++;

}

double[,]HA2=newdouble[line*2,6];

double[]aa11=newdouble[line];

double[]aa12=newdouble[line];

double[]aa13=newdouble[line];

double[]aa14=newdouble[line];

double[]aa15=newdouble[line];

double[]aa16=newdouble[line];

double[]aa21=newdouble[line];

double[]aa22=newdouble[line];

double[]