基于ENVISAT ASAR 数据的龙目海峡海表流场反演

李晴烁，范陈清，孟俊敏，张杰，柯长青

（1.南京大学 地理与海洋科学学院，江苏 南京 210023；2.自然资源部第一海洋研究所，山东 青岛 266061）

海流是指海水较为稳定的大尺度流动，代表在准定常状态下海水不同规模的运动状态，其空间运动规模可达几百千米甚至全球范围，时间跨度可为月尺度、季节尺度甚至多年尺度[1-2]。海流决定了海洋中的热能和物质（如营养盐、污染物、泥沙等）在不同海区之间的输送和分布，包含了海洋中各种物化过程，对海洋的气候变化、海洋生物的生存和地形变迁都会造成直接影响。此外，海洋溢油、海上运输线路的优化、海上搜救和海洋环境监测等均需考虑流场的影响。因此，观测海洋流场的变化并监测分析其特性与规律，对海洋科学的研究具有重大意义[1,3]。

海表流场的观测手段主要包括传统现场观测和卫星遥感观测，传统现场观测主要包括浮标观测等，其观测精度较高，但现场测量的仪器安放危险性高、观测范围有限、观测成本高，具有较大的限制性[2]。卫星遥感可以对海流进行大范围的观测，利用该手段不仅可以获得完整的二维流场结构，还可以分析海面流场随时间变化的情况。卫星遥感主要有3 种观测手段：①雷达高度计[4]。利用雷达高度计测得海面高度参数，结合地转平衡关系获取流场的流速与流向，该方法具有观测面积大、精度高和时间准同步等优点，但仅适用于测量大尺度的地转流变化，且存在重复观测周期较长的问题[4-6]。②高频地波雷达。基于多普勒效应和布拉格散射理论对海表流场进行监测，该方法可连续获取大面积海流且不受天气状况影响，但目前业务化应用较少，且剖面流场随深度的不均匀性和流场剖面的时变情况会引起径向流速的测量误差[7]。③合成孔径雷达（Synthetic Aperture Radar，SAR）。它具备高分辨率、不受天气状况影响等优势[8]，其反演流场的方法主要有2 种：一种是顺轨干涉测量方法（Along Track Interferometry,ATI），该方法利用影像的相位差获取海表运动信息，空间分辨率高，但成像条件受限且数据成本较高，适用范围小；另一种是本文所采用的多普勒频移技术（Doppler Centroid Anomoly，DCA），该方法利用SAR 与地球相对运动、地球自转、海表面运动产生的回波信号频率偏移获取海表流场，分辨率较高、应用范围较广，适用于大范围的海流探测[9-10]。

目前已有许多国内外学者将多普勒频移技术应用于海表流场反演。国内方面：李慧敏[2]利用ENVISAT ASAR IMS（Image Mode SLC）数据，采用符号多普勒法反演了长江口近岸海域的海表径向流速，反演结果与潮汐表的结果较为一致。杨小波[11]基于ENVISAT ASAR WSM（ Wide Swath Mode）数据，提出了方位向和距离向多普勒中心模型，采用多普勒频移法反演长江口及杭州湾近海港口海表面流速，反演结果与模式结果的差值为0.1 m/s。候富城[12]基于ENVISAT ASAR IMS 数据，采用多普勒频移法反演海表流速并提取内波，反演得到的内波致海表面流的方向与内波传播方向在雷达距离向上的分量一致。宋小霞等[13]利用ENVISAT ASAR WSM 数据，采用多普勒频移法反演厄加勒斯和黑海的海表流场，流场结果与AVISO 流场的相关系数为0.84。国外方面：Johannessen等[14]基于ENVISAT ASAR WSM 数据，分析了从多普勒频移中反演得到的海流速度组成成分，并由此反演了厄加勒斯海表流速，反演结果与高度计所测结果较为一致。Rouault等[15]利用ENVISAT ASAR WSM 数据，采用多普勒频移法计算了厄加勒斯海表平均径向流速，并以陆地覆盖影像为参照校正频移误差，计算结果与AVISO 数据匹配度较高。Hasen等[16]利用ENVISAT ASAR WSM数据，基于多普勒频移法对挪威海沿岸海域进行了流场反演研究，并校正了多普勒中心计算中存在的卫星姿态运动导致的误差，反演得到的流场结果在10 m 分辨率下误差小于5 cm/s。这些研究结果表明，SAR 多普勒技术可有效用于海表流场探测。

