CN105043413A - 一种基于地球形体特征和大气参数的太阳位置的校正方法 - Google Patents

Abstract

本发明提供的一种基于地球形体特征和大气参数的太阳位置校正方法，通过太阳位置算法获取基准太阳位置，在基准太阳位置的基础上对地心视差和蒙气差进行逆向校正，使用校正后的数据对进行校正；地心视差逆向校正的方法为：当阳光照射到遮挡物上产生阴影时，计算遮挡物低端到地心的距离，根据这个距离、太阳高度角实现对地心视差校正，对蒙气差逆向校正时，对计算蒙气差时的每一层产生的光的折射过程进行再次进行细化，计算太阳光进入每一层的蒙气差，将此蒙气差带入到正向校正中得出的总的蒙气差的积分公式中，得出经过逆向校正的蒙气差公式。本发明基于地球形体及大气参数的校正能有效减小SPA原始计算误差，提高太阳位置的计算精度。

Description

一种基于地球形体特征和大气参数的太阳位置的校正方法

技术领域

本发明涉及太阳位置的计算，尤其涉及一种基于地球形体特征和大气参数的太阳位置校正方法。

背景技术

日照分析的基础算法是太阳位置算法(SolarPositionAlgorithm，简称SPA)，它主要根据分析时间和地理位置计算出某一时刻的太阳位置。

目前，在日照分析研究中主要存在以下几方面的问题：

(1)日照分析的研究对象方面

高低起伏的地形、具有一定高度和形体体征的植被、公用标志建筑、住宅周边的亭台楼阁，甚至被文化保护的古树、古塔等与居住建筑物同样对日照光线能产生遮挡效应。传统建筑日照分析研究中，仅以人居建筑为研究对象，忽略了真实条件下诸多复杂的地物要素，计算出的结果并不能逼真表达真实场景下的日照时数。

我国山地(包括丘陵)面积约占国土面积的三分之二(孙汉群,2005)，山地城镇约占全国城镇总数50％，山地居住的人口占全国人口50％。因此，人居环境周边的丘陵、山地、植被、标志建筑等地物要素对日照产生的遮挡效应不容忽视。

(2)日照分析所采用的太阳位置算法方面

在复杂地理场景中进行日照分析时，SPA原始计算误差将表现出更为显著的“倍级”放大效应。主要体现在：

①忽略了大地坐标与天文坐标二者的差异对太阳位置计算的误差影响；

②忽略地理坐标参照的地球形体特征对SPA计算精度影响。由坐标系参照的旋转椭球定义可知：地面点至地心的距离与大地坐标经纬值及参照的旋转椭球参数紧密相关。且由视差的产生机理可知：地面点至地心的距离是SPA视差校正的关键参数。因此，忽略SPA视差校正和坐标系参照的地球形体特征对SPA计算精度影响，势必导致日照分析计算结果的偏差增大；

③忽略了大气折射对SPA计算精度影响。由蒙气差产生机理可知，大气折射对日照光线方向上的影响效应与对应时刻下的太阳高度角紧密相关，且随高度角减小呈增大趋势。因此，忽略蒙气差校正将导致SPA计算误差进一步增加；

④针对多种SPA(包括TEP-SPA和NFP-SPA)对日照时数的计算精度影响性能，尚未开展过多研究，迄今还未发现精度适合于地学研究的SPA算法基准；

⑤视差和蒙气差的研究主要集中在天文和大气领域，且研究角度是从“观测位置”(即产生地物阴影的对应太阳位置，以下简称“有效太阳位置”)到“实际位置”(即SPA计算出的理想位置，以下简称“计算位置”)的正向校正，并未过多开展适用于日照分析的“逆向纠正”研究，即：“计算位置”到“有效位置”的校正研究。

因此，为尽可能减少复杂地理场景中日照分析计算误差，逼真地表达真实场景下的日照时数，在顾及地球表面复杂多样的地形要素之前提下，势必开展针对地理坐标***的最佳选择，以及坐标系参照椭球形体参数及大气折射对SPA的校正研究。

发明内容

本发明要解决的技术问题在于提供一种对地心视差和蒙气差的逆向校正的方法，以达到更为精确的计算太阳位置的目的。

