海面风场是海洋上层运动的主要动力来源, 与海洋中多数的物理过程密切相关, 在海洋动力学中, 它是形成海面波浪的直接动力, 也是区域和全球海洋环流的重要动力来源。因此, 快速、准确地获取海面风场信息, 有助于对海洋环境乃至全球气候变化做出更科学的判断和预报。

目前的海面风场信息观测方式有岸基观测站、船只、浮标、散射计及合成孔径雷达。其中, 岸基观测站、船只、浮标的观测范围有限且是单点分布, 维护成本较高, 无法实现大范围观测; 随着卫星遥感技术的发展, 散射计被应用到海洋风场探测领域, 可以在短时间内获取大面积的海面风场数据, 是比较成熟的风场信息观测方式, 但其空间分辨率较低且近海岸风场反演的精度较差^[1]。随着合成孔径雷达技术的发展, 星载合成孔径雷达可以实现全天候、全天时观测, 具有较高的空间分辨率(可达数米至数十米), 弥补了散射计的不足, 合成孔径雷达是目前乃至未来实现海面风场大范围观测的主要手段^[2]。

地球物理模型函数(geophysical model function, GMF)是比较成熟且使用较多的提取SAR风场信息的方法, 首先被应用于散射计。散射计利用不同的方位角对同一海面单元进行多次观测, 根据最大似然估计建立基于GMF的代价函数, 获取可能的风速、风向组合, 筛选出最优解^[3]。合成孔径雷达的观测方位角是固定的, 无法在短时间内对同一海面单元进行多次观测, 所以无法使用与散射计相同的方法从SAR数据中获取风场。但是, 如果将风向作为初始化信息, 则地球物理模型函数中的入射角、相对风向(风向与雷达观测方向的夹角)确定, 从而可以确定风速^[4]。风向的获取方法至少有3种: 1) 利用SAR图像中的线性纹理信息获取海面风向; 2) 利用数值预报模式获得风向; 3) 通过散射计获取风向信息。

早期大量的观测和研究表明, 大气边界层涡旋的轴线方向与海面风矢量的方向基本一致, 在SAR图像上表现为具有线性纹理特征的风条纹^[5]。对图像进行二维傅里叶变换可以得到风向, 但存在180°模糊。并且, 在获得的海面SAR图像中, 具有风条纹信息的SAR图像大约占60%^[6]。

Monaldo^[7]利用海军实用全球大气预测系统(NOGAPS)与SeaWinds散射计的风向作为初始风向提取SAR数据的风速, 前者反演风速与SeaWinds风速相比的均方根误差(root mean square error, RMSE)为1.78 m/s, 后者为1.36 m/s。在使用数值预报模式或散射计的风向数据提取风速时, 风向数据的空间分辨率较低, 需要将风向插值到SAR图像的各像元上, 忽略了许多小尺度上的风向变化, 从而影响风速的精度。同时, 风向数据的获取时间可能与SAR图像的获取时间不同, 时间差越大, 对风速精度的影响越大^[7-8]。

在使用地球物理模型函数从SAR图像中提取风场信息时, 需要输入风向数据, 而上述风向被用于风场反演时存在各种问题。针对风场反演中存在的问题, 本文提出了组网SAR卫星的一种工作模式, 在该模式下多颗星载SAR同时以不同入射角观测同一海面。目前在轨的组网卫星无法实现多颗SAR同时以不同入射角观测同一海面, 因此, 本实验通过SAR成像仿真的方法获取了组网SAR卫星数据, 并建立了适用于组网SAR卫星数据的风场优化反演方法。该方法可以在不输入风向的条件下提取风速信息, 但是在中高海况下的风速反演精度较低。在一定条件下, 将海浪方向作为参考风向用于筛选结果, 可以获取较准确的风向信息, 并进一步提高风速的精度。

1 数据源

目前, 在轨的组网卫星无法实现多颗SAR同时以不同入射角观测同一海面。针对这种情况, 本文采用SAR成像仿真技术模拟了组网SAR卫星数据, 以此作为后续研究的基础数据。

