《1 前言》

1 前言

海洋表面风矢量是影响海浪、海流、水团的活跃因子和海洋动力学的基本参数，在提高全球大气、海洋动力学预报模式的准确性等研究中有着重要的价值[1，2 ] 。同时，海面风矢量是影响航海、海上作业、渔业生产等的主要因素，是优化航线、航路保证、避免台风、搜索和救援工作的关键，因此，海面风矢量的观测具有重要意义。

卫星遥感以大面积的同步测量、获取速度快、覆盖范围大、时空分辨率高、可连续观测等优点而成为目前全球风矢量观测的主要手段，其中星载微波散射计以其能够在晴空和有云条件下全天候提供海面风矢量（风速和风向）观测数据等特点成为目前为止获取全球海洋风矢量最主要的微波遥感器[1] 。自 1978 年美国发射了第一个搭载业务化运行散射计的卫星SeaSat以来[3] ，已有多个星载散射计投入业务化运行，包括ERS1/2 散射计[4] 、NSCAT[5] 、QuikSCAT 与 ADEOSⅡ[6] 以及 Metop[7] 上搭载的微波散射计。

我国于20世纪80年代开始了微波散射计的研究计划，并于2002年发射的神舟四号飞船上搭载我国第一个试验用星载散射计，2011 年8月发射了海洋二号（HY-2A）卫星，搭载有我国首个可业务化运行的卫星散射计。本文将结合国外同类型微波传感器成熟的业务化风场反演算法以及HY-2A卫星散射计自身传感器的特点，建立我国自主的HY-2A 卫星业务化运行风场反演算法。

《2 HY-2A卫星散射计简介》

2 HY-2A卫星散射计简介

2011年8月发射的HY-2A卫星，搭载有我国第一个可业务化运行的卫星散射计（HY2-SCAT）。 HY2-SCAT主要用于全球海面风场观测，测风风速范围为 4～24 m/s，风速精度为 2 m/s 或 10 %；风向测量范围为 0°～360°，风向精度为±20°。HY2- SCAT 工作频率为 13.256 GHz，采用笔形波束圆锥扫描方式，通过笔形波束以固定仰角围绕天底方向旋转，在卫星平台顺轨方向的运动中形成一定的地面覆盖刈幅（见图1a）；散射计系统包括 VV 和 HH两个极化方式，分别以不同入射角进行观测，在平台的运动过程中对同一分辨单元可获取不同极化方式，不同入射角度的多次后向散射系数（σ 0 ）测量结果（见图1b），以克服海面风场方向反演的多值模糊问题。其中内波束采用HH极化方式，入射角为41°，对应地面刈幅宽度约为1 350 km。外波束采用VV极化方式，入射角为 48°，对应地面刈幅宽度约为 1 700 km。散射计的主要仪器参数如表1所示。作为对比，表1中同时给出了QuikSCAT相应的系统参数[8] 。

《图1》

图1 HY2-SCAT观测几何图

Fig.1 HY2-SCAT measurement geometry

《表1》

表1 HY2-SCAT主要系统参数与QSCAT主要系统参数对比

Table 1 Comparison of system parameters between HY2-SCAT and QSCAT

《3 风场反演方法》

3 风场反演方法