龙目海峡是印尼贯穿流的重要出口通道之一，系统分析其海流流速和季节变化有助于分析印尼贯穿流在气候和海洋中的作用，前人研究多基于数值模拟分析班达海和弗拉瑞斯海[17-19]，针对龙目海峡流速和流向的研究较少。近年来，SAR 多普勒技术的出现为系统分析龙目海峡海流变化提供了新的手段[1-3,20]。本文基于2003—2004 年ENVISAT ASAR WSM 数据和ENVISAT ASAR IMS 数据，利用多普勒技术反演了龙目海峡海表流场，并分析了龙目海峡海表流场的变化规律，以期为研究龙目海峡流场的动态变化提供有力的技术和数据支撑。

1 研究区域与数据来源

1.1 研究区域

龙目海峡（114°24′～116°42′E，6°30′～10°00′S）位于印度尼西亚群岛的巴厘岛和龙目岛之间，是地处印度洋和太平洋之间的海上交通要冲，其南北全长80.5 km，水深呈北深南浅分布，最大水深可达1 306 m，受强烈海流的侵蚀，近些年来水深仍在增加。龙目海峡是世界级的海运通道，同时也是印尼贯穿流的重要出口通道之一，局地潮致与海气相互作用混合使得该海域混合效应增强，其不仅对海上输运和交通意义重大，而且对军事活动也具有重要的影响[19]。

研究区域为龙目海峡部分海域，覆盖范围为（114°48′～116°00′E，8°30′～9°48′S）（图1 红色框区域），为了便于分析对比多普勒频率的差异，本文选取影像包含陆地区域。

图1 研究区域Fig.1 Study area

1.2 数据来源

1.2.1 ASAR 数据

本文采用欧洲航天局（European Space Agency,ESA）ENVISAT 卫星搭载的ASAR WSM[21]和IMS[22]产品作为数据源。ENVISAT 卫星运行时间为2002—2012 年，轨道高度为800 km，轨道倾角为98°，工作波段为C 波段，重访周期为35 d[21-23]。其中，WSM 数据影像幅宽约为420 km，最小时间间隔为3 d，空间分辨率为150 m，有利于观测大范围海表特征[13,21-23]。IMS数据含有相位信息和斜距信息，影像幅宽为100 km、空间分辨率为4 m，有利于小尺度海洋现象观测[2,12-13,23]。选 取2003 年5 月 至10 月 的ENVISAT ASAR WSM 产品和2003 年11 月至2004 年6 月的 ENVISAT ASAR IMS 产品对龙目 海峡研究区域海流进行反演，数据有关参数如表1 所示。

表1 ENVISAT ASAR WSM 和IMS 图像参数Table 1 ENVISAT ASAR WSM and IMS image parameters

1.2.2 AVISO 数据

AVISO 数据因其覆盖范围广、数据质量较高和易于获取等优势，在验证ASAR 流场精度方面发挥了重要的作用[14,16]，本文选取AVISO 的GLORYS12V1 版本流场产品[24]作为验证数据。GLORYS12V1 在欧洲CMEMS（Copernicus Marine Service）框架下构建，由ERA-Interim（ECMWF Re-Analysis Interim）和 ERA5 再分析数据融合所得，主要包括流场、温度、盐度、海平面、混合层深度和冰参数的日平均和月平均数据。时间间隔为1 d，空间最小分辨率为8.3 km，产品预测的时间长度为1993 年1 月1 日至2019 年12 月31 日。GLORYS12V1 的流场数据与潜标观测数据较为一致，大部分流场数据与潜标数据的均方根误差小于0.25 m/s；在西太平洋吻合度最高，与潜标数据的均方根误差小于0.20 m/s；在赤道潜流处吻合度最低，与潜标数据的均方根误差小于0.30 m/s；随着深度的增加，GLORYS12V1 的流场数据与潜标数据的相关性逐渐减小[25]。

2 方法

2.1 多普勒质心异常计算

多普勒质心的计算是合成孔径雷达（SAR）成像的重要部分。WSM 数据的多普勒质心可由其多普勒网格数据插值得到（图2），由WSM 数据插值得到的多普勒中心频率fDc在方位向分辨率约为8 km，距离向分辨率最高约为9 km，最低约为 3.5 km[11,26]。为方便计算，后文WSM 数据的多普勒质心和多普勒质心异常均以类似图2a 的时间多普勒网格坐标成像。