1.1 SAR成像仿真的基本原理

SAR成像仿真由海浪谱模拟、海面模拟、SAR回波模拟、SAR海面成像等步骤组成, 具体实现流程如图 1。

SAR成像仿真的第一个工作是海浪谱模拟, 海浪谱模拟是海面模拟的前提。常用的海浪谱有PM谱^[9]、JONSWAP谱^[10]、E谱^[11], 本文选用E谱进行海浪谱模拟, E谱描述的是全波数谱。E谱的低波数部分(B_l)如下:

$ {B_{\rm{l}}} = \frac{{{\alpha _{\rm{p}}}}}{2} \cdot \frac{{{c_{\rm{p}}}}}{c} \cdot {F_{\rm{p}}}, $

(1)

其中, α_p是长波的平衡区域参数, c_p是与峰值波数对应的相速度, c代表相速度, F_p是长波边缘效应函数。E谱的高波数部分(B_h)如下:

$ {B_{\rm{h}}} = \frac{{{\alpha _{\rm{m}}}}}{2} \cdot \frac{{{c_{\rm{m}}}}}{c} \cdot {F_{\rm{m}}}, $

(2)

其中, α_m是短波的平衡区域参数, c_m代表最小波相速度, c代表相速度, F_m是短波边缘效应函数。全波数谱(S)如下:

$ S = {k^{ - 3}}\left( {{B_{\rm{l}}} + {B_{\rm{h}}}} \right), $

(3)

其中, k代表海浪波数。实际海浪除了沿主波方向传播外, 还向其他方向扩散, 因此, 第一项工作是需要借助方向函数描述不同传播方向的海浪。将一维海浪谱乘以归一化的方向函数即可得到二维海浪方向谱, 然后根据重力色散关系, 把二维海浪方向谱转换为二维波数方向谱。本文选用的方向函数(f)如下:

$ f\left( {k, \varphi } \right) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\left[ {1 + {\rm{\Delta }}\left( k \right)\cos \left( {2\varphi } \right)} \right], $

(4)

其中, φ代表与风向相对的海浪方向, 即风向与浪向的夹角。波数方向E谱(E)可由全波数谱和方向函数表示, 如下:

$ E\left( {k, \varphi } \right) = S \cdot f\left( {k, \varphi } \right)/k. $

(5)

第二项工作是海面模拟, 海面模拟为海面后向散射系数的计算提供背景场。之后, SAR才能对海面发射电磁波并接收回波, 从而结合电磁散射模型计算海面后向散射系数。本文采用蒙特卡罗方法即线性滤波法模拟二维海面^[12]。

第三项工作是SAR回波模拟。在模拟海面回波数据的过程中, 有两个需要解决的关键问题: 海面后向散射系数的计算、海面回波信号的生成。本文选用双尺度电磁散射模型计算海面后向散射系数, 规定SAR的波段为C波段, 极化方式为VV极化, SAR的具体系统参数如表 1所示。模拟回波信号的算法有时域算法和频域算法两种。频域算法计算量相对较小, 但在生成海面回波信号的过程中无法引入速度聚束效应, 生成的海面回波信号精度受限。时域算法模拟的是SAR接收机的真实工作过程, 能够获取原始回波信号, 虽然计算量大, 但生成的海面回波信号最为准确^[13]。由于对精度的需求, 本文选用时域算法, 假定SAR工作方式为正侧视, 模拟海面回波信号。

最后一个工作是SAR海面成像。本文采用距离多普勒(range Doppler, RD)成像算法对海面回波信号进行成像处理, 得到海面的SAR图像。RD算法包括3个关键步骤: 距离向压缩、距离徙动校正、方位向压缩。

通过上述流程获取的仿真SAR图像符合SAR探测海面的原理, 且可用于实验研究^[14]。本文使用的实验数据如表 2所示。在优化风速反演模型时, 本文使用了仿真的单星SAR数据。为保证优化模型对仿真SAR数据的适用性, 本文模拟了多种条件的SAR数据, 条件如下: 入射角为20°~60°、相对风向为0°~360°、风速为1~30 m/s。在进行多入射角风场反演时, 本文模拟了3颗星载SAR同时以不同入射角观测同一片海面, 观测海面时的风速为1~30 m/s、相对风向为0°~180°, 3颗SAR的入射角分别为: 31°~40°、41°~50°、51°~60°。