为达到上述目的，本发明采用以下技术方案：

一种基于地球形体特征和大气参数的太阳位置校正方法，通过太阳位置算法获取基准太阳位置，在基准太阳位置的基础上对地心视差和蒙气差进行逆向校正，使用校正后的数据对太阳位置进行校正；

其中，

地心视差逆向校正的方法为：当阳光照射到高度为H的遮挡物上产生阴影时，从地心和遮挡物顶端两个地方观测到的太阳的高度角分别为h和h’，计算出h和h’之间的参数方程，再结合大地测量理论，计算遮挡物低端到地心的距离，将这个距离、现有太阳位置算法计算出的太阳高度角带入到上述的h和h’之间的参数方程中，实现对地心视差h’的校正，地心视差h’的校正公式为：

$h^{\′}=\arctan \left(\tan \left(\arcsin \left(\cos \mathrm{\ψ}\cos \mathrm{\δ}\cos t+\sin \mathrm{\ψ}\sin \mathrm{\δ}\right)\right)-\frac{OM}{OS\cos \left(\arcsin \left(\cos \mathrm{\ψ}\cos \mathrm{\δ}\cos t+\sin \mathrm{\ψ}\sin \mathrm{\δ}\right)\right)}\right)$ ；

OM为遮挡物MM′上高度为H的顶点M与地心之距离，OS为某一时刻t日地距离，OS为某一时刻太阳质心至地球质心距离，Ψ为当地经度(度)，δ为太阳 光线垂直照射的地点与地球赤道所夹的圆心角；

蒙气差逆向校正的方法为：使用太阳实际位置对应的高度角计算蒙气差；具体来说，对蒙气差进行正向校正时，将蒙气差正向校正中对大气层分割成密度均匀的许多层后，对每一层中的光线线路进行分段线性处理，得到每一层的蒙气差，进而得到总的蒙气差的积分公式，而对蒙气差逆向校正时，则对每一层产生的光的折射过程进行再次进行细化，即根据折射原理，计算太阳光进入每一层的蒙气差，将此蒙气差带入到正向校正中得出的总的蒙气差的积分公式中，得出经过逆向校正的蒙气差公式：

$\mathrm{\ζ}=\sin V{\∫}_1^L\left(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V}}+\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{\left({\mathrm{\μ}}_i^2-{\sin }^{2}V\right)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^2}{8\sqrt{{\left({\mathrm{\μ}}_i^2-{\sin }^{2}V\right)}^5}}+...\right)\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i},$ V为计算出的太阳天顶角，即V＝90°-h，L为地面折射率，S为观测到的太阳在天空中的位置，μ_i为第i大气层的折射率。

所述获取基准太阳位置的太阳位置算法为基于天文算法的SPA算法，所述的太阳位置指的是方位角和高度角。

所述地心视差逆向校正的公式中遮挡物顶端到地心的距离OM为遮挡物的高度H加上遮挡物底端到地心的距离OM′，其中遮挡物底端到地心的距离OM′满足：

${\mathrm{OM}}^{\′}=R_{OP}=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}$

遮挡物顶端到地心的距离OM满足：

$OM=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}+H,$ H为遮挡物MM′高度，e为椭球曲

率，N为P点卯酋圈曲率半径，B为大地纬度。

所述某一时刻太阳质心到地球质心的距离OS为：

OS≈15210万千米(远日点，4月16日-10月15日)，

OS≈14710万千米(近日点，前一年10月16日至4月15日)。

蒙气差正向校正时的蒙气差微分公式为:

$\mathrm{\ζ}={\mathrm{\μ}}_1{\mathrm{sinv}}_1{\∫}_1^{{\mathrm{\μ}}_1}\left(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1}}-\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^4}{8\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^5}}+...\right)\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i},$

其中为观测太阳在天空中的位置，r_i为第i层和第i-1层的球面边界的半径；r_i-1为第i-1和i-2层的球面边界半径，v₁为地面观测太阳天顶角，μ_i为积分项对应的第i层大气折射率，μ₁为地面折射率。

对蒙气差的逆向校正公式的简化公式为：

