CN113960642A - 青藏高原区域光学遥感卫星影像大气折射定位误差修正方法 - Google Patents

Abstract

本发明涉及卫星遥感技术。本发明的目的是提供一种青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，首先，根据青藏高原区域历史地面观测气象资料，对青藏高原区域大气的温度、气压及水汽压强参数按大气高度变化进行数值拟合，从而建立青藏高原区域球形大气模型；其次，根据预设层差对大气进行分层，并根据拟合后的数值计算每一层的大气折射率；然后，计算初始高度的折射角，并以此依次向下逐层计算出最底层大气的折射角；并根据最底层大气折射角计算出地面大气折射误差；再根据地面大气折射误差修正卫星影像定位。能适用于青藏高原区域的特殊气候，能根据大气高度变化，计算不同层高的大气折射率，从而精准确定卫星定位误差，并进行修正。

Description

青藏高原区域光学遥感卫星影像大气折射定位误差修正方法

技术领域

本发明涉及卫星遥感技术，特别涉及光学遥感卫星影像定位误差修正技术。

背景技术

卫星遥感技术在政治、经济、军事和社会等众多领域都有极大的应用价值，是目前各国争相发展的热门技术。遥感卫星平台在地球大气层以外，光线传播透过地球大气层时，由于空气密度变化，在地表附近产生弯曲，称为大气折射效应。特别是在卫星侧视角较大的情况下，这种折射对光学遥感卫星定位产生较大影响。由于地球大气层具有不均匀性，其密度，分布情况随时间、空间、所处地理位置等诸多因素的影响而变化，该误差具有极大的不确定性。

光学遥感卫星影像的几何定位基于摄影测量中的共线方程原理，通过确定影像的内、外方位元素来恢复摄影光线，建立严密的影像几何模型来计算影像上每个像点的物方空间位置。大气折射会改变地物反射光线在星地间的直线传播方向，破坏地物点、成像载荷透视中心、CCD传感器三者的共线条件，造成基于共线方程的严密几何定位结果存在一定的大气折射误差。随着目前光学超高分辨率卫星技术的发展，分辨率0.5m以下卫星越来越多，在侧视角较大的情况下，大气折射误差远远超过1个像素，必须进行纠正。