2 组网SAR卫星风场优化反演 2.1 理论背景

目前, 地球物理模型函数是C波段SAR数据海面风速反演的常用方法。该模型是通过统计大量的雷达后向散射系数与相应位置的浮标或者数值预报模式资料而建立的经验模型, 它描述了雷达入射角、风速、相对风向和后向散射系数之间的关系。目前常用的地球物理模型函数主要有CMOD4、CMOD-IFR2、CMOD5和CMOD5.N。CMOD4是欧洲中期天气预报中心(European centre for medium-range weather forecast, ECMWF)根据ERS-1卫星上搭载的C波段散射计结合数值预报模式风场数据拟合得到的风场反演模型, 后来被证实同样可用于C波段同极化SAR数据, 是目前最具代表性的地球物理模型函数^[15]; CMOD-IFR2是法国海洋开发研究院根据NOAA浮标数据和ERS系列数据开发的适用于C波段同极化SAR数据的地球物理模型函数^[16]; Hersbach等在CMOD4基础上开发得到更适用于高风速海况的CMOD5, 该模型在台风、飓风等高海况下表现出较好的风速反演效果^[17]; CMOD5.N是对CMOD5的28个可调系数进行重新拟合得到的模型, 一方面修正了CMOD5存在的0.5 m/s的低估偏差, 另一方面因其针对中性风设计, 更能代表海表状况, 避免了大气分层可能带来的误差^[18]。

本文选择CMOD5.N作为研究的基础模型, CMOD5.N的函数形式如下:

$ {\sigma ^0} = {B_0}{\left( {1 + {B_1}\cos \phi + {B_2}\cos 2\phi } \right)^{1.6}} $

(6)

其中, σ⁰为后向散射系数, ϕ为相对风向, B₀、B₁和B₂为风速v和入射角θ的函数。

在SAR风场反演中, 直接使用CMOD需要输入风向数据。因此, 风速反演精度受风向精度的影响。散射计能以不同的方位角对同一海面进行多次有效观测, 因此, 可以利用极大似然估计技术同时提取出风速和风向。利用这一原理, 将ENVISAT ASAR图像分割为多个25 km×25 km的单元, 相邻单元之间的入射角差大于1°。在已知其他相关参数的情况下, 根据相邻单元入射角的差异建立代价函数, 通过一系列求解过程获取风矢量^[19]。该方法为SAR风场反演提供了新的思路。但该方法仅适用于风场分布均匀且变化相对缓慢的海况, 且需要外部风向数据来确定唯一解。将该方法应用于星载SAR数据时, 需要保证相邻子图像的入射角度不同, 相邻子图像的范围相对较大, 风场反演的分辨率较低。

利用组网SAR卫星观测海面, 获得同一海域不同入射角下的SAR数据, 可以保证不同SAR数据之间海面风场的一致性。因此, 将多入射角风场反演方法应用于组网SAR卫星数据时, SAR风场反演的分辨率不受入射角差异的限制, 且可处理的SAR数据不受海况约束。

2.2 星载SAR真实数据的风速反演

为了确定GMF对仿真SAR数据的适用性, 本文先验证了将GMF应用于真实SAR数据时的风速反演精度, 使用CMOD5.N提取了RADARSAT-2数据和Sentinel-1A数据的风速信息, 在提取过程中输入的风向数据为ECMWF风向。最后, 使用ECMWF风速验证反演风速的准确性。在海面风速为3~13 m/s的海况下, 风速反演的RMSE为0.86 m/s, 满足海面风速反演的精度要求。反演风速和ECMWF风速之间的对比如图 2所示。

2.3 星载SAR仿真数据的风速反演