其中，μ₁为地面折射率，h为地心视差。

在任何温度、气压条件下，地面太阳可见光的大气折射率为：

其中，T为开氏温度，t为摄氏温度，T＝t+273.15，P为大气压强。

本发明的有益效果：本发明基于地球形体及大气参数的SPA校正能有效减小SPA原始计算误差，提高太阳位置的计算精度。

附图说明

图1为天球模型示意图。

图2为地心视差原理示意图。

图3为有效太阳高度角与计算太阳高度角之间的地心视差示意图。

图4为空间直角坐标系下地面P点坐标。

图5为大气层分层折射表达。

图6为有效太阳高度角与计算太阳高度角之间的蒙气差。

图7为球面大气层入射角与折射角的关系。

图8为平行平面大气层蒙气差。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明。

本发明使用天文算法即《AstronomicalAlgorithms》(Meeus,1998)公布的太阳位置计算方法(即TEP-SPA的代表算法)作为基准SPA算法，通过天文算法计算太阳方位角A、高度角h。

概念解释：

“观测位置”或者“有效太阳位置”：即产生地物阴影的对应太阳位置，本 发明中是指针对日照分析而言，包含视差和蒙气差影响的太阳有效位置；

“实际位置”或“计算位置：”即SPA计算出的理想位置；

“逆向纠正”即“计算位置”到“有效位置”的校正。

符号含义如表1-1所示：

表1-1本发明中使用的符号的参数含义。

如图1所示为天球模型示意图，根据天球坐标系之间的映射关系，可以构建出地平与时角坐标关系式：

$\left\{\begin{array}{c}\cos A=\frac{\sin \mathrm{\ψ}\cos \mathrm{\δ}\cos t-\cos \mathrm{\ψ}\sin \mathrm{\δ}}{\cosh }\\ \sin A=\frac{\cos \mathrm{\δ}\sin t}{\cosh }\\ \sinh =\cos \mathrm{\ψ}\cos \mathrm{\δ}\cos t+\sin \mathrm{\ψ}\sin \mathrm{\δ}\end{array}\right.---\left(1\right)$

其中t＝15°×(Tr-12)，Ψ为当地理经度(东经取“+”，西经取“-”)。

在地心视差的校正中，如图2所示，当观测者在地球表面之上观测太阳高度角时，其观测值是以“在地面上”而言的，而天体坐标系常采用地心作为坐标系原点，因此，欲建立二者的对应，必须把地面上观测值归化至地球中心之相当数 值，由此而引起的高度角差值即为地心视差(简称“视差”)，此校正过程即为视差校正。

S为观测太阳在天球中的位置，Z′为在地球面上一点M观测太阳的天顶角(距)，Z为在地球中心所见S的天顶角(距)，因此，视差P即为：P＝Z′-Z，根据三角形ΔMOS的正弦定理可得：

$\frac{\sin P}{OM}=\frac{{\mathrm{sinZ}}^{\′}}{OS}$

其中，OM为观测点至地心的距离，OS为某一时刻日地距离。如果太阳在观察点地平上，其视差为P₀，则有：

${\mathrm{sinP}}_0=\frac{OM}{{\mathrm{OS}}_0}$

由于天球为正球体，因此则有：OS₀＝OS

因此，

sinP＝sinP₀sinZ′

又由于P及P₀均为极小角，sinP≈P,sinP₀≈P₀，所以：

P＝P₀sinZ′＝P₀cosh′

其中h′为地面M点观测的太阳高度角，与天顶距Z'互余，即：Z'＝90°-h'。

如图2所示，在复杂地理环境中，产生日照阴影的太阳高度角h′(天顶距Z′的余角)是基于地面而言，根据SPA计算出的太阳高度角h(天顶距Z的余角)是基于地心而言，需要对“计算位置”下的太阳高度角进行视差反向校正，即本发明中对地心视差的逆向校正过程。

如图3所示，S为太阳在天球中的实际位置，h为SPA计算出的太阳高度角，是基于地心O而言；在地面上，当阳光实际照射到遮挡物MM′产生遮挡阴影时，太阳高度角即为h′。假设遮挡物MM′高度为H，Z′和Z分别是高度角h′和h对应的太阳天顶距。

由三角函数关系可知： $Z^{\′}=Z+P=\frac{\mathrm{\π}}{2}-h^{\′};Z=\frac{\mathrm{\π}}{2}-h,$ 即：

P＝h-h′

根据三角形ΔOMS正弦定理可知：

$\frac{\sin P}{OM}=\frac{{\mathrm{sinZ}}^{\′}}{OS}=\frac{{\cosh }^{\′}}{OS}$

即：

$\sin P=\frac{OM}{OS}{\cosh }^{\′}$

因此有：

$sin(h-h^{\′})=\frac{OM}{OS}{\cosh }^{\′}$

即：

$h^{\′}=arctan(\tanh -\frac{OM}{OS\cosh })---\left(2\right)$

其中，OM为遮挡物MM′上高度为H的顶点M与地心之距离，OS为某一时刻日地距离。因此，如果求出遮挡物MM′地面点M′距地心距离，即可求出OM值。

如图4所示，在地心空间直角坐标系下，如果地球表面上某点P对应上述地面点M′，根据大地测量理论，地面点P坐标(X，Y，Z)与P(B，L)关系存在以下关系：

$\left\{\begin{array}{c}X=N\cos B\cos L\\ Y=N\cos B\sin L\\ Z=N(1-e^2)\sin B\end{array}\right.$

其中N为P点卯酋圈曲率半径，(a为椭球长半轴)，e为椭球曲率。如果把P点的大地纬度B转换为球心纬度为φ，则有：

tanφ＝(1-e²)tanB

此时P点至椭球中心距离R_OP为：

即：

$R_{OP}=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}$

正如前述，根据前述假设三，在WGS84坐标系下，椭球的质心即为椭球几何中心，所以：

${\mathrm{OM}}^{\′}=R_{OP}=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}$

