《1 前言》

1 前言

HY-2 卫星是一种获取海洋动力环境信息的专用的对地遥感卫星[1],其主要使命是监测和调查海洋环境,获得包括海面风场、浪高、海流、温度、海上风暴和潮汐等海况的多种重要动力参数,是海洋防灾减灾的重要监测手段,可直接为灾害性海况预警报,为国民经济建设和国防建设服务,为海洋科学研究提供实测数据。

HY-2 卫星的原始数据由海洋卫星地面站负责接收,经解包等预处理后向数据处理中心提供 0 级数据产品和其他辅助产品(轨道和姿态文件等)。数据处理中心在接收到 0 级数据产品后,将其制作成 1 级数据产品,并进行产品存档与分发服务;同时负责在 1 级产品的基础上制作成 2 级数据产品,并进行产品存档与分发服务。

笔者开展了微波散射计的风场反演技术研究,并利用现有的微波散射计数据模拟建立了 HY-2 卫星微波散射计风场反演算法和模型,利用神舟 4 号(SZ - 4)飞船多模态微波遥感数据[2]和国外卫星散射计数据分别进行了应用验证。

《2 HY-2 卫星平台与微波散射计》

2 HY-2 卫星平台与微波散射计

HY-2 卫星装载有微波散射计、雷达高度计和微波辐射计。微波散射计的频率为 13.256 GHz;雷达高度计的频率为 13.58 GHz 和 5.25 GHz;微波辐射计为 5 个频率 9 个通道的扫描辐射计,其中 5 个频率为 6.6,10.7,18.7,23.8 GHz 和 37 GHz,其中除了 23.8 GHz 是 V 极化外,其余频率都是双极化。HY-2 卫星的研制基于 ZY 卫星成熟平台,将继承ZY 系列卫星的大部分设计,并在卫星总体方面进行优化调整。 卫星重量小于 1500 kg,计划利用 CZ - 4B 火箭发射上天,卫星工作设计寿命为 3 年。卫星轨道为太阳同步近圆形轨道,高度为 963 km 和 965 km 。卫 星数 传系 统 载频 选 为 X 波 段, 并有 120 Gbits 的记录/回放数据记录器。卫星轨道设计为晨昏轨道,卫星发射状态包络尺寸为 3.35 m × 4.55 m , 符合 CZ - 4B 火 箭可 提供 的 运载 空间。 HY-2 卫星平台主要技术指标见表 1 。

《表1》

表1 HY-2 卫星平台主要技术指标

Table1 The specifications of HY-2 satellite’s plat

微波散射计作为海洋动力环境卫星的重要有效载荷,主要完成获取海面风场的测量等功能。微波散射计的 σ° 测量精度要求为 0.5 dB,采用双点波束圆锥扫描体制,线性调频脉冲提高空间分辨率,支持内定标技术。HY-2 卫星微波散射计的主要技术指标见表 2 。

《表2》

表2 HY-2 卫星微波散射计主要技术指标

Table2 The specifications of the scatterometer on HY-2 satellite

《3 散射计反演风矢量场算法与模型》

3 散射计反演风矢量场算法与模型

《3.1 风矢量场反演》

3.1 风矢量场反演

微波散射计反演海面风矢量场分为两个主要步骤。a. 利用风矢量与海面后向散射系数的经验关系(即地球物理模型函数),获得多个风矢量解;b. 风矢量多解的模糊性消除。地球物理模型函数描述了海面风矢量(风速和风向)与雷达后向散射系数之间的关系。风矢量反演算法主要是通过地球物理模式函数以及海面风矢量单元的不同方位角的观测获得海面的风矢量解。由于地球物理模式函数本身的特征以及散射计的各种测量噪声,一般会获得多个风矢量解。风矢量多解消除就是从一系列的多解风矢量中选择唯一的风矢量解。地球物理模型函数的一般形式为[3~5]

式(1) 中, σ代表散射计测量的后向散射系数;w 为风速; χ 为风向的相对方位角; f 为散射计的工作频率; p 为极化方式; θ 为天线的入射角。

