CN104715143A - 基于卫星图像数据的临边像素点的地球切线高度估计方法 - Google Patents

Abstract

本发明公开了基于卫星图像数据的临边像素点的地球切线高度估计方法。该方法包括：获取卫星位置信息；从卫星遥感数据中获取临边像素点探测器的指向信息；根据卫星位置信息和临边像素点探测器的指向信息，计算临边像素点的位置；计算临边像素点的切线高度。本发明所提供的基于卫星图像数据获得临边像素点的地球切线高度的方法，不受临边星观测的限制，并且考虑了地球的椭球形状，因此计算的结果更为准确。

Description

基于卫星图像数据的临边像素点的地球切线高度估计方法

技术领域

本发明涉及大气红外辐射技术，特别涉及一种基于卫星图像数据估计临边像素点的地球切线高度的方法。

背景技术

大气层既是地球气候与环境的主要载体和活动舞台,又是空间天气与环境的重要组成部分。大气层下连地球海洋及生物圈,上接太阳活动重大影响区域磁层和电离层,上下发生着重要的相互作用。对大气层进行遥感探测从而发现和理解其中的整体行为和相互作用,一直是地球和空间科学界共同的追求目标。

根据星载探测器、探测目标(地球表面或大气)和信号源(太阳辐射、地球/大气辐射等)位置的几何关系来进行划分，卫星探测一般有三种主要探测方式：下视对地观测、掩星(掩日/月)观测和临边观测。其中，临边观测方式是一种新的观测方式，它既具有天底观测方式(即下视对地观测)的高空间覆盖范围，又同时具有掩星观测方式的高垂直分辨率，并且能够有效避开复杂地表和对流层云和雨等强信号的干扰。因此，开展临边大气红外临边辐射特性研究，了解其光谱与垂直分布规律,对于大气成分的监测等具有重要意义。

传统的测量临边高度的方法一般是利用星光折射角获得。当星光穿过地球附近大气层时，星光会被折射，利用星敏感器的观测结果和星表可以计算得到星光折射角；再通过大气层对称假设和光的传输路径大气的透过率模型以及星与观测点之间几何关系，可以计算得到切线高度。但是，上述的计算过程中受到可观测星的位置限制和大气透过率模型的限制。

发明内容

有鉴于此，本发明提供一种基于卫星图像数据获得临边像素点的地球切线高度的方法，从而可以有效地提高计算遥感数据中每个地球临边像素所处的地球切线高度的准确度。

本发明的技术方案具体是这样实现的：

一种基于卫星图像数据获得临边像素点的地球切线高度的方法，该方法包括：

获取卫星位置信息；

从卫星遥感数据中获取临边像素点探测器的指向信息；

根据卫星位置信息和临边像素点探测器的指向信息，计算临边像素点的位置；

计算临边像素点的切线高度。

较佳的，所述获取卫星位置信息包括：

从预报星历或精密星历中读取所需的卫星位置信息。

较佳的，所述获取卫星位置信息包括：

根据卫星轨道根数计算探测器成像时刻对应的平近点角；

根据卫星轨道根数和探测器成像时刻对应的平近点角，计算探测器成像时刻对应的偏近点角；

根据卫星轨道根数和探测器成像时刻对应的偏近点角，计算卫星在地球中心惯性系中的位置矢量。

较佳的，所述卫星轨道根数中包括：

历元时刻t₀、轨道半长轴a₀、历元时刻t₀对应的平近点角m₀、轨道偏心率e₀、轨道倾角i₀、升交点赤经Ω₀和近地点幅角ω₀。

较佳的，根据如下的公式计算探测器成像时刻对应的平近点角m_j：

m_j＝m₀+(GM/a₀ ³)^1/2(t_j-t₀)      (1)

其中，GM为地球引力常数，t_j为探测器成像时刻。

较佳的，所述计算探测器成像时刻对应的偏近点角包括：

令E_j ⁽⁰⁾＝m_j；；

根据公式：E_j ⁽ⁿ⁺¹⁾＝E_j ⁽ⁿ⁾–(E_j ⁽ⁿ⁾-e₀*sin E_j ⁽ⁿ⁾-m_j)/(1-e₀*cos E_j ⁽ⁿ⁾)进行反复迭代，直到满足条件：|E_j ⁽ⁿ⁺¹⁾-E_j ⁽ⁿ⁾|≤σ(σ＝1e-6)；

若n≥20仍不满足上述条件时，则取：E_j＝m_j+e₀*sin(m_j)+e₀ ²sin(2m_j)；

其中，e₀为轨道偏心率；n为迭代次数；σ为迭代阈值。