多普勒中心频率fDc包含由海表运动产生的多普勒质心异常fDca和由卫星与地球相对运动产生的卫星多普勒频移fSat，为获得fDca需去除fSat影响[27]。根据星载SAR 多普勒频率理论模型，WSM 产品的成像时刻基带中心多普勒频率fDc分别在距离向和方位向具有近似线性变化的特点[11,13,28-30]。沿距离向选取截面A 和B，沿方位向选取截面C 和D，其中截面A 和D 含有陆地区域，截面B 和C 不包含陆地区域（图2a）。由图2b至图2e 可见，fDc在距离向和方位向均呈线性变化，截面A、B 与距离向时间的拟合优度分别为0.94、0.95，截面C、D 与方位向时间的拟合优度分别为 0.96、 0.91。拟合优度R2是回归直线对观测值的拟合程度，R2值越接近1 说明回归直线对观测值拟合程度越好，也说明多普勒质心在距离向和方位向上呈线性变化，因此分别在2 个方向上采用均值滤波和线性拟合的方法去除卫星多普勒频移fSat影响是可取的，且该方法也适用于含有陆地的WSM 影像[13]。

图2 2003 年6 月29 日WSM 数据的多普勒中心网格及网格截面在距离向、方位向的变化趋势Fig.2 Doppler centroid grid and the trend of the grid section in range and azimuth directions of WSM data on June 29,2003

频谱能量均衡法通过分析回波的多普勒谱计算fDc。由于多普勒谱与能量谱的特性相似，根据能量均衡理论可将频率谱划分为两部分，若这两部分的能量均衡对称[31]，则中心均衡点的频率即为多普勒质心频率。基于以上理论，采用频谱能量均衡法计算IMS 产品的fDc，具体处理过程为：首先，对影像的局部窗口进行FFT 变换，同时考虑多普勒中心频率的精确度与其分辨率大小，将窗口大小设置为512×128。然后，将变换结果与 [−fs/2，fs/2] 对应，采用正弦曲线拟合方位向上的多普勒谱能量，则该窗口的多普勒中心频率即为拟合得到的曲线中心点频率。整景影像的多普勒中心频率fDc可通过滑动窗口得到，能分辨的最小单元为100 m[12]。采用Tie-point 网格数据中的斜距时间和多普勒系数计算IMS 产品的卫星频移[12]，计算公式为：

式中：fcoef为多普勒系数；t为斜距时间；t0为标准斜距时间。

在多普勒质心频移fDc中减去卫星多普勒频移fSat影响后可计算得到fDca，即：

选取2003 年6 月29 日 WSM（图3）数据和2004 年4 月7 日 IMS（图4）数据，分别针对这2 组数据计算多普勒质心，结果显示：在WSM 多普勒质心频率影像中（图3b），方位向多普勒网格13 处频移误差较大，这可能是由海洋和陆地的电磁波散射特性差异较大所致[13]。fDc图像在朝向卫星方向上海洋和陆地相交处含有明显的透视收缩，而在背离卫星方向上只有在海洋和陆地相交的区域才有阴影，由于图像透视收缩及阴影的范围仅出现在海陆交界处，其对反演海流的影响可以忽略[11，13]，本文在海洋和陆地相交处剔除掉由透视收缩和阴影产生的影像异常值，最终得到2003 年6 月29 日 WSM 整景影像的多普勒频率（−50～150 Hz，图3d）。2004 年4 月7 日 IMS 整景影像的多普勒频率为−30～15 Hz，在IMS 的多普勒质心异常图像（图4b）上可观察到明显的海浪波纹。对比WSM 数据fDca（图3d）和IMS 数据fDca（图4b）发现，相较WSM 数据，IMS数据的多普勒质心异常图像分辨率更高，更适合分析小范围的海表面运动。

图3 2003 年6 月29 日 WSM 影像数据、fDc、剔除异常值后的fDc和fDca 结果Fig.3 WSM image data,results of fDc,fDc after eliminating outliers and fDca on June 29,2003

图4 2004 年4 月7 日 IMS 数据fDc与fDca 结果Fig.4 Results of Doppler Centroid fDc and Doppler Centroid Anomaly fDca of IMS data on April 7,2004

2.2 海表面信息反演