海面每一个风矢量单元,相互独立的后向散射系数测量次数以及它们之间的相对方位角变化是由雷达参数决定的。每一个风矢量单元有两个以上不同方位角的后向散射系数数据,理论上可以唯一的确定风速和风向,但由于模式函数呈双调和函数的特点,因此可以有多达 4 个甚至 4 个以上的风矢量解,这就是风向模糊性问题。由于地球物理模式函数对于风向是非线性关系,而对于风速是准线性关系,因此对于两个独立的观测,风速基本是唯一的。此外,每一个风矢量单元后向散射系数测量的独立程度是由相互之间的方位角差异决定的,最理想的方位角差异为 90° 。如果每一个风矢量单元有超过 3 个相互独立的后向散射系数测量数据,在噪声比较小的情况下,可以确定一个风矢量解。如果其中一个后向散射系数测量采用另外一种极化方式,那么模糊性大大降低。

微波散射计测量的后向散射系数主要是风速、风向、方位角、极化方式和入射角的函数。在入射角、风向、方位角一定的情况下,模型函数给出的后向散射系数随着风速的增加而增加。对于给定的风速,模式函数所描述的后向散射系数与风向、方位角之间呈双调和函数的特点。风向对后向散射系数的调制使得卫星散射计可以测量风向信息,在固定的风速条件下,后向散射系数在逆风观测时最大,顺风其次,而横风最小。 逆风和顺风对后向散射系数的调制以及多个不同天线方位角和入射角的后向散射系数原则上可以唯一确定风矢量信息。此外,风向对后向散射系数的调制在不同极化方式条件下有所不同。笔者采用最大似然估计作为目标函数进行优化反演,设目标函数的形式为[6,7]

式(2)中, σoi 为卫星观测后向散射系数的测量值;σm为后向散射系数的模型计算出来的后向散射系数;N 表示后向散射系数的独立测量次数; + βσ + γ =  为散射测量方差,α,β,γ 为常数;负号表示目标函数的局部极大值所对应的风速和风向为可能的风矢量解。反演流程如图 1 所示。

《图1》

图1 散射计风矢量场反演流程图

Fig.1 The flow chat about scatterometer’s wind vector field retrieved

《3.2 地球物理模式函数》

3.2 地球物理模式函数

式(2)中的 σm  模型函数形式有很多种,由于 HY-2 卫星与目前在轨运行的 SeaWind 散射计的频率很接近,笔者采用 SeaWind 模型函数[4,5]

式(3)中, w 为风速; θ 为波束入射角; p 为极化方式; χ 为相对方位角。由此可见,在其他参数不变时,后向散射系数随相对方位角的变化,以 2 为周期,且关于  对称。

《3.3 风向多解消除》

3.3 风向多解消除

利用上述目标函数求极值的方法可以得到多个海面风矢量解, 为了得到唯一的风矢量解,需要采用风向多解消除算法[8 ~10]。风向多解消除算法利用矢量中值滤波技术从一系列的多解风矢量中选择唯一的风矢量解。

矢量中值滤波的滑动窗口大小为 N × NN 为奇数),窗口中心的网格点坐标为( i j ),其多解风矢量为 Vmn 为窗口内网格点(mn)上最大似然估计值所对应的最可能风矢量解,即 Vij = h = INT( / 2 ) -1, Wm'n' 为窗口权重函数( m' m - i' = n - j ),表示网格点(ij)上第 k 个风矢量解所对应的最大似然估计值,p 为权重系数。对于每一个窗口计算窗口中心( ij ) 的滤波函数值  ,用最小的 所对应的风矢量解   代替 ,这样重复计算滑动窗口,直到   不再变化。 这里滑动窗口的大小为 5 × 5 或 7 × 7,选择目标函数 J 最小值对应的第一风场矢量解作为初始场。如图 2 所示为风向模糊排除的流程图。

《图2》

图2 圆中数滤波流程图

Fig.2 The flow chat about circle median filter

《3.4 利用神经网络训练 HY-2 卫星微波散射计物理模型》

3.4 利用神经网络训练 HY-2 卫星微波散射计物理模型

文章同时开展了基于神经网络获得地球物理模型函数的方法以及反演算法,该方法可直接应用于我国未来的 HY-2 卫星微波散射计风场反演,为我国未来的 HY-2 卫星微波散射计的地球物理模式函数研究打下一定的基础。