较佳的，根据如下的公式计算卫星在地球中心惯性系中的位置矢量R：

R＝(x_p，y_p，z_p)＝a₀(cosE_j-e₀)P+a₀(1-e₀ ²)^1/2sinE_jQ；

其中，P＝(cosω₀cosΩ₀-sinω₀sinΩ₀cosi₀，cosω₀csinΩ₀+sinω₀sinΩ₀cosi₀，sinω₀sini₀)，

Q＝(-sinω₀cosΩ₀-cosω₀sinΩ₀cosi₀，-sinω₀csinΩ₀+cosω₀cosΩ₀cosi₀，cosω₀sini₀)；

其中，i₀为轨道倾角，Ω₀为升交点赤经，ω₀为近地点幅角。

较佳的，所述临边像素点探测器的指向信息为临边像素点探测器的方向矢量L：

L＝(cosαcosβ,cosαsinβ,sinα)＝(m，n，p)。

其中，α为方位角，β为俯仰角。

较佳的，根据如下的公式计算临边像素点N的位置：

(x_n，y_n，z_n)＝(x_p+mt，y_p+nt，z_p+pt)；

其中，(x_n，y_n，z_n)为临边像素点N的位置，t＝-(mx_p+ny_p+pz_p)。

如上可见，在本发明中的基于卫星图像数据获得临边像素点的地球切线高度的方法中，可以利用卫星遥感数据中卫星的位置和指向等信息，在地球坐标系中计算遥感数据中每个地球临边像素所处的切线高度。与现有技术中的方法相比，本发明中所提出的方法不受临边星观测的限制，并且考虑了地球的椭球形状，因此计算的结果更为准确。

附图说明

图1为本发明实施例中的基于卫星图像数据获得临边像素点的地球切线高度的方法的流程示意图。

图2为本发明实施例中的计算临边像素点的位置的示意图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举实施例，对本发明进一步详细说明。

本实施例提供了一种基于卫星图像数据获得临边像素点的地球切线高度的方法。

图1为本发明实施例中的基于卫星图像数据获得临边像素点的地球切线高度的方法的流程示意图。如图1所示，本发明实施例中的基于卫星图像数据获得临边像素点的地球切线高度的方法主要包括如下所述的步骤：

步骤101，获取卫星位置信息。

在本发明的技术方案中，首先可以从卫星遥感数据中获取卫星位置信息。

较佳的，在本发明的具体实施例中，所述卫星位置信息可以是卫星在地球中心惯性系中的坐标(x_p,y_p,z_p)。

另外，在本发明的技术方案中，可以通过多种方式来实现上述的步骤101。以下将以其中的两种具体实现方式为例，对本发明的技术方案进行介绍。

第一种实现方式：从预报星历或精密星历中读取所需的卫星位置信息。

在预报星历或精密星历中，一般都有卫星位置信息，因此，可以直接从预报星历或精密星历中读取所需的卫星位置信息。例如，直接从预报星历或精密星历中读取卫星在地球中心惯性系中的坐标(x_p,y_p,z_p)。

第二种实现方式：根据卫星遥感数据获取卫星位置信息。

在卫星遥感数据(例如，卫星图像数据)中，一般都会提供一些辅助信息，该辅助信息中包含有卫星瞬时轨道根数，因此，可以根据卫星瞬时轨道根数计算得到卫星位置信息。例如，较佳的，在本发明的具体实施例中，可以通过如下所述的步骤来计算卫星位置信息：

步骤11，根据卫星轨道根数计算探测器成像时刻对应的平近点角。

较佳的，在本发明的具体实施例中，所述卫星轨道根数中可以包括：历元时刻t₀、轨道半长轴a₀、历元时刻t₀对应的平近点角m₀、轨道偏心率e₀、轨道倾角i₀、升交点赤经Ω₀和近地点幅角ω₀。因此，可以根据上述卫星轨道根数计算探测器成像时刻对应的平近点角。

较佳的，在本发明的具体实施例中，可以根据如下所述的公式计算探测器成像时刻对应的平近点角m_j：

m_j＝m₀+(GM/a₀ ³)^1/2(t_j-t₀)     (1)

其中，上述公式(1)中的m₀为历元时刻t₀对应的平近点角，GM为地球引力常数(一般取值为3.986006e14)，a₀为轨道半长轴，t_j为探测器成像时刻，t₀为历元时刻(包括年月日时分秒的信息)。

步骤12，根据卫星轨道根数和探测器成像时刻对应的平近点角，计算探测器成像时刻对应的偏近点角。

较佳的，在本发明的具体实施例中，可以根据如下所述的步骤计算探测器成像时刻对应的偏近点角E_j：