由OM＝OM′+MM′，因此：

$OM=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}+H---\left(3\right)$

至于某一时刻太阳质心至地球质心距离OS的计算相对比较复杂，目前比较精确的数据主要来自天文观测数据，也可采用开普勒第一、第二定律计算。而针对于日照分析计算模型而言，在进行视差校正时，地球半径远小于日地距离，因此，可以采用日地平均距离代替某一时刻日地距离。日地距离最大值为15210万千米(北半球而言，地球处于远日点，每年7月4日前后)；最小值为14710万千米(北半球而言，地球处于近日点，每年1月3日)；平均值为14960万千米，也称为一个天文单位(AU)。1976年国际天文学联合会IAU确定为149597870千米。因此，当计算日期处于某一范围时，可采用对应的日地平均距离代替OS，即：

OS≈15210万千米(远日点，4月16日-10月15日)

OS≈14710万千米(近日点，前一年10月16日至4月15日)(4)

由SPA的基本公式(1)可知：

sinh＝cosψcosδcost+sinψsinδ

也就是：

h＝arcsin(cosψcosδcost+sinψsinδ)

进行视差逆向校正后的地心视差h′即为：

$h^{\′}=arctan\left(tan\left(arcsin\left(cos\mathrm{\ψ}cos\mathrm{\δ}\cos t+sin\mathrm{\ψ}sin\mathrm{\δ}\right)\right)-\frac{OM}{OScos\left(arcsin\left(cos\mathrm{\ψ}\cos t+sin\mathrm{\ψ}sin\mathrm{\δ}\right)\right)}\right)$ ，

其中，OM和OS可分别由公式(3)和公式(4)计算得出。

接下来是是对蒙气差的逆向校正过程。

在真空中或在均匀的透明介质中光遵循直线传播。但是实际上，由于地球的大气是不均匀的，太阳光通过大气层时便受到曲折，因而我们所观测到的太阳方向与其视方向之间会形成一个小的夹角，这种现象称为大气折射，而由大气折射形成的夹角称为蒙气差。

在蒙气差的正向校正中，为方便，如图5所示，把大气层分割为密度均匀的N+1层，其中第N+1层为外太空，其折射率μ_N+1＝1；接近地平的大气层为第1 层，其折射率为μ₁。