利用微波散射计测量的后向散射系数与风速、风向以及入射角等参数的经验关系,这里将获得风场反演模式函数的神经网络结构设计(见图 3)。网络输入层具有 5 个神经元,分别为风速,风向方位角的正、余弦值以及入射角的正、余弦值。输出层为后向散射系数,训练数据集采用利用时间—空间三维插值算法获得的 SeaWind 散射计和欧洲中期天气预报模式风场的配对数据集。 时间—空间三维插值算法的表达式如下[11]

式(5)中, , 

《图3》

图3 地球物理模式函数的神经网络结构

Fig.3 The neural network framework of geophysical model

另外,SPD_Met 为采用上式插值算法计算的与卫星轨道数据对应的模式风场数据, di  为权重系数,与模式风场数据和卫星数据的经纬度有关。time1,time0 为模式风场数据的时间;PRD_time 为卫星数据的测量时间;Zon _Meti ( time1 ),Zon _ Meti ( time0 )为 time1,time0 时刻欧洲中期天气预报模式风场的水平分量; Mer _ Meti ( time1), Mer _ Meti(time0)为 time1,time0 时刻欧洲中期天气预报模式风场的垂直分量;Met_Lati ,Met_Loni 为插值点 i 的纬度和经度,Mes_Lat,Mes_Lon 为测量点的纬度和经度。基于上述方法,我们可以获得约 40 万对 Sea - Wind 散射计和欧洲中期天气预报模式风场匹配的数据集。

这里将利用神经网络为 SeaWind 散射计的每一根天线都获得一个相应的模式函数,然后综合所有的模式函数分别获得 SeaWind 散射计 VV 极化和 HH 极化 的地 球 物 理 模 式 函 数 NN - GMF 和 HH - GMF。图 4 表示神经网络输出的后向散射系数与 SeaWind 散射计测量的后向散射系数的对比,相应的均值偏差、标准偏差以及相关系数如图中的左上角所示。图 5 表示不同的风速条件下 NN - GMF 描述的后向散射系数与 SeaWind 模型函数的对比。图 6 表示不同的风速条件下 HH - GMF 描述的后向散射系数与 SeaWind 模式函数的对比。 由图可以看出,NN - GMF 和 HH - GMF 所描述的后向散射系数与 SeaWind 模式函数的结果基本一致,并且与风向相对方位角呈简谐函数的特点,在逆风(0°) 和顺风(180°)观测时,后向散射系数达到极大值,在顺风(90°或 270°)观测时为极小值。

《图4》

图4 SeaWind 测量的后向散射系数与神经网络获得的后向散射系数的比较

Fig.4 The comparing results about the backscatter coefficient between the SeaWind’s measurement and NN computation’s

《图5》

图5 VV 极化的地球物理模型函数(σ)

Fig.5 The geophysical model function (σ) about VV polarize

《图6》

图6 VV 极化地球物理模型数据(σ)

Fig.6 The Geophysical Model Function (σ) about HH polarize

《4 风矢量场产品误差分析》

4 风矢量场产品误差分析

利用 QUIKSCAT 卫星的 SeaWind 散射计的数据根据 HY-2 卫星的参数转换成模拟的 HY-2 卫星原始散射计数据,作为文章的数据反演原始后向散射系数和对应的几何关系。

本节采用实测浮标的海面风矢量数据对卫星遥感反演得到的风矢量进行验证。同步测量数据集所采用的现场实测数据为 TAO 浮 标 实 测 海 面 的 10 min 平均风矢量数据,其站位分布如图 7 所示。对比的散射计模拟数据取自美国 RSS(remote sensing system)制作的 0.25°× 0.25° 面元的数据,取时间窗口匹配为 5 min 内,距离窗口为 0.25° 面元与 TAO 浮标同步测量数据集。对 2006 年共配对获得了 23252 组无雨条件下的同步测量数据,散射计反演的风速对 TAO 浮标风速散点图如图 8 所示,风向对 TAO 浮标风向散点图如图 9 所示。从图中可以看出,散射计反演所得风矢量与 TAO 浮标观测资料具有很好的一致性,其中 SeaWind 反演所得风速均方根误差为 1.1 m/s,风向均方根误差为 17.6° 。图 10 与图 11 则分别表示了散射计风速误差与风向误差的空间分布。从图中可以看出,卫星遥感风速反演误差东部要高于西部,特别是在 165°E 线上,误差相对较大,这很有可能是由于东部的大气与洋流活动较西部地区更强所造成的。

