小秦的《大地测量学》笔记1

        俺是有普普通通有点点憨憨的小秦同学。一整个星期都莫得更新惹！辣么为什么捏？
       因为小秦在整理各种专业课的笔记，偶尔敲敲代码，运动运动，摸摸鱼，时间就这样悄悄过去辣◑﹏◐

所以，这次不来点干货可就真过意不去辣！下面和小秦一起开启《大地测量学》的学习之路叭！

PS：小秦简单介绍后正篇直接从教材的第4章开始

在开始之前先简单介绍下我们用的教材哦。

武大的《大地测量学基础》(第二版)

1.大地测量学学了什么东东？

学习知识的第一步当然是了解我们在学什么那么大地测量学到底要研究虾米捏？~(￣▽￣)~*

我们首先康康书上的高大上定义：

大地测量学(英语 Geodesy,德语 Geodasie)是在一定的时间-空间参考系统中,测量和描绘地球及其他行星体的一门学科。

大地测量学是地球科学中的一个分支,而且是发展最活跃、最具有重要地位的一个分它的最基本的任务是测量和描绘地球并监测其变化,为人类活动提供关于地球等行星体的间信息。因此,从本质上讲,它是一门地球信息学科,既是基础学科,又是应用学科。

这乍一看可把小秦高兴坏了，是不是精通了这门课小秦幼儿园想成为大科学家的理想是不是就能实现了！

理想和现实果然还是差了亿点点

第四章开始关于地球椭球及其数学投影变换的基本理论中繁复的空间解析几何公式推导和拗口难懂的专有名词可真是让人头秃。

不过小秦认为4.1-4.5也就是下面的正篇不过是那些大科学家们试图认识地球这个大胖小子长啥样而想出的各种方法，嘿嘿。

下面小秦将用Python代码的形式定义和计算地球椭球的各个数学要素。

2.确定一个地球椭球——地球椭球的基本参数

我们小学一年级就学过，地球不是个球是个赤道略宽两极略扁的且凹凸不平的球体。那我们想要用数学的方法研究地球就学要将它用规则的数学模型表示。

在控制测量中，用来代表地球的椭球叫做地球椭球(简称椭球)，是地球的数学代表。 

 

而具有一定几何参数，定位及定向的用以代表某一地区的大地水准面的椭球叫做参考椭球。

地球椭球的长半轴为a，短半轴为b，椭球是绕短轴旋转而成的几何体。

包含旋转轴的平面与椭球面相截所得椭圆叫做子午圈，垂直于旋转轴的平面与椭球面相截的圆叫平行圈，赤道是最大的平行圈，南北极是最小的平行圈。

任何椭球都由子午椭圆的五个基本元素决定：

椭球长半轴a、短半轴b、扁率α、第一和第二偏心率。

 为满足不同地区的不同建设及研究需要，全球建立了参数各异的椭球参数，来最大程度的使建立的椭球贴近自己所处区域的地球表面。

同时椭球的各个参数间也存在着它们间的数学关系。

我们可以用程序定义保存椭球各要素，形成自己的大地测量计算模块。

下面为第一部分:椭球实例及椭球要素

import math
import numpy as np
from jzh import dfm2hdz as d2h,hdz2dfm as h2d

 首先使导入我们要运用的函数库，math库为简单数学计算库，内置三角函数等常用算法。Numpy库为功能强大的数据处理库，小秦的代码将用到Numpy的矩阵计算功能，而jzh(角转换)库是小秦自己编写的测量中常用的各种角度变换功能，此处d2h()为度分秒转弧度制，h2d()则相异。

class Ellipsoid():
    # 引用类时声明椭球名称，如 Ellipsoid('CGCS2000')
    def __init__(self,name,a=0,b=0):
        self.name = name
        if name == 'CGCS2000' or name == 'cgcs2000':
            # 长半轴a
            self.a = 6378137
            # 短半轴b
            self.b = 6356752.3141
            # 极曲率半径c
            self.c = self.a**2/self.b
            # 扁率α
            self.alpha = (self.a-self.b)/self.a
            # 第一偏心率e1平方
            self.e1_2 = (self.a**2 - self.b**2)/self.a**2
            # 偏心率e2平方
            self.e2_2 = (self.a**2 - self.b**2)/self.b**2
         elif name == 'WGS84' or name == 'wgs84':
            pass
         else:
            # 长半轴a
            self.a = a
            # 短半轴b
            self.b = b
            # 极曲率半径c
            self.c = self.a**2/self.b
            # 扁率α
            self.alpha = (self.a-self.b)/self.a
            # 第一偏心率e1平方
            self.e1_2 = (self.a**2 - self.b**2)/self.a**2
            # 偏心率e2平方
            self.e2_2 = (self.a**2 - self.b**2)/self.b**2    

这里我们首先创建椭球类Ellipsoid()并初始化参数name椭球名称及可选参数椭球的长短半轴长a和b。

初始化函数可通过name判断椭球类型获取相应的椭球参数，这里pass部分省略了小秦定义的其他椭球的参数。

而最后的else后可通过更改长短半轴和椭球的名称定义自己想要的椭球参数。

Ellipsoid1 = Ellipsoid('CGCS2000')
# Ellipsoid1.a 可获取实例的属性

3. 椭球面上的常用坐标系及相互转换

为了表示椭球面上点的位置，必须要建立相应的坐标系。这里我们介绍几个常用的。

