《1 前言》
1 前言
海洋二号(HY-2)海洋动力环境卫星是一颗获取海洋动力环境信息的微波遥感卫星,其主要任务是监测和调查海洋动力环境,包括海面风场、海浪、海流、海面温度、海上风暴和潮汐等,为海洋环境监测与预报、海洋调查与资源开发、海洋污染监测与保护、海洋权益维护、军事应用、全球变化和海洋科学研究等方面提供遥感数据[1] 。HY-2卫星有效载荷包括雷达高度计、微波辐射计、微波散射计、校正辐射计和数据传输分系统。此外,HY-2卫星还具有精密定轨功能,需要在地面对定轨数据进行处理和应用。
微波辐射计(RM)天线安装在卫星载荷舱中板上,如图1所示,其电轴指向与卫星+Z 轴呈40º。因此,RM坐标系原点为RM安装面中心,其-Y 轴方向为从安装面中心指向精测镜安装方向,-Z 轴垂直安装面指向天线,X 轴与Y 轴、Z 轴呈右手系。名义上,卫星主动段时RM各坐标轴方向与整星坐标系方向一致。RM探测头部转动的0º位置定义为卫星飞行方向(+X 向),顺时针为角度增加方向,有效观测角度为+X 方向的±70°范围。
《图1》
图1 RM对地观测示意图
Fig. 1 Schematic diagram of RM earth observation
RM主要技术指标要求如表1所示。
《表1》
表1 RM主要技术指标
Table 1 RM main technical indicators
注:V—垂直极化;H—水平极化
RM 分系统具有两种基本工作模式:正常工作模式和应急工作模式。正常工作模式是指RM在轨期间处于持续工作状态,观测天线进行圆锥扫描,旋转一周约3.6 s,依次观测地球、冷空和定标热源,并形成4个子包数据。RM扫描一周系统工作模式示意图如图2所示,+X 方向为0°,每周观测的70°~ 290°(标称值,具体数值可从遥感包中读出,下同)为对地观测范围,152°~168°为冷空观测范围,211°~ 235°为热源观测范围。获取的遥感数据和定标数据组包后下传。应急工作模式是指RM出现故障或整星需要时,关闭部分或全部设备。
《图2》
图2 RM扫描一周系统工作模式示意图
Fig. 2 Schematic diagram of RM scanning work mode
《2 RM海洋参数反演算法》
2 RM海洋参数反演算法
《2.1 RM海洋遥感原理》
2.1 RM海洋遥感原理
在晴空、有云和不大于2 mm/h的轻度降雨的情况下,在6~37 GHz RM频段,当不考虑大雨滴和云颗粒散射影响的时候,以地表为下边界、以冷空背景为上边界的大气辐射传输模式可以用吸收-辐射模式近似。Wentz等[1,2] 对SSM/I的降雨观测结果研究表明,雨率大于2 mm/h的情况仅为所有海洋降雨的 3 %,因此吸收-发射模式将可以应用于 97 %的海洋辐射计观测。
对于被动微波遥感来说,最重要的量是单色辐照度 L (θ,ϕ ) ,单位为Js-1 m-2 sr-1 Hz-1 ,其定义为:在球坐标系下,沿 (θ,ϕ ) 方向传播的单位频率、单位立体角的辐射通量密度。此外,还要考虑辐射的极化方式。通常把信号分解成为水平极化和垂直极化两部分,分别用H和V表示,在以后的表达式中统一用p 表示。根据Beer定律,沿路径ds 的辐射衰减(消光)与沿路径的质量密度ρ 成正比,即
式(1)中, = + , 是质量消光系数, 是质量吸收系数, 是质量散射系数,当使用体消光、吸收和散射系数时,分别乘以密度即可。与辐射消光相对应的是辐射的热发射和散射。根据Kirchhoff定律,物质的热辐射和吸收系数与Planck函数B(t)成正比
式(2)中,t 是热力学温度,c 是光速,ν 为频率,h 是Planck常量,k 是Boltzmann常量。由于散射引起的辐射增强可以表示为单次散射反照率 ω= ks /ke 和辐照度L对立体角的加权平均值的乘积,权为r(θ,ϕ ;θ′, ϕ ′)。通过以上定义,完整的辐射传输方程(RTE)可以表示为
对于微波辐射,应用Rayleigh-Jeans近似,L可以用亮温T B 表示
同时,使用吸收-辐射近似,不考虑散射影响,ω=0 , = ,RTE简化变为
由于 = ,故用 α 代替e,表示大气吸收系数。在大多数情况下(降雨除外),考虑微波在大气中的辐射传输时,仅考虑平面平行大气情况就已足够,也就是说, α 和T 都是高度h 的函数(忽略地球曲率)。这时,ds 可用 dh/ μ 代替, μ=cosθ ,则下行亮温T BD可表示为
式(6)中,T BC 是宇宙背景亮温,为 2.7 K; 是大气透过率。同理,大气顶部观测到的上行亮温T BU 的贡献为
从空间观测到的亮温包括式(6)、式(7)的贡献,还包括地球表面反射及发射成分。最简单的情况是平面反射,表达式为
式(9)中,ep 是在极化方式p下的平面发射率;TS 是表面热力学温度。然而,在大多数情况下,表面反射不是一个严格的平面发射,而是漫反射,因而
式(10)中, rp (μ,ϕ;μ′,ϕ′ ) 是双站散射系数,表示下行 (θ,ϕ ) 方向的辐射亮温散射到向上的卫星辐射计 方向 (θ′,ϕ′) 的角度权重函数。 其中,μ′=cosθ′ 。漫反射表面发射率 ep 定义为
可见,无雨情况下,RM从空间观测到的亮温不但取决于大气的温度廓线T (h )和吸收系数廓线 α (h ),还取决于表面温度 TS 和双站反射系数 rp( μ,ϕ;μ′,ϕ′ ) 。其中,α (h )是与观测频率相关的量, rp与观测的极化方式和海面状况,尤其是海面风场有关。当这些参数能够通过模型模拟后,就可以通过积分计算获得观测亮温的模拟结果。
《2.2 亮温的模拟方法和模拟结果》
2.2 亮温的模拟方法和模拟结果
2.2.1 海洋和大气模式
为了模拟RM不同通道的亮温,使用欧洲数值预报中心(ECMWF)的 FASTEM3 模式作为海洋发射率和反射率模式,其考虑了大尺度和小尺度粗糙表面的发射率贡献,同时考虑了各风向对微波辐射信号的调制,以及大气稳定性对海面散射的影响。大气模式分别考虑了氮气、氧气、水汽和液态水的贡献。本文分别采用MPM93中的氧气和氮气模式,使用PWR98的水汽模式以及PWR04的液态水模式[3] 。
2.2.2 海面和大气数据集
使用了EUMETSAT/ECMWF提供的全球60层大气廓线,其中海面廓线6 813条,本文选取了表面温度大于 270 K 的 6 154 条廓线用于海面亮温模拟。这些廓线包括大气的温度、压力、水汽密度、液态水密度、云盖等廓线,同时包括海面温度、风速和风向等信息。另外,为了分析降雨的影响,还将液态水含量大于0.2 mm的廓线去掉,得到3 962条无雨大气廓线,单独进行了亮温模拟。
2.2.3 不同通道亮温的模拟结果
图3~图6分别给出在47.7°入射角下对6.6 GHz 和37 GHz垂直极化和水平极化的全球亮温的模拟结果(其他通道类似,由于篇幅关系这里省略)。其中分别给出了包括降水在内的所有点(6 813 条廓线)和无雨数据点(3 962 条廓线)的模拟结果的比较。由图3~图6可以看出,当降雨出现时,亮温有所升高,特别是37 GHz,亮温升高比较明显。
《图3》
图3 6.6 GHz垂直极化亮温分布
Fig. 3 Brightness temperature at 6.6 GHz vertical polarization
《图4》
图4 6.6 GHz水平极化亮温分布
Fig. 4 Brightness temperature at 6.6 GHz horizontal polarization
《图5》
图5 37 GHz垂直极化亮温
Fig. 5 Brightness temperature at 37 GHz vertical polarization
《图6》
图6 37 GHz水平极化亮温分布
Fig. 6 Brightness temperature at 37 GHz horizontal polarization
《2.3 反演算法的建立》
2.3 反演算法的建立
星载RM接收的亮温包括两部分:一部分是没有泡沫的海面贡献,另一部分是泡沫覆盖区的辐射贡献。对式(10)简化,得到卫星高度接收的垂直极化亮温TBV和水平极化亮温TBH为[4]
式(12)和式(13)中,Ep 、Rp (p=V,H)表示海面发射率和反射率。假设泡沫的覆盖率为f,泡沫的发射率为1,那么
式(14)和式(15)中, 是与风向无关的表面发射率, 、 表示风向引起的发射率变化。
式(16)和式(17)中, 是平静海面发射率; 、 是风成海面发射率。
对于反射率,为了简化可表示为
(18)和式(19)中, Ω V 、 Ω H 分别表示对粗糙表面引起反射率变化的校正系数。如果风向引起的发射率变化表示为下面的形式
那么,星载微波辐射计接收的垂直和水平极化亮温中的风向信号( ΔTV 和 ΔTH )分别为
因此,反演算法就要建立平静海面发射率 、风成海面发射率 的反演算法,另外还要确定校正系数 、不同频率的系数 。其中p,q=V,H表示极化方式。同时为了反演大气参数,如水汽含量V 和液态水含量L,还要建立这些大气参数与大气上行亮温 T BU 、下行 亮温 T BD 、大气透过率 之间的关系。
2.3.1 海面算法
1) 平静海面反射率算法。海面可以近似看成空气和海水两种介质的作用面。因此,确定入射辐射被海面反射的问题就变成以一定的入射角和极化方式入射的平面电磁波与边界面作用的问题。
反射率R 是入射角 θ 和极化方式的函数,同时也是空气的复介电常数 ε 和海水的复介电常数 ε′λ 的函数。由 Fresnel 反射定律可以计算出反射系数 R0。根据Ellison等[5] 提出的介电常数模式模拟海面温度 范围为-2~32 ℃、盐度 s 范围为 29 ‰~37 ‰,利用Fresnel反射定律,计算平静海面反射率[4] ,这里从略。
2) 风成海面发射率算法。海面泡沫覆盖率采用Monahan等[6] 提出的方法计算,即
的算法形式为
表 2 给出了 5 个频率下水平和垂直极化通道的系数。
《表2》
表2 入射角为47.7º的风成海面发射率系数
Table 2 Emissivity coefficient of windy sea surface at 47.7º
海面反射率的校正系数 Ωp ,对于垂直极化和水平极化,表示为如下形式
式(26)中,A表示大气光学厚度, A=-lnτ secθ ; σ 表示海面坡度方差,是风速的线性函数,而且随着频率ν (这里单位为 GHz)的变化而变化,其表达式为 [7]
式(27)中, W 10 表示10 m高度风速,m/s,以后如果不特殊说明,风速都指10 m高度测量的风速,用 W 表示。式(26)的系数见表 3,生成系数的数据来自RTTOV8.5软件的反射率系数文件 [8] 。
《表3》
表3 反射率校正算法的系数
Table 3 Coefficients of reflectivity correction
3) 风向算法。风向引起发射率的变化采用下面的形式确定系数 的表达式
表4给出式(28)和式(29)中不同频率和极化下的多项式系数。这些系数是根据RTTOV8.5软件模拟的结果采用多元回归分析产生的。表4最后一列给出拟合出现的最大偏差的绝对值。可见这项引起的海面亮温变化可以忽略。
《表4》
表4 风向算法系数
Table 4 Coefficients for wind direction
2.3.2 无雨和有雨两种情况下的大气算法
如果假设垂直路径(θ=0º)时下行亮温和上行亮温分别为T BD0 和T BU0 ,那么任何其他路径上的T BD 和T BU 都可以用T BD0 和T BU0 的线性方程表示为
式(30)和式(31)中的系数、b、c 和 d 可以表示为secθ 的二阶多项式形式(见表5和表6),即(以为例,b、c 和d 形式相同)
《表5》
表5 TBD 线性方程系数的二阶多项式系数
Table 5 Coefficients for TBD linear equation
《表6》
表6 TBU 线性方程系数的二阶多项式系数
Table 6 Coefficients for TBU linear equation
表5和表6中给出算法的残差的定义为拟合值和测量值之差的平方和的算术平方根。TBD0 和TBU0之间的关系同样可以用直线方程近似描述为
式(33)中,在入射角分别为0º和47.7º下的系数m 和n 见表7和表8。
《表7》
表7 无雨大气上行和下行辐射亮温之间的关系系数
Table 7 Coefficients for TBD and TBU equation (without rain)
《表8》
表8 有雨情况下在47.7°入射角时大气上行和下行辐射亮温之间的关系系数
Table 8 Coefficients forTBD and TBU equation (with rain,incidence angle is 47.7°)
由表7和表8可见,23.8 GHz和37 GHz的残差虽然大一些,但是由于离散大的点主要体现在大气的异常辐射区,如高层水汽区域,所以对于常规大气条件,这样拟合的离散误差基本可以忽略。从表 7和表8也可看出,TBD0 和TBU0 非常接近。利用直线拟合可以用一个参数替代另一个参数。但是通过对无雨算法的残差分析可以发现,残差随着水汽含量的增大而发生变化,具有明显的二次曲线的关系,因此把这个关系加入到无雨模型中去,则有
表9给出了系数及残差。
《表9》
表9 在0°入射角时大气上行和下行辐射亮温之间的关系系数
Table 9 Coefficients for TBD and TBU equation (incidence angle is 0°)
垂直路径大气透过率 可以用大气水汽含量 V(mm)和液态水含量L(mm)表示为
通过对无雨和有雨的大气剖面分别进行模拟的结果发现,这个线性方程能够较好地拟合透过率。表 10和表11给出式(35)中的系数、b、c。那么在任何路径下,大气透过率 可以表示为
《表10》
表10 无雨情况下0°入射角观测大气透过率算法的系数
Table 10 Coefficients for atmospheric transmissivity (without rain,incidence angle is 0°)
《表11》
表11 有雨情况下0°入射角观测大气透过率算法的系数
Table 11 Coefficients for atmospheric transmissivity (with rain,incidence angle is 0°)
上文给出了利用大气透过率 和下行亮温TBD表示的大气算法,为了达到反演水汽含量V和液态水含量L的目的,笔者建立了一个简单的大气透过率 、下行亮温TBD 和大气参数V、L 的反演算法。根据模拟的大气下行亮温TBD 和大气参数V、L 直接进行回归分析,获得如下形式的算法
式(37)中,系数~e 在不同频率、47.7º入射角下的拟合结果见表 12~表 15。可以看出,算法的均方差(RMS)随着频率增大有增大的趋势。
《表12》
表12 无雨情况下47.7º入射角观测下行大气亮温算法的系数
Table 12 Coefficients for TBD (without rain,incidence angle is 47.7°)
《表13》
表13 无雨情况下47.7º入射角观测上行大气亮温算法的系数
Table 13 Coefficients for TBU (without rain,incidence angle is 47.7°)
《表14》
表14 有雨情况下47.7º入射角观测下行大气亮温算法的系数
Table 14 Coefficients for TBD (with rain,incidence angle is 47.7°)
《表15》
表15 有雨情况下47.7º入射角观测上行大气亮温算法的系数
Table 15 Coefficients for T BU (with rain,incidence angle is 47.7°)
2.3.3 算法误差分析
前面的海洋算法和大气算法是对理论微波辐射传输方程的近似,目的是用来反演所关心的海洋和大气参数。对于海洋算法,输入的参数包括海面温度(SST)、风速(WS)、风向 ϕ 和海水盐度 s ;大气算法的输入参数包括大气的水汽含量V 和液态水含量L。大气的上行亮温TBU 、下行亮温TBD 和透过率 都是这两个大气参数的函数。但是由于它们之间的函数关系非常复杂,这里只给出了简单的线性模式,所以如果用V 和L 作为风场反演的输入参数,必然会引入额外的误差。因此首先是用TBU 、TBD 、 作为大气算法的输入参数。另外由于大气的TBU 和TBD之间的差值很小,存在很好的线性关系,因此用TBD代替 TBU ,这样大气的输入参数就变为两个:TBD 、 。表16给出了用大气算法的两种输入方式和4个海洋参数(其中盐度 s 设定为35 ‰)算法计算的亮温与理论计算的亮温之间的RMS(单位为K)。可见用V 和L 算法引入的误差随着频率增高而增大,并且垂直极化的误差小于水平极化。而用TBD 和 作为大气参数的算法误差很小,基本上可以完全忽略。
《表16》
表16 使用不同大气参数反演算法的RMS
Table 16 RMS by using different atmospheric retrieval algorithm
K
《3 海洋和大气参数反演方法》
3 海洋和大气参数反演方法
《3.1 海面参数SST、WS、V和L 反演方法》
3.1 海面参数SST、WS、V和L 反演方法
根据第2.3节的算法,海面参数SST、WS、V 和L 可以表示为
式(38)中, ci 是系数; Fi (i=1,2,3,...,9)分别对应于 6.6V,6.6H,10.7V,10.7H,18.7V,18.7H,23.8V,37V,37H这9个通道的亮温数据,定义为对于23.8V(i=7)以外的通道, Fi =TB ch_i -150 ,对于23.8V(i=7)通道, F 7 =-ln(290-TB ch_7 ) 。
式(39)中, ci 是系数; Fi (i=1,2,3,...,9)分别对应于 6.6V,6.6H,10.7V,10.7H,18.7V,18.7H,23.8V,37V,37H这9个通道的亮温数据,定义为对于6.6 GHz和10.7 GHz的4个通道(i=1,2,3,4), Fi =TB ch_i ,此外的通道(i=5,6,7,8,9), Fi =-ln(290-TB ch_i ) 。
式(40)中, ci 是系数; Fi (i=1,2,3,...,9)分别对应于 6.6V,6.6H,10.7V,10.7H,18.7V,18.7H,23.8V,37V,37H 这 9 个通道的亮温数据,定义为对于6.6 GHz 和 10.7 GHz 的 4 个通道(i=1,2,3,4),Fi =TB ch_i ,此外 的通 道 (i=5,6,7,8,9),Fi =-ln(290-TB ch_i ) 。
表17和表18分别给出了无雨和有雨条件下的反演系数。无雨条件下,SST、WS、V 和 L 的算法RMS分别为0.88 K,0.41 m/s,0.61 mm和0.02 mm;而有雨条件下的算法RMS分别为10.9 K,0.55 m/s,0.85 mm和0.05 mm。
《表17》
表17 无雨条件下算法系数表
Table 17 Algorithm coefficients without rain
《表18》
表18 有雨条件下算法系数表
Table 18 Algorithm coefficients with rain
《3.2 海洋和大气中间参数的反演》
3.2 海洋和大气中间参数的反演
根据式(6)的算法模式进行18.7 GHz和37 GHz的大气透过率 和上行亮温TBU 的反演。算法的输入参数是海面的发射率和反射率,其通过前面的反演结果和算法模型计算得到。此外,为了进行最小二乘法估值,需要输入 TBD19 、 TBD37 、 和 初始值,作为待反演这4个参数的输入条件。这4个参数的初始值通过的大气算法结合 V 和 L 的反演值产生。
利用算法结合实测亮温,用最小二乘法反演TBD19 、 TBD37 、 和 。
《3.3 液态水贡献》
3.3 液态水贡献
在应用第2.3节进行反演的时候,反演的参数包括TBD19 、 TBD37 、 和 。为了进行降雨的估算,需要利用Wentz的AMSR-E算法剥离大气氧气和水汽的贡献,得到大气液态水贡献的测量值,这样才能进行反演。因此,其中的关于大气氧气和水汽透过率的算法是必需的模型。当SST、WS和V确定以后,根据Wentz [1] 的AO 和AV 模式获得 和 ,进一步获得液态水的透过率 、 ,以及各自的上行大气辐射亮温。
式(41)和式(42)中的系数见表19。
《表19》
表19 算法系数
Table 19 Coefficient
因此,
对应于18.7 GHz和37 GHz分别有
《3.4 降雨判断方法》
3.4 降雨判断方法
如果 ≥0.8 ,则可认为大气无降雨,雨率(RR)=0,利用无雨算法系数反演海面参数。如果 <0.8 ,开始降雨反演。或者等效地认为,如果0.04≤A L37 ≤0.15 则表示降雨出现;若 A L37 =0.04 则表示有毛毛雨出现;若 A L37 =0.15 则表示雨率大约为2 mm/h,则利用有雨算法系数反演海面参数。
《4 海洋参数反演结果和验证》
4 海洋参数反演结果和验证
《4.1 海洋参数反演结果》
4.1 海洋参数反演结果
图7~图10给出了用2011年12月31日RM升轨数据反演的全球参数。其结果表明,算法合理可行,能够反演出大气海洋参数。
《图7》
图7 升轨海面温度的反演结果
Fig. 7 Retrieved SST (ascending)
Zhao [9] 和Huang等 [10] 根据第2节中的算法,利用RM数据,对海面温度和海面风速进行了反演,并与其他星载辐射计产品及海面浮标产品进行了比对验证。
《图8》
图8 升轨海面风速的反演结果
Fig. 8 Retrieved WS (ascending)
《图9》
图9 升轨海面水汽含量的反演结果
Fig. 9 Retrieved V (ascending)
《图10》
图10 升轨海面液态水含量的反演结果
Fig. 10 Retrieved L (ascending)
《4.2 海面温度反演结果及验证》
4.2 海面温度反演结果及验证
11 给出了 RM 反演的海面温度和 NDBC 现场测量结果及 ARGO 浮标结果的比对。数据匹配的空间尺度为25 km,时间窗口为30 min。图11中的红点在计算偏差和标准偏差时排除在外,红线表示线性拟合线,而蓝线表示置信区间。两个验证结果都表明海面温度的初始反演精度优于1.7 K[9] 。
《图11》
图11 RM反演的海面温度与NDBC及ARGO 数据比对结果[9]
Fig. 11 Scatter plots of RM versus NDBC mooring and ARGO SST [9]
图12给出了2012年6月,RM反演的海面温度与 WindSat 反演结果平均偏差,网格精度 0.25º× 0.25º。结果表明,在2012年上半年,RM与WindSat 反演结果的RMS为1.8 K[9] 。
《图12》
图12 2012年6月,全球RM与WindSat反演的海面温度偏差[9]
Fig. 12 The mean bias between RM/HY-2 and WindSat SST in June 2012 [9]
《4.3 海面风速反演结果及验证》
4.3 海面风速反演结果及验证
图 13 给出了 RM 反演的海面风速与 NDBC 浮标测量结果的比较,空间尺度为25 km,时间窗口为 5 min。其结果显示,RMS 为 1.89 m/s,偏差为 0.67 m/s[10] 。
《图13》
图13 RM反演的风速与NDBC风速比对结果[10]
Fig. 13 Scatter plots of comparison between HY-2 RM WS observations and buoys measurements [10]
图14给出了2012年1月—6月RM反演结果与 WindSat 反演结果平均偏差全球分布,网格大小为 0.25º×0.25º。从图14中可以看出,风速差异出现在纬度较高沿海地区,特别是欧洲和北非地区差异较大,接近2 m/s;而在低纬度热带海洋,差异在-0.5~ 0.5 m/s;在中纬度温带地区,在1~2 m/s。两者比较的偏差和RMS分别为0.23m/s和1.8m/s[10] 。
《图14》
图14 2012年1月—6月RM和WindSat风速反演结果差异全球分布[10]
Fig. 14 Global map of averaged WS difference between HY-2 RM and WindSat gridded data between January 2012 and June 2012 [10]
《5 结语》
5 结语
1)本文探讨了用RM亮温数据反演海洋大气参数的方法原理,建立了反演的算法,并且根据RM在轨运行的数据进行了这些参数的初步反演。反演结果表明,算法合理,能够反演海面和大气参数。
2)分析了 RM 反演的海面温度和海面风速与 NDBC、AGRO浮标及WindSat产品的比对结果。结果显示,海面温度与 NDBC 及 AGRO 浮标比对, RMS优于1.7 K,而与WindSat比对,RMS为1.8 K;海面风速与NDBC浮标比较,RMS为1.89 m/s,而与 WindSat比较,RMS为1.85 m/s。
目前,其他参数的反演产品的检验还在进行。由于RM反演算法是基于模拟数据建立的,存在一定偏差,反演算法还要根据对辐射计测量值定标的深入而进一步优化。
《致谢》
致谢
感谢黄骁骑,赵毅力提供 RM 反演结果与 NDBC及WindSat产品的比对结果。