多普勒质心异常fDca中包含了布拉格波[11-13]、风场[16,32-34]和海表面流速等的贡献，所以fDca的计算公式为：

式中：fB为布拉格散射导致的多普勒频移；fw为风场影响导致的多普勒频移；fs为海表流场引起的多普勒频移。根据式(3)去除布拉格波和风场的多普勒频移影响后，即可获得海表流场的多普勒信息。

通常将海面与雷达入射电磁波产生共振的波浪定义为布拉格波，其导致的多普勒频移fB通常表示为[33]：

式中：g为重力加速度；Kr为雷达入射电磁波波数，对于ASAR 数据Kr取112 m/s；θi为入射角。布拉格波导致的多普勒频移的正负是由海流的传播方向与雷达视向方向相同或相反决定的[32]。

多普勒质心异常中去除布拉格波产生的频移后，即可计算出雷达视向的海表运动速度vslc，即：

为了去除海面风场引起的多普勒频移，本文利用ASAR 数据反演海面风场信息。通常认为，海面上方10 m高度风速与雷达径向风速vw间有如下关系[27]：

式中：γ是与风速大小有关的系数，在真实海洋中，γ通常取0.03[11-12]；U10为雷达径向方向上海面10 m 高度的风速。

本文采用CMOD 5（C-band models 5）模型反演海面风场，并结合ECMWF[35]风场的风向信息解决反演模型的180°模糊问题[12,34]。在实际SAR 图像中，分析窗口的大小对风向反演结果具有重要意义：窗口大小为20 km的反演结果可以体现出风向局部细节，有利于小尺度风场反演；窗口大小为40 km 的风向反演结果具有相对综合性，在大尺度上可体现出风向的一致性，在小尺度上也能较好地体现出风场的局部细节，有利于中尺度风场反演；窗口大小为60 km 的反演结果可充分展示大范围风场结构，有利于大尺度风场反演[11-12]。

对于风速数据的反演，本文将VV 极化的归一化雷达散射截面（Normalized Radar Cross Section,NRCS）数据应用到CMOD 5 模型中，但此模型不适用于HH 极化数据。所以，需要将 HH 极化的 NRCS 转换为相应的VV 极化的NRCS[34]，即：

多普勒质心异常中去除风场引起的多普勒频移后，即可计算出在雷达视向海表流速的真实值vs（风速的正负由风向与雷达视向方向相同或相反决定）[32]：

式中：vslc为雷达视向的海表运动速度；vw为雷达径向风速。

为进一步表示雷达径向的海表流场速度，根据雷达视线方向海流运动速度与海表运动产生的多普勒频移之间的线性关系，可计算出垂直于卫星飞行方向的海表径向流速vd[12-13]：

综合海表面运动信息，利用SAR 多普勒方法可以反演得到沿雷达视向的海表径向流速。

3 结果与讨论

将WSM 宽幅数据裁剪至研究区，利用式（3）～式（9）计算出径向海表流速，4 个季节各选取1 景影像反演径向速度（图5），4 景影像均为降轨数据，朝向雷达径向的流速用正值表示，背离雷达径向的流速用负值表示[13,34]。为了更好地验证流向精度，将AVISO 流场数据转为径向流场（图5 中箭头）。结果表明，反演得到的径向流场与AVISO 径向流场较为一致，其中在近岸海域，IMS 反演的流向与AVISO 流向较WSM 反演的流向匹配度更高，这是因为IMS 影像分辨率更高、幅宽更小，有利于监测和分析小范围的海面流场特征，而WSM 为宽幅影像，更有利于观测较大范围的海面流场特征。

图5 ASAR 数据反演流速（填色区）与AVISO 流场（箭头）分布Fig.5 The retrieved sea surface velocities field from ASAR (color shading) and AVISO (arrows)

为进一步验证反演的流速精度，将反演得到的流速网格分辨率重采样至8.3 km，使其与AVISO 流场分辨率一致，将反演的整幅影像流场与径向AVISO 流场求差值，对流速和流向的误差分布进行分析（图6）。其中WSM 数据反演流速主要分布于−0.5～1.0 m/s，与AVISO 流速误差主要集中于−0.35～0.45 m/s，均方根误差为0.28 m/s，相关系数为0.95。IMS 数据反演流速主要分布于−0.2～0.4 m/s，与AVISO 流速误差主要分布于−0.15～0.55 m/s，均方根误差为0.37 m/s，相关系数为0.73。由图6a 和图6c 可知，与AVISO 流场相比，ASAR 反演流速整体更大，这也符合前人研究所述[2,13]。可见，本文采用的流场反演方法可有效提取海表面信息，且反演精度较高，适用于近岸海域。