如图5所示，O为地心，M为地球上观测点，Z为天顶，K为太阳光线入射大气层的入射点，P′为太阳视位置，即太阳在地球上看到的位置；P为太阳实际角度，即太阳在天球上的实际位置；μ为大气层每一层对应的折射率；Δξ为由于大气层对阳光折射所导致的每层角度损耗，即为每一大气层对应蒙气差。

由于在对大气层分层时假设每一层的大气密度是均匀的，因此太阳光线在每一层遵循直线原则。自M点做折线KM的各小段的平行线，这样就把总蒙差角ξ分为许多小角Δξ₁，Δξ₂，Δξ₃…Δξ_N，其中每一个Δξ就是光线在对应大气层边界上所发生的曲折角度，因此：

ξ＝Δξ₁+Δξ₂+Δξ₃+…+Δξ_N

即：当n→∞时，则可用定积分表示蒙气差ξ为：

$\mathrm{\ξ}={\∫}_1^{N}d\mathrm{\ξ}---\left(5\right)$

根据大气折射入射角和折射角关系以及Snell定律，对公式(5)展开可以得到：

$\mathrm{\ζ}={\mathrm{\μ}}_1{\mathrm{sinv}}_1{\∫}_1^{{\mathrm{\μ}}_1}(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1}}-\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^4}{8\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^5}}+...)\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}---\left(6\right)$

其中，r_i为第i层和第i-1层的球面边界的半径；r_i-1为第i-1和i-2层的球面边界半径，v₁为地面观测太阳天顶角，μ_i为积分项对应的第i层大气折射率，μ₁为地面折射率。公式(6)即为蒙气差基本微分公式。

但是，SPA计算出的太阳高度角为太阳在天球中的“实际位置”(即文中“计算位置”)，而实际产生阴影的太阳高度角(即文中“有效位置”)将受到蒙气差(ξ)的影响。然而，在天文学研究中，主要基于星体的“有效位置”下高度角(h′)计算蒙气差(ξ)(图6所示)。而在日照分析中，则需要由“计算位置”对应的高度角(h)计算蒙气差，因此，此校正过程称为蒙气差的逆向校正。

在蒙气差的正向校正过程中，已经得出了蒙气差的基本微分公式(公式6)，在你逆向正中，为了方便，对某一层大气产生的折射过程再次进行细化(如图7 所示)。

如图7中所示，圆弧表示第i层大气层的球面边界，太阳光线在A点产生折射后沿着AC方向前进，AZ是折射面的法线，v，f分别是入射角和折射角，dζ为此大气层蒙气差。

可得

v＝f+dζ

设μ_i，μ_i-1分别为入射空气层和折射空气层的折射率，即：

μ_i-1＝μ_i+dμ

根据笛卡尔定律(也称为Snell定律)可以推导出：

$\frac{\sin v}{\sin f}=\frac{{\mathrm{\μ}}_{i-1}}{{\mathrm{\μ}}_i}=\frac{\sin v}{sin(v-d\mathrm{\ζ})}=\frac{{\mathrm{\μ}}_i+d\mathrm{\μ}}{{\mathrm{\μ}}_i},$ 即： $\frac{sin\mathrm{\ν}}{\sin v\cos d\mathrm{\ζ}-\cos v\sin d\mathrm{\ζ})}=1+\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}$

根据微分公式可知，当dζ为小量时，sindζ＝dζ；cosdζ＝1，上式即变化为：

$\frac{1}{1-d\mathrm{\ζ}\cot v}=1+\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i},$ 即：μ_idζ+dμdζ＝dμtanv

如果略去二阶小量dμdζ，则：μ_idζ＝dμtanv，即：由把公式(7)代入公式(5)可得到：

在公式(8)中，积分下限是大气层外(真空)的折射率μ_N为1,上限是地面空气层的折射率μ₁，用L代替，v为第i层大气层的阳光入射角，μ_i为第i层大气层的折射率。不妨假设用v_i表示v，用f_i表示f。于是蒙气差ξ的求取就变为求取v_i和μ_i的关系，因此公式(8)可描述为：

$\mathrm{\ζ}={\∫}_1^L\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}{\mathrm{tanv}}_i---\left(9\right)$