现有技术中，主要是把地球大气简化为单层球形大气的单层模型(Noerdlinger1999；Oh and Lee 2011)，以及把地球大气简化为对流层和同温层的双层模型两种(Yan etal.2016)。认为在单层模型和双层模型中，同层大气具有相同的折射率。根据斯涅尔定律(Snell's Law)，从卫星到地面计算观测光线的折射角，最后计算大气折射误差。

首先，大气折射率主要受到大气密度变化的影响，而大气密度又与大气的温度、压强、水汽参数等密切相关。同时，这些参数往往是和大气层高度密切相关的。目前的计算方法，只是简单的把大气按单层，或双层看待，并且假设在同一层中，大气折射率是固定的。但是，我们注意到，在光线传播过程中，大气高度的变化是连续的，因此，大气折射率的变化也应该是连续的。而目前的计算方法并不能反映出这种特点，存在一定的模型误差。

其次，目前的计算方法，一般是按中纬度标准大气的气象条件对大气进行建模。而在我国青藏高原地区(介于北纬26°00′～39°47′，东经73°19′～104°47′之间)，由于其特殊的气候和大气环境，大气状态与标准大气模式存在较明显的差异。如仍沿用标准大气模型来计算青藏高原地区的光学遥感影像大气折射误差，从理论上来说很难取得令人满意的结果。

发明内容

本发明的目的是要提供青藏高原区域的光学遥感卫星影像大气折射定位误差修正方法，使其能适用于青藏高原区域的特殊气候，能根据大气高度变化，计算不同层高的大气折射率，从而精准确定卫星定位误差，并对误差进行修正。

本发明解决上述技术问题采用的一种技术方案是，青藏高原区域光学遥感卫星影像大气折射定位误差修正方法,包括以下步骤：

步骤1、根据青藏高原区域历史地面观测气象资料，对青藏高原区域大气的温度、气压及水汽压强参数按大气高度变化进行数值拟合，从而建立青藏高原区域球形大气模型；

步骤2、根据预设层差对大气进行分层，并根据拟合后的数值计算每一层的大气折射率；

步骤3、计算初始高度的折射角，并以此依次向下逐层计算出最底层大气的折射角；

步骤4、根据最底层大气折射角计算出地面大气折射误差；

步骤5、根据地面大气折射误差修正卫星影像定位。

具体的，所述步骤1中，对青藏高原区域大气的温度、气压及水汽压强参数按大气高度变化进行数值拟合，具体如下：

当大气高度H在0-11(km)时：

T(K)＝-5.85473×H+305.549；

P(Pa)＝9740×exp(-((H+41.29)/27.35)²)；

e(Pa)＝194.2×exp(-((H+11.88)/8.654)²)；

当大气高度H在11-17(km)时：

T(K)＝-6.317868×H+308.321；

P(Pa)＝2501×exp(-((H+21.3)/21.3)²)；

e(Pa)＝1183×exp(-0.8605×H)；

当大气高度H在17-30(km)时：

T(K)＝2.54×H+156.4；

P(Pa)＝1985×exp(-0.2286×H)+537.7×exp(-0.1313×H)；

e(Pa)＝-0.000001378×H³+0.00009408×H²-0.002092×H+0.01566；

当大气高度H在30-100(km)时：

T(K)＝0.00008621×H⁴-0.01916×H³+1.429×H²-41.91×H+651.7；

P(Pa)＝3.956420×[214.65/(214.65-2×(h-71))]^(34.1632/-2)；

e(Pa)＝0.1393×exp(-0.1968×H)；

其中，T表示温度、H表示大气高度、P表示大气压强，e表示水汽压强。

具体的，所述步骤2中，任意一层的大气折射率的计算公式如下：

其中，

T表示温度，H表示大气高度，P表示大气压强，e表示水汽压强，均可从数值拟合表中得知，其中，N_ph为群折射度、n_ph为大气折射率、N_sph为相折射度、n_sph为相折射率、λ为光学波长。

具体的，所述步骤2中，预设层差为100米。

具体的，所述步骤3中，包括以下步骤：

步骤31、假设从卫星传感器位置，发射一束光线，在高度H₁达到距离卫星0.5Δh处产生折射，侧视角为θ，则入射光线单位向量为：

入射光线所在的直线方程为：

y₁＝kx₁+b；

其中，k＝-tanθ,b＝R+H_s,H_s表示卫星高度；

步骤32、建立以地球球心为圆心，R₁为半径的圆方程：

x₁ ²+y₁ ²＝R₁ ²；

其中，R₁＝R+H_s-d_h，d_h表示相邻大气分层的高度间隔,R表示地球半径；

联立式，可得：

步骤33、计算入射光线的折射角，令向量

则单位向量为

向量

与向量

的夹角为θ₁，则有

则：

步骤34、计算下一层入射光线的折射角，根据折射定律，有：

sinθ₁×n₁＝sinθ₂×n₂；

其中，n₁与n₂分别表示第一层大气折射率与第二层的大气折射率；

步骤35、计算下一层大气的折射角，将θ₂作为下一层入射角，重复以上步骤，依次计算每一层大气的折射角；当光线与地球表面相交时，计算出最底层大气的折射角。

进一步的，所述步骤4中，根据最底层大气折射角计算出地面大气折射误差，具体包括：

假设，P(X_P,Y_P)为光线到达地面的实际位置，P_n(X_pn,Y_pn)为不受大气折射影响，光线沿直线路径与地面相交的位置；大气折射在地面产生的误差为d，地球近似为半径为R的圆球，建立二维直角坐标系，从地球球心到卫星为y轴，天顶方向为正，与y轴垂直方向为x轴，以东方向为正，N为星下点，则

对应的圆心角为β₀，

对应的圆心角为β_n，则大气折射误差大小为：

d＝R×Δα＝R×(β₀-β_n)。

再进一步的，所述步骤5中，设影像上P点在地心地固坐标系，坐标系中的坐标为(x,y,z)，卫星方位角为α，根据地面大气折射误差修正卫星影像定位计算公式如下：

本发明的有益效果是，在本发明方案中，根据青藏高原区域历史地面观测气象资料，对青藏高原区域大气的温度、气压以及水汽压强参数按大气高度变化进行数值拟合，从而建立青藏高原区域球形大气模型，可根据拟合数值实时计算出不同层高的大气折射率。按100m高程间距，对大气层进行连续分层。由于100m的高差变化幅度小，可认为在100米的竖直剖面范围内，大气折射率相同。通过分层计算不同层高的大气折射率可以更加精准的确定入射光线在每一层的折射角，并最终计算出光线与地面的折射角，从而精准确定卫星定位的误差，并对其进行精准修正。

附图说明

图1是本发明实施例1中多层大气中的光线矢量折射示意图。

图2是本发明实施例1中影像坐标值的折射误差改正图。

图3为本发明实施例2中不同侧视角下的大气折射误差结果。

图4为本发明实施例3中的影像定位误差图。

具体实施方式

下面结合实施例及附图，详细描述本发明的技术方案。

本发明所述的青藏高原区域的光学遥感卫星影像定位误差修正方法，首先，根据青藏高原区域历史地面观测气象资料，对青藏高原区域大气的温度、气压以及水汽压强参数按大气高度变化进行数值拟合，从而建立青藏高原区域球形大气模型；其次，根据预设层差对大气进行分层，并根据拟合后的数值计算每一层的大气折射率；然后，计算初始高度的折射角，并以此依次向下逐层计算出最底层大气的折射角；并根据最底层大气折射角计算出地面大气折射误差；再根据地面大气折射误差修正卫星影像定位。现有技术中，一般是按中纬度标准大气的气象条件对大气进行建模。而在我国青藏高原地区，由于其特殊的气候和大气环境，大气状态与标准大气模式存在较明显的差异。如仍沿用标准大气模型来计算青藏高原地区的光学遥感影像大气折射误差，从理论上来说很难取得令人满意的结果。本发明中，提出一种适用于青藏高原区域的，根据大气高度变化的多层大气折射误差精确计算方法。本方法在大气建模上考虑青藏高原特殊的气候和大气环境的特点，结合实际气象资料与标准模型进行联合建模；同时，在大气折射误差几何模型构建上，本发明按固定高差对大气层进行连续分层，更能反映大气折射率的连续变化，理论上更严密。通过本发明的应用，可以直接计算影像点在地面的大气折射改正值，从而直接提高光学遥感卫星在青藏高原地区的影像定位精度。

实施例1

本实施例主要包括两个部分：第一部分，建立青藏高原区域的温度、压强、水汽压强随高度变化的拟合模型；第二部分，建立基于大气连续分层的光线矢量折射误差严密几何改正模型。

第一部分，青藏高原区域大气温度、压强随高度变化的拟合模型：

由于青藏高原在地形、气候方面存在较明显的区域特征，大气状态与标准大气模式存在较明显的差异，如仍沿用标准大气模型来计算青藏高原地区的光学遥感影像大气折射误差，从理论上来说很难取得令人满意的结果。建立区域大气模型来计算大气折射误差，提高卫星影像定位精度，是非常有必要的。本例根据青藏高原地区11个地面观测站十年的地面观测气象资料，并参考中纬度与热带大气模式，对青藏高原区域大气的温度、气压以及水汽压强参数按大气高度变化进行数值拟合，得到如下公式：

表1青藏高原区域温度、气压及水汽压强数值拟合表

大气密度计算采用公式：

ρ＝P/(R_Specific×T)(kg/m³) (1)

R_Specific是气体常数，大小为R_Specific＝287.053J·kg^-1·K^-1。

本例在可见光波段和近红外波段，在指定标准大气环境为温度T＝273.15K，大气压强P＝1013.25hPa，CO₂含量x＝0.0375％，大气群折射度计算公式为：

e为大气中水汽压强，单位hPa。

其中,大气相折射度N_sph计算公式为：

T表示温度，H表示大气高度，P表示大气压强，e表示水汽压强，均可从数值拟合表中得知，其中，N_ph为群折射度、n_ph为大气折射率、N_sph为相折射度、n_sph为相折射率、λ为光学波长。

本例采用公式(2)结合拟合数值，计算每一层大气折射系数。

第二部分，基于大气连续分层的光线矢量折射误差严密几何改正模型：

本例按100m高程间距，对大气层进行连续分层。由于100m的高差变化幅度小，可认为在100米的竖直剖面范围内，大气折射率相同。

本例计算大气折光的方法，如图1所示，当光线进入大气层时，随着空气密度逐渐增大，大气折射系数也越来越大，光线传播路径变成弯曲的连续折线。P(X_P,Y_P)为光线到达地面的实际位置，P_n(X_pn,Y_pn)为不受大气折射影响，光线沿直线路径与地面相交的位置。可以看出，大气折射在地面产生的误差为d,d是地球表面的一段弧线。地球近似为半径为R的圆球，建立二维直角坐标系，从地球球心到卫星S为y轴，天顶方向为正，与y轴垂直方向为x轴，以东方向为正，N为星下点，则

对应的圆心角为β₀，

对应的圆心角为β_n，则大气折射误差大小为：

d＝R×Δα＝R×(β₀-β_n) (4)

其中,如图1所示，

X坐标通过光线矢量与地球球面相交来确定，球面定义为半径为R+H_i的圆球面，H_i为光线发生折射时的高程面，则H_i

H_i＝H_S-0.5×Δh (5)

Δh为相邻高程折射面之间的高差，本例按每100米高差计算大气折射率。在100米高差范围内，可认为大气折射率不变。

光线在多层大气中的传播过程如图1所示，多层大气折射的严密几何改正模型具体计算技术流程为：

1.确定入射光线的直线方程

假设从卫星传感器位置，发射一束光线，在高度H₁达到距离卫星0.5Δh处产生折射，侧视角为θ，则入射光线单位向量为：

入射光线所在的直线方程为:

y₁＝kx₁+b； (7)

其中，k＝-tanθ,b＝R+H_s,H_s表示卫星高度。

2.求解入射光线与同心球面的交点位置

建立以地球球心为圆心，R₁为半径的圆方程为半径的圆方程

x₁ ²+y₁ ²＝R₁ ² ； (8)

其中，R₁＝R+H_s-d_h，d_h表示相邻大气分层的高度间隔,R表示地球半径。

联立式(7)、(8)，解得：

3.计算入射光线的折射角

令向量

则单位向量为

向量

与向量

的夹角为θ₁，则有

则：

4.计算入射光线的折射角

根据折射定律(斯涅尔定律)，有

sinθ₁×n₁＝sinθ₂×n₂ (10)

5.高程降低Δh，计算下一个高程面的折射角。

把θ₂作为下一个入射角，重复步骤1-4，依次计算每个高程面的折射角θ_i。

6.当光线与地球表面相交时，计算最终的折射角。

重复1-4步，直到H＝0为止，最后得到θ_n，带入公式(4)，计算出地面大气折射误差d。最后，对高分辨率卫星影像定位进行修正。

7.设设影像上P点在地心地固坐标系(Earth-Centered,Earth-Fixed，简称ECEF)坐标系中的坐标为(x,y,z)，卫星方位角为α，则大气折射改正值对计算如图2所示，计算公式为公式12：

实施例2

为了确定本发明的实用性与优越性，本例针对不同侧视角的大气折射误差计算。

首先，在卫星轨道高度为650km，地球平均半径为6371km条件下，以中心波长为0.5μm的蓝光计算大气折射系数，使用式(1)-式(10)以5°间隔计算15°至45°侧视角范围内大气折射产生的几何定位偏差，并与文献(Noerdlinger 1999)单层大气，文献(Ming etal.2016)双层大气计算的大气折射几何偏差进行比较，结果如图3所示。

从图3实验结果，可以分析如下：

1、在侧视角小于30°时，采用本发明算法，用标准大气模型的计算结果与参考文献中单层法，及双层法基本一致。

2、在侧视角大于30°时，三种方法计算结果的差值逐渐变大，本发明算法结果更加准确。双层法计算的大气折射误差最大，单层法最小，而本发明方法居于两者之间。分析原因如下，本发明从卫星传感器视线出发，按大气连续分层计算计算大气折射系数及其对视线的折射偏差；而双层法也是从卫星传感器视线出发进行计算，按照加权平均算法分别计算对流层和同温层折射系数。实际上，在对流层和同温层内部，大气折射系数的变化并不是线性的。加权平均法并不能准确反映这种变化。特别是在高度大于11km时，对流层顶层过渡到平流层，此时大气折射系数加速变小。双层法过度估计了大气折射误差，对于单层法，是把对流层和同温层统一为单层大气模型来计算大气折射偏差。这种假设又低估了大气折射系数的实际变化。而本发明的算法，是按每100米的高度变化值对大气层进行分层，连续计算每层大气折射系数和对应的折射角，从光线传播的几何原理上看，本发明是一种严密的计算方法，计算结果更加精确。

实施例3

本例采用青藏高原地区高分辨率卫星影像实验。

本例对本发明的实际使用效果，采用四幅青藏高原区域的WorldView-2卫星0.5米分辨率全色影像进行试验，影像覆盖区域位于东经95°41′17.19″～96°01′45.24″，北纬36°07′46.54″～36°19′57.60″，行政区划隶属于青海都兰县境内，地处柴达木盆地南缘，昆仑山中段北坡，地貌上属深切割高山区，海拔高程在3300米—4500米间，相对高差在300—800米左右。影像的基本信息如表3：

表2青藏高原地区高分辨率卫星影像实验基本信息

Image ID	Ⅰ	Ⅱ	Ⅲ	Ⅳ
					行×列	27,164×34,584	24,356×35,180	22,716×35,180	28,048×34,445
扫描方向	Forward	Forward	Forward	Forward
					平均地面分辨率(m)	0.516	0.553	0.649	0.518
卫星方位角(°)	350.000	210.500	40.300	97.400
					卫星侧视角(°)	17.800	24.500	33.100	17.600

在上述影像覆盖区域，人工非均匀布设了80个地面检查点，采用GPS技术在野外进行了实地测量，检查点平面精度和高程精度优于5cm。

按四种方案进行影像定位精度实验，第一种采用影像提供商的原始有理函数模型进行计算，没有进行大气折射改正；第二种采用单层法；第三种采用双层法；第四种采用本发明计算方法。

根据遥感影像定位模型，计算检查点的模型计算位置P(X,Y,Z)与已知坐标P₀(X₀,Y₀,Z₀)的差值，e表示在WGS84空间直角坐标系下的点位误差，则e计算公式如下：

根据每一幅影像中每个检查点的误差，计算出误差均值及中误差。

从图4中可以看出，对于实验区影像，应用本发明的方法影像定位精度最高。其中，实验区的四幅影像定位精度平均提高了18％；而对于目前的计算方法，本发明计算精度平均提高了11％。目前的结果仅仅是侧视角小于32°的情况下获得的，随着侧视角的增大，本发明计算方法上的优势将会更加明显。

综上，光学高分辨率卫星遥感是对地观测技术的重要组成部分，在国民经济社会多个方面具有重要的应用价值，而大气折射误差对光学遥感卫星定位精度产生较大影响。现有技术中的计算方法对大气层分层模型描述过于简单，忽视了大气折射率连续变化的特点，本发明采用每100米的高差变化计算大气折射率和折射角，更能精确的描述观测光线的实际传播路径；同时，针对青藏高原区域的特殊气象条件，采用区域气象模型进行大气折光系数的计算，计算结果更准确，提高了光学高分辨率卫星影像在该区域的影像实际定位精度。

Claims (7)

1.青藏高原区域光学遥感卫星影像大气折射定位误差修正方法,其特征在于，包括以下步骤：

步骤1、根据青藏高原区域历史地面观测气象资料，对青藏高原区域大气的温度、气压及水汽压强参数按大气高度变化进行数值拟合，从而建立青藏高原区域球形大气模型；

步骤2、根据预设层差对大气进行分层，并根据拟合后的数值计算每一层的大气折射率；

步骤3、计算初始高度的折射角，并以此依次向下逐层计算出最底层大气的折射角；

步骤4、根据最底层大气折射角计算出地面大气折射误差；

步骤5、根据地面大气折射误差修正卫星影像定位。

2.根据权利要求1所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤1中，对青藏高原区域大气的温度、气压及水汽压强参数按大气高度变化进行数值拟合，具体如下：

当大气高度H在0-11(km)时：

T(K)＝-5.85473×H+305.549；

P(Pa)＝9740×exp(-((H+41.29)/27.35)²)；

e(Pa)＝194.2×exp(-((H+11.88)/8.654)²)；

当大气高度H在11-17(km)时：

T(K)＝-6.317868×H+308.321；

P(Pa)＝2501×exp(-((H+21.3)/21.3)²)；

e(Pa)＝1183×exp(-0.8605×H)；

当大气高度H在17-30(km)时：

T(K)＝2.54×H+156.4；

P(Pa)＝1985×exp(-0.2286×H)+537.7×exp(-0.1313×H)；

e(Pa)＝-0.000001378×H³+0.00009408×H²-0.002092×H+0.01566；

当大气高度H在30-100(km)时：

T(K)＝0.00008621×H⁴-0.01916×H³+1.429×H²-41.91×H+651.7；

P(Pa)＝3.956420×[214.65/(214.65-2×(h-71))]^(34.1632/-2)；

e(Pa)＝0.1393×exp(-0.1968×H)；

其中，T表示温度、H表示大气高度、P表示大气压强，e表示水汽压强。

3.根据权利要求2所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤2中，任意一层的大气折射率的计算公式如下：

其中，

T表示温度，H表示大气高度，P表示大气压强，e表示水汽压强，均可从数值拟合表中得知，其中，N_ph为群折射度、n_ph为大气折射率、N_sph为相折射度、n_sph为相折射率、λ为光学波长。

4.根据权利要求1或2或3所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤2中，预设层差为100米。

5.根据权利要求4所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤3中，包括以下步骤：

步骤31、假设从卫星传感器位置，发射一束光线，在高度H₁达到距离卫星0.5Δh处产生折射，侧视角为θ，则入射光线单位向量为：

入射光线所在的直线方程为：

y₁＝kx₁+b；

其中，k＝-tanθ,b＝R+H_s,H_s表示卫星高度；

步骤32、建立以地球球心为圆心，R₁为半径的圆方程：

x₁ ²+y₁ ²＝R₁ ²；

其中，R₁＝R+H_s-d_h，d_h表示相邻大气分层的高度间隔,R表示地球半径；

联立式，可得：

步骤33、计算入射光线的折射角，令向量

则单位向量为

向量

与向量

的夹角为θ₁，则有

则：

步骤34、计算下一层入射光线的折射角，根据折射定律，有：

sinθ₁×n₁＝sinθ₂×n₂；

其中，n₁与n₂分别表示第一层大气折射率与第二层的大气折射率；

步骤35、计算下一层大气的折射角，将θ₂作为下一层入射角，重复以上步骤，依次计算每一层大气的折射角；当光线与地球表面相交时，计算出最底层大气的折射角。

6.根据权利要求5所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤4中，根据最底层大气折射角计算出地面大气折射误差，具体包括：

假设，P(X_P,Y_P)为光线到达地面的实际位置，P_n(X_pn,Y_pn)为不受大气折射影响，光线沿直线路径与地面相交的位置；大气折射在地面产生的误差为d，地球近似为半径为R的圆球，建立二维直角坐标系，从地球球心到卫星为y轴，天顶方向为正，与y轴垂直方向为x轴，以东方向为正，N为星下点，则

对应的圆心角为β₀，

对应的圆心角为β_n，则大气折射误差大小为：

d＝R×Δα＝R×(β₀-β_n)。

7.根据权利要求6所述的青藏高原区域光学遥感卫星影像大气折射定位误差修正方法，其特征在于，所述步骤5中，设影像上P点在地心地固坐标系，坐标系中的坐标为(x,y,z)，卫星方位角为α，根据地面大气折射误差修正卫星影像定位计算公式如下：

Images