C# 双像空间前方交会

C# 双像空间前方交会

1 简介

首先，以txt文件的形式读取立体像对的外方位元素和同名像点的像空间直角坐标系坐标。然后，依据双像前方交会公式计算出左右相片投影系数、同名像点的像空间辅助坐标系坐标及对应的地面摄影测量坐标系坐标。最后，将中间及最终计算结果用txt文件保存。

2 算法设计与流程图

以txt文件形式存放相片的内、外方位元素及同名像点坐标，数据格式如图1所示。本程序首先是读取内、外方位元素及同名像点坐标；然后，分别依据公式（1）、（2）和（3）计算同名像点的空间辅助坐标系坐标、左右相片的投影系数N1、N2和同名像点对应的地面摄影测量坐标系坐标；最后，将中间结果、最终结果输出到txt文件。算法流程如下图所示。

3 主要函数和变量

3.1 ExternalElement类

存放当张像片的六个外方位元素和焦距，包含Xs, Ys, Zs, o, p, q, f 7个变量，变量均为double类型。

3.2 ImagePoint类

存放像点像点坐标相关属性，计算空间辅助坐标系坐标u，v，w。7个变量分别是id、x、y、u、v、w、element，其中id为string类型，为像点点号；element为ExternalElement类，放置该像片的外方位元素和焦距。

3.3 TerrainPoint类

地面点坐标，包括计算投影系数和地面点坐标X、Y、Z。包括变量id、N1、N2、X、Y、Z

4 部分代码

publicclassExternalElement//外方位元素{publicdouble Xs, Ys, Zs, o, p, q, f;//六个外方位元素publicExternalElement(double Xs,double Ys,double Zs,double o,double p,double q,double f){this.Xs = Xs;this.Ys = Ys;this.Zs = Zs;this.o = o;this.p = p;this.q = q;this.f = f;}publicExternalElement(){
            Xs =-6911.42788;
            Ys =4181.15686;
            Zs =157.77319;
            o =0.34831;
            p =-0.30914;
            q =0.08136;
            f =-165.37034;}}publicclassImagePoint//像点{publicstring id;publicdouble x, y;publicdouble u, v, w;//空间辅助坐标系坐标publicExternalElement element;//外方位元素publicImagePoint(string id,double x,double y,ExternalElement element){this.id = id;this.x = x;this.y = y;this.element = element;double f = element.f;double o = element.o;double p = element.p;double q = element.q;double pI = Math.PI;double a1 = Math.Cos(o /360*2* pI)* Math.Cos(q /360*2* pI)- Math.Sin(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double a2 =-Math.Cos(o /360*2* pI)* Math.Sin(q /360*2* pI)- Math.Sin(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double a3 =-Math.Sin(o /360*2* pI)* Math.Cos(p /360*2* pI);double b1 = Math.Cos(p /360*2* pI)* Math.Sin(q /360*2* pI);double b2 = Math.Cos(p /360*2* pI)* Math.Cos(q /360*2* pI);double b3 =-Math.Sin(p /360*2* pI);double c1= Math.Sin(o /360*2* pI)* Math.Cos(q /360*2* pI)+ Math.Cos(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double c2 =-Math.Sin(p /360*2* pI)* Math.Cos(q /360*2* pI)+ Math.Cos(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double c3 = Math.Cos(o /360*2* pI)* Math.Cos(q /360*2* pI);this.u = a1 * x + a2 * y + a3 * f;this.v = b1 * x + b2 * y + b3 * f;this.w = c1 * x + c2 * y + c3 * f;}publicImagePoint(string id,double x,double y){this.x = x;this.y = y;this.element =newExternalElement();double f = element.f;double o = element.o;double p = element.p;double q = element.q;double pI = Math.PI;double a1 = Math.Cos(o /360*2* pI)* Math.Cos(q /360*2* pI)- Math.Sin(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double a2 =-Math.Cos(o /360*2* pI)* Math.Sin(q /360*2* pI)- Math.Sin(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double a3 =-Math.Sin(o /360*2* pI)* Math.Cos(p /360*2* pI);double b1 = Math.Cos(p /360*2* pI)* Math.Sin(q /360*2* pI);double b2 = Math.Cos(p /360*2* pI)* Math.Cos(q /360*2* pI);double b3 =-Math.Sin(p /360*2* pI);double c1 = Math.Sin(o /360*2* pI)* Math.Cos(q /360*2* pI)+ Math.Cos(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double c2 =-Math.Sin(p /360*2* pI)* Math.Cos(q /360*2* pI)+ Math.Cos(o /360*2* pI)* Math.Sin(p /360*2* pI)* Math.Sin(q /360*2* pI);double c3 = Math.Cos(o /360*2* pI)* Math.Cos(q /360*2* pI);this.u = a1 * x + a2 * y - a3 * f;this.v = b1 * x + b2 * y - b3 * f;this.w = c1 * x + c2 * y - c3 * f;}}publicclassTerrainPoint//地面点坐标{publicstring id;publicdouble N1, N2;//投影系数publicdouble X, Y, Z;//地面点坐标publicTerrainPoint(ImagePoint p1,ImagePoint p2){double Bu = p2.element.Xs - p1.element.Xs;double Bv = p2.element.Ys - p1.element.Ys;double Bw = p2.element.Zs - p1.element.Zs;this.id = p1.id;this.N1 =(Bu * p2.w - Bw * p2.u)/(p1.u * p2.w - p2.u * p1.w);this.N2 =(Bu * p1.w - Bw * p1.u)/(p1.u * p2.w - p2.u * p1.w);this.X = p1.element.Xs + N1 * p1.u;this.Y =0.5*((p1.element.Ys + N1 * p1.v)+(p2.element.Ys + N1 * p2.v));this.Z = p1.element.Zs + N1 * p1.w;}}

5 程序运行界面

6 完整程序

下载地址