HY-2A 散射计风矢量反演算法的目标是从一组HY-2A散射计 σ 0 测量结果中估算真实风矢量。对散射计风场反演算法，目前主要有：SOS（sum of square）法、MLE（maximum—likelihood estimation）法 、LS（1east squares）法、WLS（weighted least squares）法 、AWLS （adjustable weighted least squares）法、L1- norm 法和 LWSS（1east windspeed squares）法等[9~12] 。其中MLE方法具有反演精度高、完全独立于模型函数和后向散射测量值采用自然单位、取值范围不受限制等优点，海洋二号微波散射计风场反演算法采用MLE法。MLE的原理是寻找可以使目标函数（加权 σ 0 测量值与模型预测值的残差）取得极大值的风速与风向。由于联系 σ 0 与海面风矢量的地球物理模型是高度非线性的，因此，寻找全局最优解的风矢量反演过程也是非线性的优化过程，不存在一般封闭解。即使在理想无噪声的测量条件下，由于模型函数对风向的双余弦特征，对 σ 0 的反演仍将存在多解的情况。在无噪声的情况下，如果地球物理模型对顺风、逆风的差异足够敏感，并且 σ 0 有多于两个方位角的测量结果，那么全局最优解通常为“真解”。当测量结果存在噪声时，不仅不能预知解的数量，而且全局最优解为真解的概率将远小于 1，通常情况下是 50 %左右。 HY-2A散射计采用如下形式的目标函数：

式（1）中，， J MLE为最大似然值，为风速 U 和风向 Φ 的函数； N 为风矢量单元内不同方位角/入射角 σ 0测量结果的数量； 为对应第i 个σ 0测量结果，M 为地球物理模型预测的σ 0，对应方位角为、入射角为、极化方式为观测条件下，风速为 U 和风向为 Φ 情况下的σ 0结果。海洋二号微波散射计业务化运行采用NSCAT-2地球物理模型，该地球物理模型函数采用查找表形式，即在三维参数（风速、相对方位角、入射角）的网格节点上预先计算出 σ 0 值[13] 。由于查找表为离散值，需采用查找表插值算法获得所需任意风速和方位角对应的 σ 0 。

风矢量反演算法的目的就是要找出使得目标函数取得极大值的一组 U 和 Φ 。由于模型函数高度非线性，为提高运算效率，反演算法分为粗搜索和精搜索两个步骤，首先通过粗搜索并获得初始解，这些初始解被代入精细搜索及优化算法中对其进行修正。

《3.1 粗搜索算法》

3.1 粗搜索算法

粗搜索算法以相对较粗的风速和风向搜索间隔，在风速风向空间快速搜索目标函数的局部极大值，获得初始解。由于目标函数在风速风向空间对风速呈山脊状分布，为提高算法效率，采用先搜索目标函数山脊，再从山脊中搜索局部最大值的搜索方式。同时，对风速赋予初始估计值，引导算法从对第一个风向搜索得出的风速结果附近区域开始搜索。

粗搜索算法主要计算步骤包括：按一定的风向间隔，在 0°~360°范围内，对一组给定的风向，确定最似然“山脊”与“山脊”对应的风速。对于每个风向，取三点窗口，对应风速分别为 U0 - dU ，U0 ， U0 + dU（dU 为风速搜索步长），分别计算窗口中每点风速对应的目标函数值 J1 ，J2 ，J3 。若中点的目标函数值为局部最大值（仅在速度空间考虑），则最适风速 U 及相应的目标函数值 J 可通过牛顿插值公式求出：

式（3）中

如果三点窗口中1或3中任一端点目标函数取值为最大，那么就向相应的方向移动三点窗口，并替换中点速度，同时计算新加入点的 J 值，再检测三点的目标函数值。循环这一过程，直到三点窗口中中点对应的目标函数成为局部极大。同时保存中值速度U0 ，以用于对下一风向的搜索。

粗搜索的最后一个步骤是采用三点一组的策略，在通过前面的步骤得出的最优“脊”上，在 0°~ 360°风向范围内，检测连续三点 J(φ) 中，可以使第二点（中点）取最大值的 φ 及对应目标函数。粗搜索最终将获得2～6个局部最大值。在某些情况下，由于数据噪声的影响，将会出现病态解。通常，这种情况的特征是出现一组数量远超过预期的解。目标函数值较小的解将被假定为伪解，而目标函数值最大的4～6个解将被采纳为粗搜索的解，并在后续算法中做出进一步优化。

粗搜索完成后，得到的近似解将被传递至精搜索算法，以对近似解进行优化，并根据目标函数值的大小对最终获得的解进行排序，最终完成对单个风矢量单元的风矢量反演。