步骤121，令E_j ⁽⁰⁾＝m_j；

步骤121，根据公式：E_j ⁽ⁿ⁺¹⁾＝E_j ⁽ⁿ⁾–(E_j ⁽ⁿ⁾-e₀*sin E_j ⁽ⁿ⁾-m_j)/(1-e₀*cos E_j ⁽ⁿ⁾)进行反复迭代，直到满足条件：|E_j ⁽ⁿ⁺¹⁾-E_j ⁽ⁿ⁾|≤σ(σ＝1e-6)；

若n≥20仍不满足上述条件时，则取：E_j＝m_j+e₀*sin(m_j)+e₀ ²sin(2m_j)。

其中，e₀为轨道偏心率；n为迭代次数；σ为迭代阈值，两次迭代的差值小于σ时迭代结束。

步骤13，根据卫星轨道根数和探测器成像时刻对应的偏近点角，计算卫星在地球中心惯性系中的位置矢量。

较佳的，在本发明的具体实施例中，可以根据如下所述的公式计算卫星在地球中心惯性系中的位置矢量R：

R＝(x_p，y_p，z_p)＝a₀(cosE_j-e₀)P+a₀(1-e₀ ²)^1/2sinE_jQ；

其中，P＝(cosω₀cosΩ₀-sinω₀sinΩ₀cosi₀，cosω₀csinΩ₀+sinω₀sinΩ₀cosi₀，sinω₀sini₀)，

Q＝(-sinω₀cosΩ₀-cosω₀sinΩ₀cosi₀，-sinω₀csinΩ₀+cosω₀cosΩ₀cosi₀，cosω₀sini₀)。

其中，i₀为轨道倾角，Ω₀为升交点赤经，ω₀为近地点幅角。

由上可知，通过上述的步骤11～13，即可根据卫星瞬时轨道根数计算得到卫星位置信息。在上一具体实施例中，所述卫星位置信息即为卫星在地球中心惯性系中的位置矢量R。

步骤102，从卫星遥感数据中获取临边像素点探测器的指向信息。

卫星遥感数据(例如，卫星图像数据文件)中一般都包含每一个像素点的传感器指向信息(即探测器的指向信息)，例如，风云3号卫星和MODIS卫星的一级和二级数据产品中均包含上述信息。因此，可以从卫星遥感数据中获取临边像素点探测器的指向信息。

较佳的，在本发明的具体实施例中，临边像素点探测器的指向信息包括方位角α和俯仰角β。例如，所述临边像素点探测器的指向信息即为临边像素点探测器的方向矢量L：

L＝(cosαcosβ,cosαsinβ,sinα)＝(m，n，p)。

步骤103，根据卫星位置信息和临边像素点探测器的指向信息，计算临边像素点的位置。

较佳的，在本发明的具体实施例中，可以根据如下所述的公式计算临边像素点N的位置：

(x_n，y_n，z_n)＝(x_p+mt，y_p+nt，z_p+pt)；

其中，(x_n，y_n，z_n)为临边像素点N的坐标(即位置)，t＝-(mx_p+ny_p+pz_p)。

图2为本发明实施例中的计算临边像素点的位置的示意图。如图2所示，P为卫星位置，坐标为(x_p,y_p,z_p)，O为地球质心，PN为临边像素点探测器的指向，平面基于ABCD垂直于PN，O在ABCD平面内。

因此，PN的直线方程为：(x-x_p)/m＝(y-y_p)/n＝(z-z_p)/p；

OABCD平面的直线方程为：mx+ny+pz＝0。

步骤104，计算临边像素点的切线高度。

较佳的，在本发明的具体实施例中，之前的计算均是在地球中心惯性系中进行的，而考虑到地球并非标准的球体，因此还需要将坐标转换到地心地固坐标系中进行计算，以获得临边像素点的切线高度。

较佳的，在本发明的具体实施例中，可以利用转换矩阵将地球中心惯性系中的坐标转换为地心地固坐标系中的坐标。

具体地，利用转换矩阵F＝(cosθ_g sinθ_g 0,-sinθ_g cosθ_g 0,0,0,1)，将N点在J2000.0坐标系中的坐标转换到WGS84坐标系中，其中θ_g为格林尼治真恒星时角。

N_wgs84＝(x_nwgs84，y_nwgs84，z_nwgs84)＝FN；

N点的纬度B可以迭代计算获得；

令tanB⁽⁰⁾＝z_nwgs84/(x_nwgs84 ²+y_nwgs84 ²)^1/2，

由tanB^(n＝1)＝(z_nwgs84+Re₀ ²sinB⁽ⁿ⁾)/(x_nwgs84 ²+y_nwgs84 ²)^1/2