如图7所示，r_N为第N层和第N-1层的球面边界的半径；r_N-1为第N-1和N-2层的球面边界半径；v_i为第i层入射角；L代表为地面层大气折射率μ₁；f_i由第i层大气折射产生的折射角，根据笛卡尔定律(Snell折射定理)可知：

$\frac{{\mathrm{sinv}}_N}{{\mathrm{sinf}}_N}=\frac{{\mathrm{\μ}}_{N-1}}{{\mathrm{\μ}}_N}---\left(10\right)$

根据三角形ΔKOB的边角关系可知：

$\frac{{\mathrm{sinf}}_N}{r_{N-1}}=\frac{{\mathrm{sinv}}_{N-1}}{r_N}---\left(11\right)$

由公式(10)和公式(11)相乘得出：

μ_Nr_Nsinv_N＝μ_N-1r_N-1sinv_N-1(12)

公式(12)所表达的含义为：任意大气层折射率μ_i，球面边界半径r_i与入射角v_i的正弦的乘积为常数，这既是球面斯涅耳定律。因此，任意层的折射率μ_i，球面边界半径r_i与入射角v_i和外大气层的折射率μ_N、球心半径r_N和入射角v_N存在以下关系：

μ_ir_isinv_i＝μ_Nr_Nsinv_N＝r_NsinV

其中V为计算出的理论太阳天顶角，即V＝90°-h，r_N也可理解为外大气层至地心半径。

因此： ${\mathrm{tanv}}_i=\frac{\frac{r_N}{{\mathrm{\μ}}_i{r}_i}\sin V}{\sqrt{1-{\left(\frac{r_N}{{\mathrm{\μ}}_i{r}_i}\sin V\right)}^2}}---\left(13\right)$

由公式(9)和公式(13)可得：

$\mathrm{\ζ}={\∫}_1^L\frac{\frac{r_N}{{\mathrm{\μ}}_i{r}_i}\sin V}{\sqrt{1-{\left(\frac{r_N}{{\mathrm{\μ}}_i{r}_i}\sin V\right)}^2}}\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}$

即： $\mathrm{\ζ}={\∫}_1^L\frac{r_N\sin V}{\sqrt{{\mathrm{\μ}}_i^2{r}_i^2-r_N^2{\sin }^{2}V}}\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}---\left(14\right)$

如果令由于S是小量，有理由认为S²＝0，同时,μ_N＝1,即：

$\mathrm{\ζ}={\∫}_1^L\frac{\sin V}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V-2{\mathrm{\μ}}_i^{2}s}}\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}$

即： $\mathrm{\ζ}=\sin V{\∫}_1^L\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V-2{\mathrm{\μ}}_i^{2}s}}\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}$

令 $a={\mathrm{\μ}}_i^2-{\sin }^{2}V;b=2{\mathrm{\μ}}_i^{2}s,$ 则由二项式展开公式：

${(a+b)}^n=a^n+{\mathrm{na}}^{n-1}b+\frac{n(n-1)}{2!}{a}^{n-2}{b}^2+\frac{n(n-1)(n-2)}{3!}{a}^{n-3}{b}^3+...+b^n$

可得： $\mathrm{\ζ}=\sin V{\∫}_1^L\left(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V}}+\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{\left({\mathrm{\μ}}_i^2-{\sin }^{2}V\right)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^4}{8{\sqrt{\left({\mathrm{\μ}}_i^2-{\sin }^{2}V\right)}}^5}+...\right)\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}---\left(15\right)$

公式(14)即为蒙气差逆向校正基本微分公式,其中，V为计算出的太阳天顶角，即V＝90°-h，L为地面折射率，S与第i大气层边界至地心距离r_i与大气层外边界至地心距离r_N有关，μ_i为第i大气层的折射率。

在公式(15)中，如果仅保留第一项，则蒙气差ξ可表达为：

$\mathrm{\ζ}=\sin V{\∫}_1^L\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V}}\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i}---\left(16\right)$

令a＝sinV，则公式(16)可化简为：

又根据积分公式：