为了验证CMOD5.N对星载SAR仿真数据的适用性, 本文处理了入射角为20°~60°、风速为1~30 m/s的仿真SAR数据。根据不同的入射角将数据分成41组, 使用CMOD5.N提取每组数据的风速信息, 然后将反演风速与实际风速进行对比。部分对比结果如图 3所示。

从图 3中可以发现: 当实际风速大于10 m/s时, 反演风速与实际风速的差值逐渐增大。这是因为仿真的后向散射系数与实际的后向散射系数不完全一致, 当入射角为40°, 相对风向为45°时, 后向散射系数与风速的关系如图 4所示。

CMOD5.N是一个拟合真实SAR数据与海面风场信息的经验函数, 假设CMOD5.N中的后向散射系数与输入参数之间的关系和现实世界是一样的。从图 4中可以看到: 当风速超过10 m/s时, 仿真后向散射系数增长率逐渐小于真实后向散射系数增长率。原因是SAR数据的仿真过程没有考虑泡沫对后向散射系数的影响, 在风速小于10 m/s时, 泡沫引起的海面辐射率变化相对较小, 但随着风速和白浪覆盖面积的增大, 这种变化迅速增大^[20]。

从图 3中还可以看出: 当入射角小于27°时, 反演风速的分布无规律且波动大; 而当入射角不小于27°时, 反演风速的分布有规律且波动小。这种现象也与仿真SAR后向散射系数相关。如图 5和图 6所示, 分别展示了后向散射系数与相对风向、风速的关系, 其中, 虚线代表通过CMOD5.N计算的后向散射系数, 实线代表仿真的SAR后向散射系数。

从图 5中可以看出: 随着入射角的增大, 仿真后向散射系数随着相对风向变化的趋势越来越近似余弦形状, 也逐渐接近CMOD5.N中后向散射系数随着相对风向变化的趋势; 当入射角大于24°时, 仿真后向散射系数随着相对风向变化的趋势基本近似余弦形状, 且不存在明显的异常值。从图 6中可以看出: 随着入射角的增大, 仿真后向散射系数随着风速变化的波动越来越少, 变化趋势逐渐变得平缓; 当入射角大于25°时, 不同风速对应的仿真后向散射系数均小于通过CMOD5.N计算的后向散射系数, 不再存在明显的异常值点。因此, 为了保证数据的可用性, 将合格数据的入射角临界值设置为27°, 将入射角小于27°的仿真SAR数据视为不合格的仿真数据, 不再进行处理。

接下来, 统计了不同入射角条件下的反演风速的均方根误差, 结果如图 7所示。34组数据的RMSE在5.5 m/s至8 m/s之间, 这说明: 直接使用CMOD5.N提取仿真SAR数据的风速, 不能满足海面风速反演的精度要求。为了保证后续研究的可靠性, 有必要找到适用于仿真SAR数据的风场反演模型。

本文使用的数据是组网SAR卫星的模拟观测结果。通过更改输入参数, 可以获得各种海况下的SAR数据。首先, 使用CMOD5.N提取不同海况下的仿真SAR数据的风速; 然后, 拟合反演风速与实际风速之间的关系并建立关系数据库; 最后将拟合关系添加到CMOD5.N, 获得适用于仿真SAR数据的风场反演模型。本实验采用最小二乘法拟合反演风速与实际风速之间的关系, 获得三阶非线性关系的拟合函数。反演风速与拟合关系的分布如图 8所示。

从图 8中可以看出: 拟合关系的形状相似, 且与反演风速的拟合度高。将拟合的非线性函数添加到CMOD5.N获得优化的模型, 通过优化模型提取仿真SAR数据的风速。反演风速与实际风速之间的对比如图 9所示。

从图 9可以看出: 通过优化模型提取的反演风速与实际风速的一致性较好。接下来, 计算不同入射角条件下的反演风速与实际风速相比的均方根误差, 结果如图 10所示。

从图 10中可以看出: 在不同入射角的条件下, 通过优化模型提取风速的RMSE均小于1.5 m /s。这说明采用优化模型提取的风速精度符合海面风速反演的精度要求, 进一步证实, 优化模型适用于多数入射角和海况下的仿真SAR数据。