反复迭代直到

|tanB⁽ⁿ⁺¹⁾-tanB⁽ⁿ⁾|≤σ(σ＝1e-6)为止，其中R为N点对应地球椭球的曲率半径，最终得到B。

临边像素点的切线高度：

H＝z_nwgs84/sinB–R(1-e₀ ²)

R为N点对应地球椭球的曲率半径。

综上可知，在本发明中的基于卫星图像数据获得临边像素点的地球切线高度的方法中，可以利用卫星遥感数据中卫星的位置和指向等信息，在地球坐标系中计算遥感数据中每个地球临边像素所处的切线高度。与现有技术中的方法相比，本发明中所提出的方法不受临边星观测的限制，并且考虑了地球的椭球形状，因此计算的结果更为准确。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明保护的范围之内。

Claims (9)

1.一种基于卫星图像数据估计临边像素点的地球切线高度的方法，其特征在于，该方法包括：

获取卫星位置信息；

从卫星遥感数据中获取临边像素点探测器的指向信息；

根据卫星位置信息和临边像素点探测器的指向信息，计算临边像素点的位置；

计算临边像素点的切线高度。

2.根据权利要求1所述的方法，其特征在于，所述获取卫星位置信息包括：

从预报星历或精密星历中读取所需的卫星位置信息。

3.根据权利要求1所述的方法，其特征在于，所述获取卫星位置信息包括：

根据卫星轨道根数计算探测器成像时刻对应的平近点角；

根据卫星轨道根数和探测器成像时刻对应的平近点角，计算探测器成像时刻对应的偏近点角；

根据卫星轨道根数和探测器成像时刻对应的偏近点角，计算卫星在地球中心惯性系中的位置矢量。

4.根据权利要求3所述的方法，其特征在于，所述卫星轨道根数中包括：

历元时刻t₀、轨道半长轴a₀、历元时刻t₀对应的平近点角m₀、轨道偏心率e₀、轨道倾角i₀、升交点赤经Ω₀和近地点幅角ω₀。

5.根据权利要求4所述的方法，其特征在于，根据如下的公式计算探测器成像时刻对应的平近点角m_j：

m_j＝m₀+(GM/a₀ ³)^1/2(t_j-t₀)   (1)

其中，GM为地球引力常数，t_j为探测器成像时刻。

6.根据权利要求5所述的方法，其特征在于，所述计算探测器成像时刻对应的偏近点角包括：

令E_j ⁽⁰⁾＝m_j；；

根据公式：E_j ⁽ⁿ⁺¹⁾＝E_j ⁽ⁿ⁾–(E_j ⁽ⁿ⁾-e₀*sin E_j ⁽ⁿ⁾-m_j)/(1-e₀*cos E_j ⁽ⁿ⁾)进行反复迭代，直到满足条件：

|E_j ⁽ⁿ⁺¹⁾-E_j ⁽ⁿ⁾|≤σ(σ＝1e-6)；

若n≥20仍不满足上述条件时，则取：E_j＝m_j+e₀*sin(m_j)+e₀ ²sin(2m_j)；

其中，e₀为轨道偏心率；n为迭代次数；σ为迭代阈值。

7.根据权利要求6所述的方法，其特征在于，根据如下的公式计算卫星在地球中心惯性系中的位置矢量R：

R＝(x_p，y_p，z_p)＝a₀(cosE_j-e₀)P+a₀(1-e₀ ²)^1/2sinE_jQ；

其中，P＝(cosω₀cosΩ₀-sinω₀sinΩ₀cosi₀，cosω₀csinΩ₀+sinω₀sinΩ₀cosi₀，sinω₀sini₀)，

Q＝(-sinω₀cosΩ₀-cosω₀sinΩ₀cosi₀，-sinω₀csinΩ₀+cosω₀cosΩ₀cosi₀，cosω₀sini₀)；

其中，i₀为轨道倾角，Ω₀为升交点赤经，ω₀为近地点幅角，P和Q均为单位矢量。

8.根据权利要求7所述的方法，其特征在于，

所述临边像素点探测器的指向信息为临边像素点探测器的方向矢量L：

L＝(cosαcosβ,cosαsinβ,sinα)＝(m，n，p)。

其中，α为方位角，β为俯仰角。

9.根据权利要求8所述的方法，其特征在于，根据如下的公式计算临边像素点N的位置：

(x_n，y_n，z_n)＝(x_p+mt，y_p+nt，z_p+pt)；

其中，(x_n，y_n，z_n)为临边像素点N的位置，t＝-(mx_p+ny_p+pz_p)。