武汉大学测绘学院大地测量计算与实习报告

课程编号：1300158 课程性质：必修

大地测量计算与实习

实 习 报 告

学 院： 测绘学院

专业方向： 测绘工程A方向

实习地点： 武汉大学

班级组号： X 班第X小组

学生姓名： XXX

学号： 2011301610XXX

指导教师： 刘宗泉、苏新洲、丁士俊

向 东、史俊波 等

2013年 6月24日 至 2013年 7月12日

目 录.....................................................2

一 技术报告.............................................4

1.1.实习的目的和意义…………………………………………………………..4

1.2. 实习的任务................................................................4

二 二等精密水准测量外业观测与计算 .……………………………4

2.1. 测区概况.……………………………………………………4

2.2 已知高程.……………………………………………………4

2.3 作业依据.……………………………………………………5

2.4 水准测量详述.………………………………………………5

2.4.1踏勘、选点.………………………………………………6

2.4.2使用的仪器和仪器检验.…………………………………6

(1) i角检验..…………………………………………………6

(2) 一对水准尺零点差检验..……………………………… 6

2.4.3 外业数据采集与概算 ………………………………… 7

2.5 外业成果表格.…………………………………………… 7

三 大地测量计算课程设计.…………………………………………7

3.1 编程语言简介.………………………………………………7

3.2 编程任务、程序框图及基本数学模型……………………8

1）高斯投影坐标正反算.……………………………………8

2）实测斜距化算到高斯平面距离.…………………………10

3）大地主题正反算（高斯平均引数+白塞尔）.………………11

3.3 大地测量计算成果…………………………………………14

四 实习总结…………………………………………………………15

附录1 数字水准仪i角检验结果表…………………………………21

附录2 数字水准尺的零点差检验结果表……………………………22

附录 3 水准点点之记..………………………………………………23

附录4 自测二等水准外业观测高差与高程概算表.………………24

附录5 大地测量计算源程序.………………………………………25

附录6 国测二等水准外业观测高差与高程概算表.………………40

一、技术报告

1.1 实习的目的和意义

《大地测量计算与实习》课程，是学习《大地测量学基础》之后的后续课程，是一次重要的集中教学实习。

通过这次集中地教学实习，旨在巩固我们在课堂上学过的理论知识，锻炼和提高我们的仪器操作能力和动手能力，增强综合、灵活运用所学知识的能力；在我们自主编程计算的过程中，提高发现问题，分析问题和解决相关问题的能力；开拓我们的视野，把基本理论与实践紧密的结合在一起，对理论知识进一步消化和吸收，激化和培养我们的创新思维的能力，为今后走上工作岗位打下较为扎实的专业基础。

1.2 实习的任务

本次实习共有两项任务

1、二等精密水准测量外业观测与概算。

(1) 踏勘和选点：了解水准路线的布设和水准点位置的选择方法，根据安排的水准路线，各组独立完成本组的选点工作。

(2) 填写水准点之记：每人完成一个水准点“点之记”的绘制，且每个组的成员间不能同点。

(3) 仪器检验：每人完成数字水准仪i角测定成果和一对水准标尺零点差的测定结果各一份。

(4) 数字水准观测和记录：每人使用数字水准仪完成至少80站的观测和记簿，每小组完成珞珈山环或樱园环的往返闭合水准环线。

(5) 光学水准仪观测记录：每人使用光学精密水准仪完成至少10站的观测和记簿，每组构成闭合环。

(6) 外业观测数据检查与概算。

2、大地测量编程实习

(1) 高斯投影正反算的计算程序。

(2) 实测斜距化算至高斯投影平面边长的边长改正程序。

(3) 大地主题正反算的计算程序。

(4) 水准测量外业高差与概略高程表的编制。

二、二等精密水准测量外业观测与计算

2.1 测区概况

本次我们5班第2小组的测区是樱园环，这段水准线路线，地质为混凝土和柏油马路，部分测段穿过马路，来往车辆较多，教务部至行政楼高差起伏较大，其余路段较为平坦。

2.2 已知高程

2.3 作业依据

国家测绘局，国家一、二等水准测量规范2006-05-24 测绘出版社，2010.

2010仪器的技术指标

(1) 水准仪的i角限差为15″

(2) 标尺的零点不等差为0.10mm

(3) 测站观测限差

(4) 观测程序

(5) 闭合差

2.4 水准测量详述