2.4 多入射角风场反演

利用2.3节优化的风速反演模型建立代价函数。在建立代价函数时, 使用不同入射角的SAR数据。本研究使用了3颗卫星的SAR数据, 建立的代价函数如公式(7)所示。

$ \begin{array}{l} {J_{{\rm{cost}}}}\left( {\theta , \Phi , v} \right) = {\left[ {\sigma _1^{\rm{m}}\left( {{\theta _1}, {\Phi _1}, v} \right) - \sigma _1^0\left( {{\theta _1}, {\Phi _1}, v} \right)} \right]^2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {\left[ {\sigma _2^{\rm{m}}\left( {{\theta _2}, {\Phi _2}, v} \right) - \sigma _2^0\left( {{\theta _2}, {\Phi _2}, v} \right)} \right]^2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {\left[ {\sigma _3^{\rm{m}}\left( {{\theta _3}, {\Phi _3}, v} \right) - \sigma _3^0\left( {{\theta _3}, {\Phi _3}, v} \right)} \right]^2} \end{array} $

(7)

其中, $ \sigma _1^{\rm{m}}, \sigma _2^{\rm{m}}和 \sigma _3^{\rm{m}} $是通过优化模型计算的后向散射系数, $ \sigma _1^{\rm{0}}, \sigma _2^{\rm{0}}和 \sigma _3^{\rm{0}} $是仿真的组网SAR后向散射系数, θ为入射角, Φ为相对风向, v为风速。当代价函数的值最小即接近0时, 则通过模型计算的后向散射系数与仿真的SAR后向散射系数最接近, 此时的风速与风向就是反演结果。

为了最小化代价函数, 分别获取了代价函数关于风速和相对风向的偏导数, 如式(8)、(9)所示, 式中字母含义同式(7)。由于优化模型是一个相对复杂的函数, 因此式(7)、(8)、(9)也是复杂的函数。为了减轻计算难度, 本文使用了一种双精度搜索方法来获取目标代价函数的最小值, 从而获得海面风速和风向的解。搜索过程主要分为两个步骤: 粗搜索和细搜索。在粗搜索过程中, 将0~30 m/s的风速以1 m/s为间隔, 将0°~360°的相对风向以10°为间隔代入式(8)与式(9), 得到偏导数。在结果中, 一些相邻结果的符号不同, 说明在相邻数据之间存在一个偏导数为0的点, 即极值点。然后, 在存在极值点的区间内进行细搜索, 将风速以0.1 m/s为间隔, 将相对风向以1°为间隔代入式(8)与式(9)。将偏导数结果中最接近0的值保存, 并将其对应的风速和相对风向代入式(7)以计算代价函数, 与代价函数最小值对应的风速是风速反演结果, 由于风向存在180°模糊的问题, 故无法获取准确的风向结果。

$ \frac{{\partial {J_{{\rm{cost}}}}\left( {\theta , \Phi , v} \right)}}{{\partial v}} = \mathop \sum \limits_{{\rm{i}} = 1}^3 2\left( {\sigma _i^{\rm{m}} - \sigma _{\rm{i}}^0} \right)\frac{{\partial \sigma _{\rm{i}}^{\rm{m}}}}{{\partial v}}, $

(8)

$ \frac{{\partial {J_{{\rm{cost}}}}\left( {\theta , \Phi , v} \right)}}{{\partial {\rm{cos}}\left( \Phi \right)}} = \mathop \sum \limits_{{\rm{i}} = 1}^3 2\left( {\sigma _{\rm{i}}^{\rm{m}} - \sigma _{\rm{i}}^0} \right)\frac{{\partial \sigma _{\rm{i}}^{\rm{m}}}}{{\partial {\rm{cos}}\left( \Phi \right)}}. $

(9)

该方法可以解决SAR风速反演依赖于风向数据的问题, 但只能获得反演风速。当相对风向为45°时, 对实际风速为1~30 m/s的30组数据进行风速反演, 对比反演风速与实际风速, RMSE为2.25 m/s, 相关系数为0.975。反演风速与实际风速的对比如图 11所示。

另一种筛选结果的方法需要将参考风向引入处理流程。在计算完代价函数后, 设置合适的筛选阈值, 将筛选后的结果视为风矢量的模糊解。最后, 使用参考风向消除模糊解, 获得唯一解。前文提到了获取风向的方法, 散射计数据和数值预报模式数据的分辨率低; 基于风条纹信息获取的反演风向存在180°模糊, 且大约40%的SAR图像没有风条纹。因此, 本文使用海浪方向来消除模糊解, 从而获得唯一解。从SAR图像中可以获得不存在180°模糊的海浪方向, 将海浪方向作为参考风向可以增强SAR风场反演的独立性。

海浪包括风浪和涌浪。风浪是在区域风的直接作用下形成的, 从开始持续成长。当风速突然下降或风向突然改变时, 风浪可能会成为涌浪。当海浪是风浪时, 海浪的方向可以在某种程度上代表风场的方向。为了验证使用风浪方向代替风向进行风场反演的可行性, 本文使用ECMWF平均风浪方向作为CMOD5.N的输入风向来反演真实SAR数据的风速。实验数据是2.2节中使用的RADARSAT-2和Sentinel-1A数据。反演风速和ECMWF风速的对比如图 12所示。

如图 12所示, 当输入风向为风浪方向时, 反演风速的RMSE为0.91 m/s。这说明: 在一定条件下, 将风浪方向代替风向用于地球物理模型函数是可行的。接下来, 对比用于实验的风浪方向与风向, 风浪方向与风向相比的RMSE为3.87°, 这表明: 风浪方向与风向具有较好的一致性。对比结果如图 13所示。

根据海浪谱的能量密度可以获得海浪方向, 海浪谱的能量在海浪方向上最强。已有的海浪谱反演方法包括: MPI算法、交叉谱算法、PARSA算法和SPRA算法。MPI算法是经典的海浪谱反演方法, 但无法解决海浪方向存在180°模糊的问题^[21]。交叉谱算法解决了海浪传播方向的180°模糊问题, 但反演结果受截断波长的影响^[22]。PARSA算法是MPI和交叉谱算法的组合, 可以获取完整的二维海浪谱, 且解决了海浪方向的180°模糊问题^[23]。SPRA算法只获得涌浪谱, 并依赖于散射计信息^[24]。

本文将反演的海浪方向作为参考风向引入多入射角风场反演方法, 从仿真SAR数据中提取的海浪方向与实际风向的大小一致。对风向为45°、风速为1~30 m/s的仿真SAR数据进行处理, 在得到多个模糊解后, 将反演的海浪方向作为参考风向用于确定唯一解, 得到风速和风向的反演结果。反演风速、风向与实际风速、风向的对比结果如图 14、15所示。

从图中可以看出: 反演风速的RMSE为0.61 m/s, 相关系数为0.998;反演风向的RMSE为3.21°, 相关系数为1。这说明: 引入准确的参考风向可以提高风速反演的精度, 并使该方法可以提取精度较高的风向信息。

3 结论

针对单星SAR时间分辨率低且风场反演依赖风向的问题, 本文提出了一种卫星组网的工作模式, 针对该工作模式进行了SAR数据仿真与风场反演研究。实验与验证结果说明: 在不引入参考风向的前提下, 仅通过最小化代价函数筛选结果, 则风速反演精度在中低风速时较高, 但无法提取准确的风向信息; 在引入准确的参考风向后, 各海况下的风速、风向反演精度都较高。本文的研究说明, 通过卫星组网获取的不同入射角的SAR数据可以为海面风场反演提供更多的可能, 可以使SAR风速反演不依赖风向数据, 使基于入射角差异的风场反演不受海况与分辨率限制。当海浪以风浪成分为主时, 海浪方向可以被作为参考风向来筛选结果, 提高SAR风场反演的精度。本文的研究可以在一定程度上为未来C波段组网SAR卫星的工作模式及风场反演的研究提供一些参考。