《图7》

图7 TAO 浮标分布

Fig.7 The bouy distribution of TAO

《图8》

图8 散射计反演风速对比散点图

Fig.8 The wind speed scatter diagram about satterometer retrieved result

《图9》

图9 散射计反演风向散点图

Fig.9 The wind direction scatter diagram about satterometer retrieved result

《图10》

图10 风速误差空间分布(单位:m/s)

Fig.10 The spatial distribution about the wind speed error (Unit:m/s)

《图11》

图11 风向误差空间分布(单位:(°))

Fig.11 The spatial distribution about the wind direction error(Unit:(°))

采用与浮标验证类似的方法,可构建 HY-2 卫星模拟数据反演与船测风矢量的同步观测数据集,对 HY-2 卫星业务化运行反演算法所得风矢量进行验证。所采用的船测资料为 908 专项调查获得的风矢量资料。908 专项调查测量的站位图如图 12 所示,站点覆盖了几乎整个中国近海区域。在 2006 夏季航次与冬季航次两个航次,对船测资料与 HY-2卫星散射计(对应 SeaWind 散射计数据转换)共配对获得了 291 个同步观测数据,其中散射计风速与船测风速散点图如图 13 所示。统计分析得出风速均方根误差为 1.8 m/s 。由此看出,HY-2 卫星散射计风速与船测风速具有很好的一致性,其均方根误差满足业务化运行所要求的小于 2 m/s 的要求。由于 908 资料未能提供船的航向等信息,与船测的风向结果未经验证。

《图12》

图12 908 专项调查站位图

Fig.12 The investigate stations of 908 project

《图13》

图13 SeaWind 风速对908 专项调查船测风速散点图

Fig.13 The wind speed scatter diagram about satterometer retrieved result

《5 微波散射计数据处理、信息反演与产品制作技术流程》

5 微波散射计数据处理、信息反演与产品制作技术流程

《5.1 数据预处理流程》

5.1 数据预处理流程

HY-2 卫星微波散射计的反演与制作主要技术流程为利用 HY-2 卫星的 0 级数据产品经过地理定位(L1A)、海陆标识、物理量转换以及面元匹配(L1B)等预处理流程获得带有定位信息和相关物理量信息的 1 级数据产品,微波散射计 1 级数据产品的处理流程如图 14 所示。

《图14》

图14 HY-2 卫星微波散射计的 1 级产品处理流程

Fig.14 The production process flow chat about HY-2 satellite scatterometer

5.1.1 地理定位

微波散射计 1 级数据产品的地理定位主要是通过卫星在地心惯性坐标系的位置矢量以及姿态确定散射计面元在地球表面的位置(经纬度)。对于微波散射计的每一个测量面元(见图 15),假设面元中心在地心惯性坐标系的位置矢量是 Rr Xr YrZr ),此刻卫星位置矢量 RSXsYs Zs ) 可以通过传统的轨道报、哥达德轨道理论计算获得(下标 r 表示地球表面面元,S 表示卫星点)。

《图15》

图15 微波散射计的地面测量面元示意图

Fig.15 The sketch map about the microwave scatterometer’s ground measurement cell

式(6) 中, R 对应着卫星到面元中心的距离,ddxdydz )是卫星与面元中心连线的 3 个余弦分量。对于卫星的速度矢量 VVSXVSY VSZ )和 Vs 的单位矢量 VVsx Vsy , Vsz ),则有以下关系存在,其中 Rxs yszs )是卫星位置矢量的单位矢量。

扫描方向的单位矢量 d 相当于卫星的地心矢量 Rs 绕卫星的飞行方向矢量旋转角度 p(散射计天线的下视角),再绕卫星的地心矢量旋转角度 α,用以下关系式表示:

其中矩阵  是从星体坐标系向惯性坐标系转换的矩阵,由于面元中心位于地球表面,所以有以下关系式存在:

或:

式(12) 中,Ae 是地球赤道半径;Be 是地球极地半径。解此方程得:

式(13)中,A = ,

=  ,

C