2.4.1 踏勘、选点

2.4.2 使用的仪器和仪器检验

使用的仪器为DL200。进行了下面两项检验。

（1）i角检验

1.在较平坦的地方选定适当距离（约取40m）的两点A、B，并用尺垫确定不动。

2.置水准仪于A 点附近，如图：

此时测量A，B尺的读数a1,b1，并记录；

3.然后将水准仪置于B 点的右边，在线段AB 的延长线上，如图

读取此时A,B尺的读数a2,b2,并记录;

4.利用相关公式解算i角大小；

（2）一对水准尺零点差检验

在距离水准仪约20m—30m 的等距离处选择3 个点，使得各个点之间存在着约20cm 的高差。此项检验进行3 个测回，每一测回中，分别在3 个点上一次安置一对标尺，分别用仪器进行读数，测回间应变换仪器高。由于使用的是数字水准仪，因而设置重复测量次数为5 次，每测回每点连续观测4 次。水准尺零点差即两次读数中数之差。

检验的结果见附录。

2.4.3 外业数据采集与概算

水准测量路线图：

我观测的是行政楼到新闻与传播学院的往返测，与6月27日下午完成。往测采用手簿记录，返测为电子记录。大约有40站（加上后来补测的一些站数总数达到要求的80站）。往返测闭合差+2.0mm，具体成果见附录4。

2.5 外业成果表格

成果表格有：

1） 数字水准仪i角检验结果表(附录1)。

2） 数字水准尺的零点差检验结果表(附录2)。

3） 水准点点之记(附录3)。

4） 自测测量外业观测高差与高程概算成果表(附录4)。

三、大地测量计算课程设计

3.1 编程语言简介

本次大地测量编程实习，我采用的是C++语言，主要利用了MFC。C++语言是一种优秀的面向对象程序设计语言，它在C语言的基础上发展而来，但它比C语言更容易为人们学习和掌握。C++以其独特的语言机制在计算机科学的各个领域中得到了广泛的应用。面向对象的设计思想是在原来结构化程序设计方法基础上的一个质的飞跃，C++完美地体现了面向对象的各种特性。MFC(Microsoft Foundation Classes)，是一个微软公司提供的类库（class libraries），以C++类的形式封装了Windows的API，并且包含一个应用程序框架，以减少应用程序开发人员的工作量。其中包含的类包含大量Windows句柄封装类和很多Windows的内建控件和组件的封装类。

3.2 编程任务、程序框图及基本数学模型

本次大地测量实习主要编写三个程序，分别是高斯投影坐标正反算和实测斜距归算至高斯平面和大地主题正反算。

1） 高斯投影坐标正反算

1. 高斯投影坐标正算

高斯投影坐标正算公式

2. 高斯投影坐标反算

高斯投影坐标反算公式

主要代码见附录5.

2） 实测斜距归算至高斯平面边长

先利用高斯坐标计算大地方位角A，之后利用下面的公式求出实测斜距在椭球面上的大地线长度，

S=D*{[1-(h2-h1)/D*(h2-h1)/D]/[(1+h1/Ra)*(1+h2/Ra)]}

然后利用距离改化公式，

d=S*(1+ym*ym/2/Rm/Rm+(y2-y1)*(y2-y1)/24/Rm/Rm+ym*ym*ym*ym/24/Rm/Rm/Rm/Rm);

主要代码见附录5.

3） 大地测量主题正反算

1. 高斯平均引数正算

首先把勒让德级数在P1点展开改在大底线长度中低昂M处展开，使收敛级数减少，收敛快，精度高;其次，考虑到求中点M的复杂性，将M点用大地线两单点平均纬度及平均方位角相对应的m点来代替，并借助迭代计算，便可顺利的实现大地主题正算。

2. 高斯平均引数反算

由于已知L1、B1及L2、B2，这时经度差、纬度差及平均纬度均为已知，故可以正算公式很容易导出反算公式。

3. 白塞尔大地主题正算

（1） 计算起点的归化纬度

（2） 计算辅助函数值，解球面三角形

（3） 按公式计算相关系数A,B,C以及α,β

（4） 计算球面长度

（5） 计算纬度差改正数

（6） 计算终点大地坐标及大地方位角

4.白塞尔大地主题反算

（1） 辅助计算

（2） 用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。

计算下式,重复上述计算过程2.

（3） 计算大地线长度S