3.1 大地坐标系

哎呀，扫描的图上留下了小秦的笔记 

子午面与起始子午面的二面角L叫做大地经度，过P点椭球面的法线Pn与赤道面的夹角B叫做P点的大地经度，当然如果点不在椭球面上就需要再引入参数H——大地高。

3.2 大地空间直角坐标系

奇怪的笔记变多了 

以椭球中心O点为原点，起始子午面与赤道面交线为X轴赤道面上与X轴正交的方向为Y轴，旋转轴为Z轴。构成右手坐标系。

当然啦，常用的坐标系还有很多：子午面坐标系、地心纬度坐标系、大地极坐标系等等。

这里小秦主要介绍大地坐标系同空间直角坐标系坐标间的相互转换的实现方法(L,B,H)—(X,Y,Z)。

以下为代码实现:

1. 定义LBH，从而转换XYZ

def def_LBH(self,L,B,H=0):
        self.L = L
        self.B = B
        self.H = H
        self.L_h = d2h(L)
        self.B_h = d2h(B)
        self.LB = (L,B)
        self.LB_h = (d2h(L),d2h(B))
        # 辅助计算W
        self.W = math.sqrt(1 - self.e1_2*(math.sin(d2h(B)))**2)
        # 辅助计算V
        self.V = math.sqrt(1 + self.e2_2*(math.cos(d2h(B)))**2)
        # 法线长N，同时也是卯酉圈曲率半径N
        self.N = self.a/self.W
        # 空间直角坐标系三个参数
        self.X = (self.N+H)*math.cos(self.B_h)*math.cos(self.L_h)
        self.Y = (self.N+H)*math.cos(self.B_h)*math.sin(self.L_h)
        self.Z = (self.N*(1-self.e1_2)+H)*math.sin(self.B_h)
        self.XYZ = (self.X,self.Y,self.Z)

小秦这里在椭球类Ellipsoid()下定义了:

defdef_LBH(self,L,B,H=0)函数，H为可选参数获取椭球经纬度坐标

 

 1. 定义XYZ，从而转换LBH

 

def def_XYZ(self,X,Y,Z):
        self.X = X
        self.Y = Y
        self.Z = Z
        self.L_h = math.atan(Y/X)
        self.L = h2d(math.atan(Y/X))
        
        # 迭代求B
        B1 = Z/math.sqrt(X**2+Y**2)
        while True:
            # 辅助计算W
            self.W = math.sqrt(1 - self.e1_2*(math.sin(B1))**2)
            # 法线长/卯酉圈曲率半径N
            self.N = self.a/self.W
            B2 = math.atan((self.Z+self.N*self.e1_2*math.sin(B1)) / math.sqrt(X**2+Y**2))
            if abs(B2 - B1) <= 0.0000000001:
                break
            B1 = B2
        self.B_h = B1
        self.B = h2d(B1)
        self.LB = (self.L,self.B)

 

 发现了书上的小错误，继工程测量的第二次！

3.4  椭球面上的几种曲率半径

为了在椭球面上进行控制测量计算,就必须了解椭球面上有关曲线的性质。过椭球面上任意一点可作一条垂直于椭球面的法线,包含这条法线的平面叫做法截面,法截面同椭球面交线叫法截线(或法截弧)。可见,要研究椭球面上曲线的性质,就要研究法截线的性质,而法截线的曲率半径便是一个基本内容。

1.子午圈曲率半径

 

 2.卯酉圈曲率半径

 过椭球面上一点的法线,可作无限个法截面,其中一个与该点子午面相垂直的法面同椭球面相截形成的闭合的圈称为卯西圈。如图4-13中 PEE’即为过P点的卯西圈。卯西圈的曲率半径用N表示。

 3.平均曲率

所谓平均率半径R是指经过曲面任意一点所有可能方向上的法截线曲率半径R的算术平均值。

# 辅助计算W
        self.W = math.sqrt(1 - self.e1_2*(math.sin(d2h(B)))**2)
        # 辅助计算V
        self.V = math.sqrt(1 + self.e2_2*(math.cos(d2h(B)))**2)

 在计算前我们需要先定义辅助计算的参数方便计算

        # 椭球上的几种曲率半径属性
        # 子午圈曲率半径M
        self.M = self.N/self.V**2
        # 平均曲率半径R
        self.R = self.N/self.V
        
        # 子午线弧长arc_X：赤道开始到已知纬度点子午线上线段弧长
        self.arc_X = self.a_list[0]*self.B_h - math.sin(self.B_h)*math.cos(self.B_h)*(
            (self.a_list[1]-self.a_list[2]+self.a_list[3]
            )+(2*self.a_list[2]-16/3*self.a_list[3]
              )*(math.sin(self.B_h))**2 +16/3*self.a_list[3]*(math.sin(self.B_h))**4)
        # 平行圈半径r,同时也是子午面直角坐标x
        self.x = self.N*math.cos(self.B_h)
        # 子午面直角坐标y
        self.y = self.N*(1-self.e1_2)

 而后各曲率的计算都是十分简易的.

小秦这次的学习笔记就更新到这辣，后续还有弧长计算与大地线的知识哦。

当然摄影测量的笔记也不能少，谁让咱们摄影测量那么喜欢用矩阵计算呢，嘿嘿下次再见辣。

 

 可以在公众号阅读小秦的其他学习笔记和分享的学习资料哦。持续更新小秦的GIS学习之旅。