展开可以得到： $\mathrm{\ζ}={\[arc\sec \frac{\mathrm{\μ}}{\sin V}+C\]}_1^L=V-arcsin\left(\frac{\sin V}{L}\right)$

因为V＝90°-h，L为地面层大气折射率μ₁，因此：

$\mathrm{\ζ}=arccos\left(\frac{\cosh }{{\mathrm{\μ}}_1}\right)-h---\left(17\right)$

公式(17)可以认为是蒙气差逆向校正的微分简化公式，从公式推导过程可知，公式(17)简化了大气曲率的变换对大气蒙气差的影响，把地球大气层简化为平行平面大气层。

因此，如果直接把大气层简化为平行平面大气层很容易推导出公式(17)。

如图8所示，由Snell定理可知：μ_nsinZ_n＝...＝μ₀sinZ₀

即：sinz＝μ₀sinz'。把天顶距转换为高度角，则有：cosh＝μ₀cosh′

即太阳光线折射以后的高度角为：

则蒙气差ξ为：与上述公式(17)相同。

蒙气差逆向校正公式中具有参数地面大气折射率μ0，因此，需要知道此参数的大小，才能对蒙气差进行逆向校正。

根据现有的折射率公式，在标准大气条件下(760毫米汞柱，0℃)，折射率μ₀[标]为：

上述公式略去了水汽气压对折射率的影响，其中λ光波波长，太阳可见光波长范围为0.39-0.77μm，可取中间值0.58μm，因此，标准大气条件下，地面大气折射率在标准大气条件下μ₀[标]可近似取值为1.000292566。

根据理想气体状态方程，空气密度δ是压强P，温度T的函数，即：

(其中m为大气平均分子质量，R为大气普适常数)

而根据道尔-格拉斯通(Dale-Gladstone)定理，大气折射率μ与大气密度δ成线性关系，即：

μ-1＝Cδ(其中C为常数)

令：Δ＝(μ-1)，则有：

因此，在任何温度、气压条件下，地面太阳可见光大气折射率为：

其中，T为开氏温度，t为摄氏温度，T＝t+273.15。

当然，上述蒙气差的计算过程是基于理想的大气模型进行的，但是实际的大气情况是非常复杂的，大气的折射不仅与温度、气压、水汽等因素密切相关，而且大气的等密度层并不是同心球层或者水平平面，一定会存在大气层曲率变化，因此，在计算时，蒙气差的计算只是逼近近似。因此，也可以对蒙气差进行拟合 计算，对蒙气差进行拟合计算也会有较高的精度。

本发明中基于有理函数逼近思想，采用公式(18)函数形式，借助Mathematics工具软件进行非线性回归，最终得出p₁...p₇系数。其中，x表示“计算位置”下高度角，单位为度，p₁...p₇为待定系数，f(x)为蒙气差，单位为秒。

$f\left(x\right)=\frac{(p_1+p_3*x+p_5*x^2+p_7*x^3)}{(1+p_2*x+p_4*x^2+p_6*x^3)}---\left(18\right)$

拟合的数据选择我国《2009年天文年历》的观测数据，得出最终的蒙气差校正公式：

$\mathrm{\ξ}=\frac{\left(1819.08371242143+194.887513592849*h+1.46555397475109*h^2-0.0419553783815395*h^3\right)}{\left(1+0.409283439734292*h+0.0667313795916436*h^2+0.00008468597*h^3\right)}$

其中，h为计算太阳高度角(°)，ξ为蒙气差(″)。

Claims (7)

1.一种基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：通过太阳位置算法获取基准太阳位置，在基准太阳位置的基础上对地心视差和蒙气差进行逆向校正，使用校正后的数据对太阳位置进行校正；

其中，

地心视差逆向校正的方法为：当阳光照射到高度为H的遮挡物上产生阴影时，从地心和遮挡物顶端两个地方观测到的太阳的高度角分别为h和h’，计算出h和h’之间的参数方程，再结合大地测量理论，计算遮挡物低端到地心的距离，将这个距离、现有太阳位置算法计算出的太阳高度角带入到上述的h和h’之间的参数方程中，实现对地心视差h’的校正，地心视差h’的校正公式为：