（4） 计算反方位角

主要代码见附录5.

3.3 大地测量计算成果

3.3.1 高斯投影正反算计算结果（采用第2组第2个数据）：

3.3.2 实测斜距化算至高斯投影平面边长计算结果（采用克拉索夫斯基椭球）：

D12 = 578.8686m。

3.3.3 大地主题解算计算结果（采用5班第2组第2个数据）：

3.3.4 国测数据外业高差与概略高程表计算结果

见附录6

四．实习总结

附录1 数字水准仪的i角检验结果表

i角的检校

仪器：DL-201/2007 方法：I1ABI2 观测者：王蔷

日期：2013-06-26 标尺： 记录者：黄式春

时间：09:00 呈像：清晰稳定 检查者：尹建鹏

附录二 数字水准仪的零点差检验结果表

一对水准尺零点差测定

标尺： 观测者：陆直

仪器： DL-201/2007 记录者：揭元翔

日期： 2013-06-29 检查者：尹建鹏

附录3 水准点点之记

二 等 水 准 点 之 记

武汉大学樱园 线 点名：B02

附录4 自测二等水准测量外业观测高差与高程概算成果表

ΣR=3.47km n=8 ==1.0mm ΣΔh=0.17mm

注：01：未名湖 02：行政楼 03：新闻与传播学院 04：城市设计学院 05：田园小观园 06：大自强超市 07：生科院 08：图书馆

附录5 大地测量计算源程序

1.自定义函数

//四舍五入保留n位小数

double round(double a,int n=0)