《3.2 精搜索算法》

3.2 精搜索算法

为确保风矢量反演算法的精度，在粗搜索之后，有必要采用更精细的风速及风向间隔，在粗搜索近似解附近的窗口内进行精搜索。该算法以粗搜索结果为中心，在（U ，Φ）空间取3×3大小的精细网格窗口，搜索目标函数在该窗口内的最大值，并在该最大值不在窗口中心网格点时移动该窗口，并重复上述操作，直到目标函数最大值位于窗口的中心。该算法仅需少数几步操作即可从初始（U0 ， Φ0）解中得出优化的目标函数极大值。

设置窗口中心网格点索引为（0，0），其周围节点索引间隔为 1（见图 2）。对最大值 Jmax 出现在网格点（i，j）情况，（i，j）坐标可视为在各自方向上应移动的距离——“移动矢量”。对 9 点网格，表 2 列出了所有可能出现的情况，每种情况都以 CASE = | i | + | j |标识。

《图2》

图2 精搜索算法9点式二维网格分布图

Fig.2 Distribution of 9-point grid for optimal-search algorithm

《表2》

表2 精搜索可能出现的情况表

Table 2 Three possible cases in 9-point grid optimal-search algorithm

在不考虑计算效率的情况下，可以以初始解（U0 ， Φ0）对应的格点为起始点，外围点距离起始点（±dU，±dΦ）设置搜索窗口，计算搜索窗口内各网格节点的目标函数值J（i，j），同时确定最大值所在的位置。当最大值不在网格中心时，通过最大值所在的位置（i，j）计算移动矢量，将9点窗口（U ，Φ）平面移动（i ×dU，j ×dΦ），即新网格中心点位于前一步骤得到的最大值所在位置，并计算新的网格各节点对应的目标函数值，搜索最大值，循环上述操作直到 CASE = 0，此时中心点目标函数值为网格窗口 9个格点中的极大值。

在需要考虑运算效率的情况下，这种算法显得有些效率低下。第一，没有必要在每一次移动之后计算所有9个格点的目标函数值。由于目标函数值 J 的计算是风矢量反演中运算量最大的部分，因此将旧格点值传递给新格点，将有效地提高运算效率。第二，对风矢量解的优化采用5点形式的格点进行计算将比采用 9 点形式的格点更有效率。对 CASE = 1的情况，5点格式仅需要将格点移动到边缘点。对CASE = 2的情况，5点格式需要经过两步移动，计算量为 5 + 3n = 11（n = 2），而采用 9 点格式，计算量将为9 + 5 n =14（n = 1）。在5点格式搜索完成后，有必要再计算并检验角点，以在最后一个步骤中进行高阶插值。如果位于角点的值为局部最大，则从该点开始继续执行9点格式搜索。表3 列出了在特定移动条件（CASE值）下，5点格式以及 9点格式需要的计算目标函数值的操作数，从表3中可以看出，采用 5 点法较 9 点法在计算效率上有了一定提高。

《表3》

表3 计算目标函数值的操作数统计表

Table 3 Counts of MLE evaluations

《3.3 模糊解去除算法》

3.3 模糊解去除算法

3.3.1 矢量圆中数滤波算法

在风矢量反演的过程中，有多个风矢量解可使目标函数式取极大值，其中只有一个解是真实解，其余的称伪解或模糊解。所以在利用MLE方法求得使目标函数取得局部最大值的风矢量后，还要进行风向的多解去除，以得到真实解。常用的模糊解去除算法包括基于观测资料以及雷达数据的模糊解去除算法[14] ，借助空气动力学的约束条件场方式模糊解去除方法[15] 以及中数滤波技术等[16，17] 。海洋二号微波散射计业务化运行风向多解去除采用矢量圆中数滤波算法。中数滤波用无噪声的邻点数据替代误差点数据，特别适合于一个风矢量点与周围邻点方向相反时的所谓180°模糊问题，特别当用于大气风场时，中数滤波不会把大于所开窗口的低频特性如收敛线、气旋等在风向上变化剧烈的特性滤掉[16，17] 。用 ，…表示模糊风向（由似然值按顺序给出），则风向反演误差与脉冲噪声相类似。真实风向可作为信号值，用最大似然估计的（第一风场解）作为观测值，+π作为误差值，则脉冲模型可表达如下：