$h^{\′}=\arctan (\tan \left(\arcsin \right(\cos \mathrm{\ψ}\cos \mathrm{\δ}\cos t+\sin \mathrm{\ψ}\sin \mathrm{\δ}\left)\right)-\frac{OM}{OS\cos \left(\arcsin \right(\cos \mathrm{\ψ}\cos \mathrm{\δ}\cos t+\sin \mathrm{\ψ}\sin \mathrm{\δ}))});$

OM为遮挡物MM′上高度为H的顶点M与地心之距离，OS为某一时刻t日地距离，OS为某一时刻太阳质心至地球质心距离，Ψ为当地经度(度)，δ为太阳光线垂直照射的地点与地球赤道所夹的圆心角；

蒙气差逆向校正的方法为：使用太阳实际位置对应的高度角计算蒙气差；具体来说，对蒙气差进行正向校正时，将蒙气差正向校正中对大气层分割成密度均匀的许多层后，对每一层中的光线线路进行分段线性处理，得到每一层的蒙气差，进而得到总的蒙气差的积分公式，而对蒙气差逆向校正时，则对每一层产生的光的折射过程进行再次进行细化，即根据折射原理，计算太阳光进入每一层的蒙气差，将此蒙气差带入到正向校正中得出的总的蒙气差的积分公式中，得出经过逆向校正的蒙气差公式：

$\mathrm{\ζ}=\sin V{\∫}_1^L(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\sin }^{2}V}}+\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{({\mathrm{\μ}}_i^2-{\sin }^{2}V)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^4}{8\sqrt{{({\mathrm{\μ}}_i^2-{\sin }^{2}V)}^5}}+\mathrm{...})\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i},$ V为计算出的太阳天顶角，即V＝90°-h，L为地面折射率，S为观测到的太阳在天空中的位置，μ_i为第i大气层的折射率。

2.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：所述获取基准太阳位置的太阳位置算法为基于天文算法的SPA算法，所述的太阳位置指的是方位角和高度角。

3.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置算法的校正方法，其特征在于：所述地心视差逆向校正的公式中遮挡物顶端到地心的距离OM为遮挡物的高度H加上遮挡物底端到地心的距离OM′，其中遮挡物底端到地心的距离OM′满足：

${\mathrm{OM}}^{\′}=R_{OP}=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}$

遮挡物顶端到地心的距离OM满足：

$OM=N\sqrt{{\cos }^{2}B+{(1-e^2)}^2{\sin }^{2}B}+H,$ H为遮挡物MM′高度，e为椭球曲

率，N为P点卯酋圈曲率半径，B为大地纬度。

4.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：所述某一时刻太阳质心到地球质心的距离OS为：

OS≈15210万千米(远日点，4月16日-10月15日)，

OS≈14710万千米(近日点，前一年10月16日至4月15日)。

5.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：蒙气差正向校正时的蒙气差微分公式为:

$\mathrm{\ζ}={\mathrm{\μ}}_1{\mathrm{sinv}}_1{\∫}_1^{{\mathrm{\μ}}_1}(\frac{1}{\sqrt{{\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1}}-\frac{{\mathrm{s\μ}}_i^2}{\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^3}}+\frac{3s^2{\mathrm{\μ}}_i^4}{8\sqrt{{({\mathrm{\μ}}_i^2-{\mathrm{\μ}}_1^2{\sin }^2{v}_1)}^5}}+\mathrm{...})\frac{d\mathrm{\μ}}{{\mathrm{\μ}}_i},$

其中为观测太阳在天空中的位置，r_i为第i层和第i-1层的球面边界的半径；r_i-1为第i-1和i-2层的球面边界半径，v₁为地面观测太阳天顶角，μ_i为积分项对应的第i层大气折射率，μ₁为地面折射率。

6.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：对蒙气差的逆向校正公式的简化公式为：

其中，μ₁为地面折射率，h为地心视差。

7.根据权利要求1所述的基于地球形体特征和大气参数的太阳位置校正方法，其特征在于：在任何温度、气压条件下，地面太阳可见光的大气折射率为：

其中，T为开氏温度，t为摄氏温度，T＝t+273.15，P为大气压强。