代入公式(7),(8),(9),(12),可以得到地面面元中心在惯性坐标系中的坐标( XrYrZr ),经惯性坐标系到地心固连坐标系的转换,得到地面面元中心在地心固连坐标系的坐标(xyz ),经过地心固连坐标系到大地坐标系的转换公式,可以得到地面面元中心的地理纬度和经度。

纬度:Lat = arctan (14 -1)

经度:Lon = arctan  (14 -2)

式(14)中,e = 

5.1.2 海陆标识

利用详细的陆地地图(按经纬度匹配)确定散射计 0 级产品的海陆标识;输入为经纬度,输出为标志位,陆地标志为 1,海洋标识为 0 。散射计 0 级产品中被确定为陆地的区域不参与以后的数据处理。

5.1.3 物理量转换

HY-2 卫星数据的一个标准长度的数据包中含有 13 次回波信号采样值,每一次回波信号幅度电平采样 22 个点。物理量转换的技术流程为:

1) 先把 22 次的测量 22 个字节(每次一个字节的测量)的测量值转化为十进制电压值,转换的公式为:0X00 对应 0 V 电压,0X80 对应 5 V 电压,取值范围为 0 ~ 0X7f;

2) 然后对 22 个测量值进行平均,求出回波平均电压;

3) 计算 AGC(自动增益控制)值(在探测模式下,而在校准模式下是恒定的),该一个字节的每个 bit 意义从高位到地位为:20,10,10,6,5,4,3,2,故计算 AGC 为:K1 = bit7 × 20 + bit6 × 10 + bit5 × 10 + bit4 × 6 + bit3 × 5 + bit2 × 4 + bit1 × 3 + bit0 × 2 ; 

4) 散射系数的计算公式为:

式(15)中,CAL 为定标系数; CON 为系统常数,为散射计系统提供;VAR 为测量变量。分别由以下的式子给出:

式(16)中,K1C 为外校准时探测状态接收机增益;K2C 为外校准时内校准状态接收机增益;Vcd 为外校准时内校准状态接收机输出电压;Vc 为外校准时探测状态接收机输出电压;σc 为标准散射体散射截面积;Rc 为标准散射体到天线距离。

式(18)中,Vi 为探测时探测状态接收机输出电压;Vid 为探测时内校准状态接收机输出电压;K1 为测量时探测状态接收机增益;K2 为测量时校准状态接收机增益。

5.1.4 微波散射计的面元配准

经过上述两个关键步骤,已经将微波散射计的 0 级数据产品处理成包括地理信息等在内的按照时间序列排序的海面后向散射系数数据。微波散射计的面元配准主要目的是将时间序列排序的散射计测量的后向散射系数数据按照地理位置配准到风矢量单元(见图 16),使得每一个风矢量单元具有多个后向散射系数数据,为风矢量单元的风场反演作数据匹配对的标准 1B 产品。

《图16》

图16 圆锥扫描的散射计风矢量单元示意图

Fig.16 The sketch map about the vector wind cell of cone scan scatterometer

《5.2 微波散射计 2 级数据产品处理流程》

5.2 微波散射计 2 级数据产品处理流程

HY-2 卫星散射计的 2 级数据产品是将 1 级数据产品经过数据质量控制后,利用模型函数获得每一个风矢量单元的多解风场,然后通过消多解算法获得唯一的风矢量解,其技术流程如图 17 所示。

《图17》

注:ECMWF 为欧洲中尺度天气预报,NCEP 为美国环境预报中心

图17 HY-2 卫星散射计风场反演技术框图

Fig.17 The sketch map about the wind retrieve technical of the HY-2 satellite’s scatterometer

《6 结语》

6 结语

针对我国将要发射的 HY-2 卫星微波散射计,进行了微波散射计的数据预处理和处理技术研究,提出了 HY-2 卫星散射计反演海面风矢量场的模型算法,给出 HY-2 卫星微波散射计数据获取、预处理、风矢量场反演和产品制作技术流程,并进行了误差精度分析,初步的验证结果表明其反演的风速和风向满足 HY-2 卫星的工程技术要求,文章所建立的反演算法和技术流程未来可以用于 HY-2 卫星的地面应用系统工程建设。

致谢:文章采用的 SeaWind(QuikSCAT) L2A 数据由 NASA JPL 提供,QuikSCAT 业务化全球风矢量数据由 RSS 提供,908 现场风数据由 908 水体调查项目组提供,在此一并表示感谢。