式（6）中， = [-1，0，1]是风向反演误差模型，表示反演过程中其他随机误差。当 =0时， 表示真实风向，当 =±1， 相应于与真实风向成180°的伪解。

风向模糊排除的目标是从 （k =1，2，3，…）中选择一个风向使得与真风向最接近，换言之，该方法通过选择下标k 使最小。真风向在运算过程中未知，但可用矢量圆中数滤波法对每一个风矢量面元上进行真风向估计。首先，通过选择真风向等于面元周围窗口内风矢量的矢量圆中数，然后从模糊解中选择出接近于真风向估计的解。最后，基于这些新选择的解，重新估计真风向值。上述估计真风向的过程连续迭代直到所选风矢量不变或迭代次数超过给定的最大次数。记表示真风向的估计值（有时称之为参考风向），则第m 次迭代得到的风矢量可表示为

式（7）中，CMF表示一个N ×N 窗上的矢量圆中数滤波算子，Wij为该面元上的权，当风矢量面元上不包含任何风矢量或面元不在刈幅上，Wij =0。

3.3.2 矢量圆中数滤波初始场算法

矢量圆中数滤波技术的物理基础是风矢量面元的风向不是独立的，而是与周围风矢量面元风向具有一定的相关性，通过周围风矢量面元的风向，计算出一个中数，然后将风矢量面元中风向与中数最接近的解赋为真值，对每个风矢量面元都做同样的操作，完成一次迭代。经过多次迭代，结果稳定之后，即得到多解的模糊性消除风矢量。关于初始解的选择，可以采用两种方法：a.以最可能的风矢量解作为初始解；b.以数值天气预报模式风场最为接近的风矢量解作为初始解。

海洋二号微波散射计初始场的选择采用第二种方式，即NWP初始场优化技术。NWP初始场优化技术是基于这样的事实：在超过85 %的情况，与真实风场最接近的解是第一或第二模糊解。此外，通过这两个解可粗略地确定风场流线（但存在方向模糊性）。由于模糊解消除的能力在很大程度上取决于仪器噪声，而非天气物理条件，因此，将物理因素加入到模糊解的选择、初始场的生成，都能够极大地提高模糊解消除的性能。将目标函数值居前两位的模糊解与第三方风向（NCEP风场）比较，初始场的性能可得到有效提高。

《4 风场反演结果及精度验证》

4 风场反演结果及精度验证

利用上述算法对HY-2A卫星微波散射计下传的遥感数据进行处理，所得结果如图3所示。图3给出了HY-2A卫星微波散射计2011年10月11日观测到的全球海面风场。该结果表明HY-2A卫星微波散射计具有全球海面风场的观测能力，1 天能够覆盖全球90 %的海域，可以捕捉到全球大部分气旋。

《图3》

图3 2011年10月11日全球风场图

Fig.3 Global winds derived from HY2-SCAT （2011-10-11）

为验证海洋二号微波散射计海面风场反演精度，分别利用NDBC浮标数据和ASCAT卫星散射计风场数据对HY-2A卫星散射计海面风场反演结果进行检验。

《4.1 利用NDBC浮标观测风场数据进行验证》

4.1 利用NDBC浮标观测风场数据进行验证