图6 2003 年6 月29 日WSM 数据和2004 年4 月7 日IMS 数据反演流速与AVISO 流速对比以及二者误差对比Fig.6 Comparison of the sea surface velocity retrieved from WSM (June 29,2003) and IMS (April 7,2004) with AVISO sea surface velocity and their errors contraction and the error comparison

对WSM 和IMS 数据进行时间序列变化分析，结果如表2 所示。总体来看，龙目海峡流速较大，最大值为−2.56 m/s（负号表示海流的传播方向与雷达视向相反），流速变化明显。结合图5 分析可知：龙目海峡在2003 年6 月29 日流速较大，主要集中在（115°24′～115°48′E，8°30′～9°12′S）附近的海域，其他海域流速较小（−0.50～0.75 m/s），主要流向为西南向；在2003 年10 月15 日流速较小（−0.79～0.92 m/s），主要流向为东南向；在2003 年12 月24 日整体流速较大，出现北向逆流，主要流向转为东北向，几乎与原来相反；在2004 年4 月7 日海流转回南向流。总体来说，研究区范围内，2003 年5 月至9 月流速较大，10 月流速减小，11 月至12 月流速逐月增加，在2004 年3 月流速达到最低值（−0.58 m/s），而后在2004 年4 月至5 月流速逐月增大；在南半球夏季（12 月）和秋季（5 月）出现北向逆流，其余月份皆为南向流。

表2 海表径向流速反演结果Table 2 The retrieved results of radial sea surface velocity

4 结论

基于2003 年至2004 年ENVISAT ASAR WSM 数据和ENVISAT ASAR IMS 数据，利用多普勒技术反演了龙目海峡海表流场，并分析了龙目海峡海表流场的变化，得出以下主要结果。

1）基于2003 年5 月至10 月WSM 数据和2003 年11 月至2004 年5 月IMS 数据，利用多普勒质心偏移法获取了海面运动引起的多普勒质心异常，通过去除风和Bragg 波引起的多普勒频移，最终获得龙目海峡海表径向流速。结果表明：本文的流场反演方法可准确提取海面流速信息，海表面流速反演结果与AVISO 流场匹配度较好，其中WSM 数据反演的流速与AVISO 流场均方根误差为0.28 m/s，IMS 数据反演的流速与AVISO 流速均方根误差为0.37 m/s。

2）基于反演的径向海表流速的反演结果分析龙目海峡2003 年至2004 年流速和流向变化，结果表明龙目海峡流速较大，最大值为−2.56 m/s（负号表示流向与雷达视向相反），流场变化明显，2003 年5 月至9 月流速较大（−2.56～2.51 m/s），10 月流速减小到0.72 m/s，11 月至12 月流速逐月增加（−1.64～1.25 m/s），在2004 年3 月流速达到最低值（−0.58 m/s），而后在2004 年4 月至5 月流速逐月增大（−1.46～1.72 m/s）；在南半球夏季（12 月）和秋季（5 月）出现北向逆流，部分海流流向巴厘海，其余皆为南向流，海流流向印度洋。

3）本文对WSM 和IMS 数据分别采用了不同的方法估算多普勒质心异常：WSM 数据的多普勒质心由网格化多普勒中心估计所得，卫星频移影响采用均值滤波和线性拟合的方法去除；IMS 数据的多普勒质心由频谱拟合法结合能量均衡理论计算所得，卫星频移影响采用斜距时间和多普勒系数去除。结果表明：IMS能反演出海浪波纹信息，在近岸海域IMS 数据反演得到的流向结果与AVISO 流场的匹配度更高，而从整体上来看，WSM 数据反演得到的流速结果与AVISO 流场的匹配度更高，WSM 数据反演流速与AVISO 流速均方根误差为0.28 m/s，相关系数为0.95，IMS 数据反演流速与AVISO 均方根误差为0.37 m/s，相关系数为0.73。综上，IMS 影像幅宽较小、分辨率较高，更适合分析和监测小范围的海洋区域，而WSM 数据为宽刈幅数据，更适合观测较大范围的海表面流场特征。