{

for (int i=0,times=1;i

times*=10;

a*=times;

double b=floor(a);

if(a-b>=0.5)

return (b+1)/times;

else

return b/times;

}

//度分秒表示法转换为度表示法

double dmsTod(double a)

{

double degree = 0, minute = 0, b = 0;

degree = floor(a);

a -= degree;

a *= 100;

minute = floor(a);

a -= minute;

a *= 100;

b = degree + minute / 60 + a / 3600;

return b;

}

//度表示法转换为度分秒表示法

double dTodms(double a)

{

double degree = 0, minute = 0, b = 0;

degree = floor(a);

a -= degree;

a *= 60;

minute = floor(a);

a -= minute;

a *= 60;

b = degree + minute / 100 + a / 10000;

return b;

}

//椭球参数的获取及相关量计算

int getpara(int sign,double &a,double &b,double &e_2,double &ee_2,double &a0,double &a2,double &a4,double &a6,double &a8)

{

a=0;e_2=0;ee_2=0; //椭球参数长半轴、第一偏心率平方、第二偏心率平方

switch(sign)

{

Case 0:{a = 6378245;b=6356863.0187730473;e_2 = 0.006693421622966;ee_2=0.006738525414683;}break;

case 1:{a = 6378140;b=6356755.2881575287;e_2 = 0.006694384999588;ee_2=0.006739501819473;}break;

case 2:{a = 6378137;b=6356752.3142 ;e_2 = 0.00669437999013;ee_2 = 0.00673949674227;}break;

case 3:{a = 6378137;b=6356752.3141 ;e_2 = 0.00669438002290;ee_2 = 0.00673949677548;}break;

default:return 1;

}

double m0=0,m2=0,m4=0,m6=0,m8=0;

m0 = a * (1 - e_2);

m2 = 3 * e_2 * m0 / 2;

m4 = 5 * e_2 * m2 / 4;

m6 = 7 * e_2 * m4 / 6;

m8 = 9 * e_2 * m6 / 8;

a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;

a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;

a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;

a6 = m6 / 32 + m8 / 16;

a8 = m8 / 128;

return 0;

}

2.高斯投影坐标正反算

//高斯投影坐标正算 sign表示椭球，m63表示分带方式(0代表6度带，其他代表3度带)

int GaussPro1(int ellipsoid,int m63,double B,double L,double &X,double &Y)

{

const double pi=4*atan(1);

double a=0,b=0,e_2=0,ee_2=0; //椭球参数

X=0;Y=0; //高斯坐标参数

double a0=0,a2=0,a4=0,a6=0,a8=0;

getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

B = dmsTod(B);

B = B * pi / 180; //B转换为弧度表示

L = dmsTod(L); //L转换为度表示

//计算中央子午线经度L0和带号No

int L0=0,No=0;

if (!m63)

{

No = static_cast(ceil(L/6));

L0 = 6 * No - 3;

}

else

{

if (L < 3)

L = L + 360;

No = static_cast(round(L/3));

L0 = 3 * No;

}

double l=0;

l = L - L0;

l = l * pi / 180; //l转换为弧度表示

X = a0 * B - a2 * sin(2*B) / 2 + a4 * sin(4*B) / 4 -a6 * sin(6*B) / 6 + a8 * sin(8*B) / 8;

double sinB=0,cosB=0,sinB_2=0,cosB_2=0,t=0,eta_2=0,N=0,m=0;

sinB=sin(B);

cosB=cos(B);

cosB_2=cosB*cosB;

sinB_2=sinB*sinB;

t = tan(B);

eta_2 = ee_2 * cosB_2;

N = a * sqrt(1.0/(1 - e_2 * sinB_2));

m = cosB * l;

X=X+N*t*((1.0/2+((5-t*t+9*eta_2+4*eta_2*eta_2)/24+(61-58*t*t+t*t*t*t)*m*m/720)*m*m)*m*m);

Y=N*((1+((1-t*t+eta_2)/6+(5-18*t*t+t*t*t*t+14*eta_2-58*eta_2*t*t)*m*m/120)*m*m)*m);

Y=No * 1000000 + 500000 + Y;

return 0;

}

//高斯投影坐标反算

int GaussPro2(int ellipsoid,int m63,double X,double Y,double &B,double &L)

{

const double pi=4*atan(1);

double a=0,b=0,e_2=0,ee_2=0; //椭球参数

B=0;L=0; //大地坐标参数

double a0=0,a2=0,a4=0,a6=0,a8=0;

getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

int No=0,L0=0;

double x=0,y=0;

No = static_cast(floor(Y / 1000000));

y = Y - No * 1000000 - 500000;

x = X;

if (!m63)

L0 = 6 * No - 3;

else

L0 = (3 * No)%360;

//通过迭代求底点纬度Bf

double eps=1.0E-15,Bf=0,Bf1=0;

Bf = X / a0;

do

{

Bf1 = Bf;

Bf= (X+a2 * sin(2 * Bf1) / 2 - a4 * sin(4 * Bf1) / 4 +a6 * sin(6 * Bf1) / 6 - a8 * sin(8 * Bf1)/ 8)/a0;

} while (fabs(Bf-Bf1)>eps);

double sinBf=0,cosBf=0,sinBf_2=0,cosBf_2=0;

double tf=0,Nf=0,Mf=0,etaf_2=0,Vf=0,yN=0;

sinBf = sin(Bf);

cosBf = cos(Bf);

cosBf_2 = cosBf * cosBf;

sinBf_2 = sinBf * sinBf;

tf = tan(Bf);

Nf = a / sqrt(1 - e_2 * sinBf_2);

Mf=a*(1-e_2)/pow((1-e_2*sinBf_2),1.5);

etaf_2 = ee_2 * cosBf_2;

Vf=sqrt(etaf_2);

yN = y / Nf;

double tf_2=0,Nf_3=0,l=0;

tf_2=tf*tf;

Nf_3 = Nf * Nf * Nf;

B = Bf - tf * y * y / 2 / Mf / Nf + tf * (5 + 3 * tf_2 + etaf_2 - 9 * etaf_2 * tf_2) * y * y * y * y/24/Mf/Nf_3 - tf * (61 + 90 * tf_2 + 45 * tf_2 * tf_2) * y * y * y / Nf_3* y * y * y / 720 / Mf /Nf/Nf;

l = y / Nf / cosBf - (1 + 2 * tf_2 + etaf_2) * y * y * y / 6 / Nf_3 / cosBf + (5 + 28 * tf_2 + 24 * tf_2 * tf_2 + 6 * etaf_2 + 8 * etaf_2 * tf_2) * y * y * y * y * y / 120 / Nf_3/Nf/Nf / cosBf;

B = B * 180 / pi;

l = l * 180 / pi;

L = L0 + l;

B = dTodms(B);

L = dTodms(L);

return 0;

}

void CGaussProDlg::OnButton1()

{

// TODO: Add your control notification handler code here

UpdateData(true);

if (m_63==-1)

{

MessageBox("请选择投影带！");

return;

}

if (m_ellipsoid==-1)

{

MessageBox("请选择椭球！");

return;

}

//m_B=37.37304420;m_L=112.23191622;

GaussPro1(m_ellipsoid,m_63,m_B,m_L,m_X,m_Y);

VarR8Round(m_X,4,&m_X);

VarR8Round(m_Y,4,&m_Y);

UpdateData(false);

}

void CGaussProDlg::OnButton2()

{

// TODO: Add your control notification handler code here

UpdateData(true);

if (m_63==-1)

{

MessageBox("请选择投影带！");

return;

}

if (m_ellipsoid==-1)

{

MessageBox("请选择椭球！");

return;

}

GaussPro2(m_ellipsoid,m_63,m_X,m_Y,m_B,m_L);

VarR8Round(m_B,10,&m_B);

VarR8Round(m_L,10,&m_L);

UpdateData(false);

}

3.实测斜距化算至高斯平面边长

void CDistanceDlg::OnButton1()

{

// TODO: Add your control notification handler code here

UpdateData(true);

const double pi=4*atan(1);

double x1=m_x1,y1=m_y1,B1=m_B1,H1=m_H1+50,x2=m_x2,y2=m_y2,B2=m_B2,H2=m_H2+50,D=m_D;

double a=0,b=0,e_2=0,ee_2=0,a0=0,a2=0,a4=0,a6=0,a8=0,A=0,N=0,RA=0;

getpara(0,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

//根据高斯坐标反算大地坐标

double B11=0,B22=0,L11=0,L22=0;

GaussPro2(0,0,x1,y1+19000000,B11,L11);

GaussPro2(0,0,x2,y2+19000000,B22,L22);

//利用高斯平均引数法反算点1到点2的大地方位角A

GeoSolution2(0,0,B11,L11,B22,L22,A,a0,a2);

A=dmsTod(A)*pi/180;

B1=dmsTod(B1)*pi/180;

B2=dmsTod(B2)*pi/180;

double Bm=(B1+B2)/2,Nm=0,Vm=0,Rm=0;

Vm=sqrt(1+ee_2*cos(Bm)*cos(Bm));

Nm=a/sqrt(1-e_2*sin(Bm)*sin(Bm));

Rm=Nm/Vm;

N=a/sqrt(1-e_2*sin(B1)*sin(B1));

RA=N/(1+ee_2*cos(B1)*cos(B1)*cos(A)*cos(A));

double ym=(y1+y2)/2-500000,S=0,d=0;

S=D*sqrt((1-(H2-H1)*(H2-H1)/D/D)/(1+H1/RA)/(1+H2/RA))+D*D*D/24/RA/RA;

d=S*(1+ym*ym/2/Rm/Rm+(y2-y1)*(y2-y1)/24/Rm/Rm+ym*ym*ym*ym/24/Rm/Rm/Rm/Rm);

VarR8Round(d,4,&d);

m_d=d;

UpdateData(false);

}

4.大地测量主题正反算

//大地主题正算 ellopsoid代表椭球体,solution是解算方法(0是高斯平均引数,其他是白塞尔)

int GeoSolution1(int ellipsoid,int solution,double B1,double L1,double A12,double S,double &B2,double &L2,double &A21)

{

const double pi=4*atan(1);

double a=0,b=0,e_2=0,ee_2=0,a0=0,a2=0,a4=0,a6=0,a8=0;

getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

B1 = dmsTod(B1)*pi/180;

L1 = dmsTod(L1)*pi/180;

A12 = dmsTod(A12)*pi/180;

if (!solution)

{

double Bm=B1,Am=A12,Vm=0,Nm=0,tm=0,etam_2=0,deltaB=0,deltaL=0,deltaA=0,deltaB0=0,deltaL0=0,deltaA0=0,sinAm=0,cosAm=0;

B2=B1;A21=A12;

do

{

deltaB0=deltaB;

deltaL0=deltaL;

deltaA0=deltaA;

Bm=(B1+B2)/2;

Am=(A12+A21)/2;

sinAm=sin(Am);

cosAm=cos(Am);

Vm=sqrt(1+ee_2*cos(Bm)*cos(Bm));

Nm=a/sqrt(1-e_2*sin(Bm)*sin(Bm));

tm=tan(Bm);

etam_2=ee_2*cos(Bm)*cos(Bm);

deltaB=Vm*Vm/Nm*S*cosAm*(1+S*S/24/Nm/Nm*(sinAm*sinAm*(2+3*tm*tm*(1+etam_2))+3*etam_2*cosAm*cosAm*(-1+tm*tm-etam_2-4*tm*tm*etam_2)));

deltaL=S/Nm/cos(Bm)*sinAm*(1+S*S/24/Nm/Nm*(sinAm*sinAm*tm*tm-cosAm*cosAm*(1+etam_2-9*tm*tm*etam_2+etam_2*etam_2)));

deltaA=S/Nm*sinAm*tm*(1+S*S/24/Nm/Nm*(cosAm*cosAm*(2+7*etam_2+9*tm*tm*etam_2+5*etam_2*etam_2)+sinAm*sinAm*(2+tm*tm+2*etam_2)));

B2=B1+deltaB;

L2=L1+deltaL;

A21=A12+deltaA;

} while (fabs(deltaB-deltaB0)>1.0E-15||fabs(deltaL-deltaL0)>1.0E-15||fabs(deltaA-deltaA0)>1.0E-15);

if(A12>pi)

A21-=pi;

else

A21+=pi;

}

else

{

double sinB1=sin(B1),cosB1=cos(B1),sinA12=sin(A12),cosA12=cos(A12),W1=0,sinu1=0,cosu1=0;

W1 = sqrt(1 - (e_2 * sinB1 * sinB1));

sinu1 = sinB1 * sqrt(1 -e_2) / W1;

cosu1 = cosB1 / W1;

double sinA0=0,cosA0_2=0,k_2=0,A=0,B=0,C=0,alpha=0,beta=0;

sinA0 = cosu1 * sinA12;

cosA0_2 = 1 - sinA0 * sinA0;

k_2 = ee_2 * cosA0_2;

A = (1 + k_2 / 4 - 3 * k_2 * k_2 / 64 + 5 * k_2 * k_2 * k_2 / 256) * b;

B = (k_2 / 8 - k_2 * k_2 / 32 + 15 * k_2 * k_2 * k_2 / 1024) * b;

C = (k_2 * k_2 / 128 - 3 * k_2 * k_2 * k_2 / 512) * b;

alpha = (e_2 / 2 + e_2 * e_2 / 8 + e_2 * e_2 * e_2 / 16) - (e_2 * e_2 + e_2 * e_2 * e_2) / 16 * cosA0_2 + 3 * e_2 * e_2 * e_2 / 128 * cosA0_2 * cosA0_2;

beta = (e_2 * e_2 + e_2 * e_2 * e_2) / 32 * cosA0_2 - e_2 * e_2 * e_2 * cosA0_2 * cosA0_2 / 64;

double cotsigma1=0,sin2sigma1=0,cos2sigma1=0;

cotsigma1 = cosu1 * cosA12 / sinu1;

sin2sigma1 = 2 * cotsigma1 / (cotsigma1 * cotsigma1 + 1);

cos2sigma1 = (cotsigma1 * cotsigma1 - 1) / (cotsigma1 * cotsigma1 + 1);

double sigma0=0,sin2sigma0=0,cos2sigma0=0;

sigma0 = (S - (B + C * cos2sigma1) * sin2sigma1) / A;

sin2sigma0=sin(2*sigma0);

cos2sigma0=cos(2*sigma0);

double sin2=0,cos2=0,sigma=0,cossigma=0,sinsigma=0,delta=0;

sin2 = sin2sigma1 * cos2sigma0 + cos2sigma1 * sin2sigma0;

cos2 = cos2sigma1 * cos2sigma0 - sin2sigma1 * sin2sigma0;

sigma = sigma0 + (B + 5 * C * cos2) * sin2 / A;

cossigma = cos(sigma);

sinsigma = sin(sigma);

delta = (alpha * sigma + beta * (sin2 - sin2sigma1)) * sinA0;

double sinu2=0;

sinu2 = sinu1 * cossigma + cosu1 * cosA12 * sinsigma;

B2 =atan(sinu2 / sqrt((1 - e_2) * (1 - sinu2 * sinu2)));

double tanlambda=0,lambda=0;

tanlambda = sinA12 * sinsigma / (cosu1 * cossigma - sinu1 * sinsigma * cosA12);

lambda = atan(tanlambda);

if (tanlambda > 0 && sinA12 < 0)

lambda += pi;

if (tanlambda < 0 && sinA12 > 0)

lambda -= pi;

L2 = L1 + lambda - delta;

double tanA21=0;

tanA21 = cosu1 * sinA12 / (cosu1 * cossigma * cosA12 - sinu1 * sinsigma);

A21 = atan(tanA21);

if (sinA12 * tanA21 > 0)

A21 += pi;

if (sinA12 > 0 && tanA21 < 0)

A21 += 2 * pi;

}

B2 = dTodms(B2*180/pi);

L2 = dTodms(L2*180/pi);

A21 = dTodms(A21*180/pi);

return 0;

}

//大地主题反算

int GeoSolution2(int ellipsoid,int solution,double B1,double L1,double B2,double L2,double &A12,double &A21,double &S)

{

const double pi=4*atan(1);

double a=0,b=0,e_2=0,ee_2=0,a0=0,a2=0,a4=0,a6=0,a8=0;

getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

B1 = dmsTod(B1)*pi/180;

B2 = dmsTod(B2)*pi/180;

L1 = dmsTod(L1)*pi/180;

L2 = dmsTod(L2)*pi/180;

if (!solution)

{

double deltaL=L2-L1,deltaB=B2-B1,Bm=(B1+B2)/2,cosBm=cos(Bm);

double Vm=0,Nm=0,tm=0,etam_2=0;

Vm=sqrt(1+ee_2*cosBm*cosBm);

Nm=a/sqrt(1-e_2*sin(Bm)*sin(Bm));

tm=tan(Bm);

etam_2=ee_2*cosBm*cosBm;

double r01=0,r21=0,r03=0,S10=0,S12=0,S30=0;

r01=Nm*cosBm;

r21=Nm*cosBm/24*(1+etam_2-9*etam_2*tm*tm+etam_2*etam_2);

r03=-Nm/24*cosBm/Vm/Vm/Vm/Vm*cosBm*cosBm*tm*tm;

S10=Nm/Vm/Vm;

S12=-Nm/24/Vm/Vm*cosBm*cosBm*(2+3*tm*tm+2*etam_2);

S30=Nm/8/Vm/Vm/Vm/Vm/Vm/Vm*(etam_2-tm*tm*etam_2);

double SsinAm=0,ScosAm=0;

SsinAm=r01*deltaL+r21*deltaB*deltaB*deltaL+r03*deltaL*deltaL*deltaL;

ScosAm=S10*deltaB+S12*deltaB*deltaL*deltaL+S30*deltaB*deltaB*deltaB;

double t01=0,t21=0,t03=0,deltaA=0;

t01=tm*cosBm;

t21=cosBm*tm*(2+7*etam_2+9*tm*tm*etam_2+5*etam_2*etam_2)/24/Vm/Vm/Vm/Vm;

t03=cosBm*cosBm*cosBm*tm*(2+tm*tm+2*etam_2)/24;

deltaA=t01*deltaL+t21*deltaB*deltaB*deltaL+t03*deltaL*deltaL*deltaL;

double b=B2-B1,l=L2-L1,c=0,T=0,Am=0;

c=fabs(ScosAm/SsinAm);

if (fabs(b)>=fabs(l))

T=atan(fabs(SsinAm/ScosAm));

else

T=pi/4+atan(fabs((1-c)/(1+c)));

if (b>0&&l>=0)

Am=T;

else if(b<0&&l>=0)

Am=pi-T;

else if(b<=0&&l<0)

Am=pi+T;

else if(b>0&&l<0)

Am=2*pi-T;

else

Am=pi/2;

S=SsinAm/sin(Am);

A12=Am-deltaA/2;

A21=Am+deltaA/2+pi;

if(A21>=2*pi)

A21-=2*pi;

}

else

{

double sinB1=sin(B1),cosB1=cos(B1),sinB2=sin(B2),cosB2=cos(B2),W1=0,W2=0,sinu1=0,sinu2=0,cosu1=0,cosu2=0;

W1 = sqrt(1 - e_2 * sinB1 * sinB1);

W2 = sqrt(1 - e_2 * sinB2 * sinB2);

sinu1 = sinB1 * sqrt(1 - e_2) / W1;

sinu2 = sinB2 * sqrt(1 - e_2) / W2;

cosu1 = cosB1 / W1;

cosu2 = cosB2 / W2;

double L=0,a1=0,a2=0,b1=0,b2=0;

L = L2 - L1;

a1 = sinu1 * sinu2;

a2 = cosu1 * cosu2;

b1 = cosu1 * sinu2;

b2 = sinu1 * cosu2;

double lambda=L,delta=0,delta1=0,sinlambda=0,coslambda=0,p=0,q=0,sinA12=0,cosA12=0,sigma=0,sinsigma=0,cossigma=0,

sinA0=0,cosA0_2=0,k_2=0,alpha=0,beta=0,betab=0,x=0;

do

{

delta1 = delta;

sinlambda=sin(lambda);

coslambda=cos(lambda);

p=cosu2*sinlambda;

q=b1-b2*coslambda;

A12=atan2(p,q);

if(A12<0)

A12+=2*pi;

sinA12=sin(A12);

cosA12=cos(A12);

sinsigma = p * sinA12 + q * cosA12;

cossigma = a1 + a2 * coslambda;

sigma = atan(sinsigma / cossigma);

if (sigma < 0)

sigma += pi;

sinA0 = cosu1 * sinA12;

cosA0_2 = 1 - sinA0 * sinA0;

k_2 = ee_2 * cosA0_2;

alpha = (e_2 / 2 + e_2 * e_2 / 8 + e_2 * e_2 * e_2 / 16) - (e_2 * e_2 + e_2 * e_2 * e_2) / 16 * cosA0_2 + 3 * e_2 * e_2 * e_2 / 128 * cosA0_2 * cosA0_2;

beta = (e_2 * e_2 + e_2 * e_2 * e_2) / 32 * cosA0_2 - e_2 * e_2 * e_2 * cosA0_2 * cosA0_2 / 64;

betab = 2 * beta;

x = 2 * a1 - cosA0_2 * cossigma;

delta=(alpha*sigma-betab*x*sinsigma)*sinA0;

lambda = L + delta;

} while (fabs(delta-delta1)>1.0E-15);

double A=0,B=0,C=0,BB=0,CC=0,y=0;

sinlambda = sin(lambda);

coslambda = cos(lambda);

k_2 = ee_2 * cosA0_2;

A = (1 + k_2 / 4 - 3 * k_2 * k_2 / 64 + 5 * k_2 * k_2 * k_2 / 256) * b;

B = (k_2 / 8 - k_2 * k_2 / 32 + 15 * k_2 * k_2 * k_2 / 1024) * b;

C = (k_2 * k_2 / 128 - 3 * k_2 * k_2 * k_2 / 512) * b;

BB = 2 * B / cosA0_2;

CC = 2 * C / cosA0_2 / cosA0_2;

y = (cosA0_2 * cosA0_2 - 2 * x * x) * cossigma;

S = A * sigma + (BB * x + CC * y) * sinsigma;

double sinA21=0,cosA21=0;

sinA21 = -cosu1 * sinlambda;

cosA21 = -b1 * coslambda + b2;

A21 = atan2(sinA21, cosA21);

if (A21 < 0)

A21 += 2*pi;

}

A12 = dTodms(A12 * 180 / pi);

A21 = dTodms(A21 * 180 / pi);

return 0;

}

void CBesselDlg::OnButton1()

{

// TODO: Add your control notification handler code here

UpdateData(true);

if (m_ellipsoid==-1)

{

MessageBox("请选择椭球！");

return;

}

if (m_solution==-1)

{

MessageBox("请选择解算方法");

return;

}

GeoSolution1(m_ellipsoid,m_solution,m_B1,m_L1,m_A12,m_S12,m_B2,m_L2,m_A21);

VarR8Round(m_B2,9,&m_B2);

VarR8Round(m_L2,9,&m_L2);

VarR8Round(m_A21,9,&m_A21);

UpdateData(false);

}

void CBesselDlg::OnButton2()

{

// TODO: Add your control notification handler code here

UpdateData(true);

if (m_ellipsoid==-1)

{

MessageBox("请选择椭球！");

return;

}

if (m_solution==-1)

{

MessageBox("请选择解算方法");

return;

}

GeoSolution2(m_ellipsoid,m_solution,m_B1,m_L1,m_B2,m_L2,m_A12,m_A21,m_S12);

VarR8Round(m_A12,9,&m_A12);

VarR8Round(m_A21,9,&m_A21);

VarR8Round(m_S12,3,&m_S12);

UpdateData(false);

}

附录6 国测精密水准测量外业观测高差与高程概算成果表

大地测量计算与实习成绩评定表

本文来源：https://www.2haoxitong.net/k/doc/87d49f31f7ec4afe04a1df88.html