利用NDBC浮标现场观测数据对散射计风速、风向反演结果进行对比验证。选用浮标空间分布如图4所示。取时间匹配窗口为10 min，空间匹配窗口为25 km，即当散射计测量时间与浮标测量时间小于等于10 min，空间距离小于等于25 km时，则认为散射计与浮标测量结果为同步观测，同时采用2 σ质量控制方法，剔除无效数据。从 2011 年 11 月 1 日到 2011年11月30 日，共获得 1 043 组同步观测数据，统计分析结果如图5、图6所示，该结果表明HY-2A 卫星微波散射计风速和风向大小与 NDBC 浮标海面风场大小比较的均方根误差为1.5 m/s和19.5°。

《图4》

图4 同步观测的NDBC浮标位置分布图

Fig.4 Geographic location of the NDBC buoys which were used for HY2-SCAT winds validation

《图5》

图5 HY2-SCAT测量海面风速散点图（与NDBC浮标对比）

Fig.5 Comparison of HY2-SCAT wind speed with NDBC buoy measurements

《图6》

图6 HY2-SCAT测量海面风向散点图（与NDBC浮标对比）

Fig.6 Comparison of HY2-SCAT wind direction with NDBC buoy measurements

《4.2 利用ASCAT卫星观测风场数据进行验证》

4.2 利用ASCAT卫星观测风场数据进行验证

利用ASCAT卫星观测风场数据对HY-2A卫星微波散射计风速和风向反演结果进行对比验证。数据匹配时取时间窗口为2 h，空间匹配窗口为25 km，即当HY-2A卫星微波散射计测量结果与ASCAT卫星测量结果时间间隔小于等于2 h，空间距离小于等于 25 km 时，则认为 HY-2A 卫星微波散射计与 ASCAT卫星观测结果为同步观测，同时采用2σ质量控制方法剔除无效数据。从2011年11月1日到2011年 11 月 11 日，共获得 7 401 组同步观测数据，由于 ASCAT卫星过境时间分别为世界时1点左右（降轨）和13点左右（升轨），HY-2A卫星过境时间为世界时 10点左右（升轨）和22点左右（降轨），匹配点主要分布在高纬度地区（见图7）。统计分析结果如图8、图9 所示，统计分析结果表明，HY-2A散射计海面风速大小相对于ASCAT卫星散射计观测值的均方根误差为 1.5 m/s，HY-2A 散射计海面风向大小相对于 ASCAT卫星散射计观测值的均方根误差为17.5°。

《图7》

图7 HY2-SCAT与ASCAT观测匹配点空间分布图

Fig.7 Geographic location of the colocation point between HY2-SCAT measurements and ASCAT measurements

《图8》

图8 微波散射计测量海面风速散点图（与ASCAT观测海面风场对比）

Fig.8 Comparison of HY2-SCAT wind speed with ASCAT measurements

《图9》

图9 微波散射计测量海面风向散点图（与ASCAT观测海面风场对比）

Fig.9 Comparison of HY2-SCAT wind direction with ASCAT measurements

《5 结语》

5 结语

海洋二号微波散射计是我国第一个可业务化运行的星载微波散射计，它的成功使我国卫星海洋风场观测水平与国际先进水平缩短了30年。海洋二号微波散射计可提供风速精度为2 m/s或10 %，风向精度为±20°，刈幅宽度为1 700 km的海面风场观测产品。在 2009 年 11 月 QuikSCAT 卫星散射计失效以后，海洋二号微波散射计成为目前唯一可提供25 km×25 km分辨率，且每天覆盖全球90 %以上海域的微波散射计，对保证全球海洋风场的连续性具有重要意义。目前HY-2A微波散射计已在全球海面风场观测、海上风暴监测、海洋环境预报等方面发挥了重要作用。此外，HY-2A微波散射计可以覆盖从南纬88°~北纬88°的区域，在极地航行保障、极地海冰变化等方面具有重要应用前景。从2011年 9月28日开机获取第一轨海面风场数据以来，海洋二号微波散射计已累计获得了大量宝贵的海面风场观测数据，目前正处于其3年正常服役寿命的前期，随着数据的积累和算法的改进，海洋二号微波散射计必将为海洋观测和应用做出更大的贡献。