王佳明 许 可 蒋茂飞
①(中国科学院国家空间科学中心中科院微波遥感技术重点实验室 北京 100190)
②(中国科学院大学 北京 100049)
1998年,Raney[1]提出一种参考合成孔径雷达(Synthetic Aperture Radar, SAR)沿方位向进行合成孔径处理以提高测高精度和方位向空间分辨率的高度计,用于对海冰和近岸地区的观测,称为延迟多普勒高度计,又称合成孔径高度计。合成孔径高度计从理论提出到搭载卫星在轨运行,已走过20余年的发展历程,其应用场景也不再局限于最初设计针对的海冰、近岸地区,而是实现对全球海洋的观测。在这个发展过程中,CryoSat-2卫星的观测数据对合成孔径高度计数据处理的发展做出了重要贡献,Sentinel-3A/B卫星实现了合成孔径高度计对全球海洋的覆盖观测,2020年发射的Sentinel-6/Jason-CS任务将对合成孔径高度计的新技术和应用开展新的探索[2]。
在高度计的数据处理中,通过地面重跟踪处理求解高度计回波模型的参数,进一步从回波波形中提取精确的回波时延、有效波高、后向散射系数等参数值,从而实现精确测高。文献[3]按照Cryosat-2高度计的系统参数推导了合成孔径高度计在有无热噪声的两种情况下海面回波重跟踪距离测量的克拉美-罗下界(Cramer-Rao Bound, CRB),文献[4]进一步比较合成孔径高度计和传统高度计重跟踪得到的距离、有效波高、后向散射系数参数的CRB,证明合成孔径高度计的理论测量精度优于传统高度计的同时,还通过仿真比较了最大似然估计、最小二乘估计和加权最小二乘估计3种重跟踪估计器的均方根误差和克拉美-罗下界的平方根 (root of CRB)随有效波高的变化关系,证明最大似然估计和加权最小二乘估计的重跟踪估计器性能优于最小二乘重跟踪估计器,但是由于使用的半解析模型[4,5]需要进行卷积运算,模型计算效率低,基于半解析模型的加权最小二乘重跟踪估计器不利于进行业务化应用。
本文使用由欧洲航天局(European Space Agency,ESA)发起的合成孔径高度计研究与应用(SAR Altimetry MOde Studies and Applications, SAMOSA)项目公布的合成孔径高度计回波模型[6],简称SAMOSA模型,设计一种加权最小二乘重跟踪估计器,并结合多视回波的统计特性提出了一种新的简化的合成孔径高度计的加权方法,与非加权形式相比,采用新的加权方法可以提高参数估计精度,相比于半解析模型,重跟踪参数求解速度也得到极大提升,使其可以应用业务化数据处理。本文的另一项工作是利用该估计器和Sentinel-3A公开的合成孔径高度计数据进行验证和对比,实验结果证明该估计器在开放海域的海况下能够提升重跟踪的参数估计精度,且重跟踪的距离和有效波高精度在大部分海况下要优于Sentinel-3A的L2级数据产品的精度。
本文的结构如下:第2节介绍加权最小二乘重跟踪估计器的实现方法和原理,第3节利用Sentinel-3A卫星的L1A级数据对重跟踪估计器的性能进行验证,并与Sentinel-3A的L2级数据产品进行比较,证明该加权最小二乘估计器能够提高高度计参数估计精度,第4节对全文进行总结。
重跟踪本质是调整回波模型中的参数使模型能够与回波波形达到最佳匹配,故其核心要素是回波模型和匹配算法。文献[7]总结了其他基于Brown模型发展而来的合成孔径高度计回波模型,但由于仍然保有卷积形式,在实际的数据处理中重跟踪处理的计算时间往往远大于卫星数据的覆盖时间,模型的应用范围有限;在参数匹配算法上,常用的方法有最大似然估计和最小二乘估计,但由于合成孔径高度计的回波模型解析式复杂,难以求解最大似然函数,在目前的数据处理方法中,仍以最小二乘法为主[8]。综上所述,本文采用SAMOSA模型,设计了基于加权最小二乘方法的重跟踪估计器。
由于合成孔径高度计经延迟校正处理后得到的波形是经多视处理后,各个多普勒单元方位向的回波之和,故最终的回波模型表示为式(3),其中Nb表示多视视数
在不同的有效波高下,式(2)、式(3)计算得到的SAMOSA模型如图2所示,与传统高度计的回波模型相似,波形的前沿斜率随有效波高的增加而减缓,峰值随有效波高的增大而减小。
图1 SAMOSA模型示意图
图2 不同有效波高的SAMOSA模型
本文使用Sentinel-3A卫星Cycle 034-Cycle 036的L1A级数据对重跟踪估计器性能进行验证,数据来自于欧洲气象卫星应用组织(european organisation for the exploitation of meteorological satellites)的数据发布网站(https://eoportal.eumetsat.int)。Sentinel-3A卫星运行在轨道高度815 km,轨道倾角98.62°的太阳同步轨道,地面轨道重复周期为27 d,选取的验证数据的3个Cycle覆盖的时间长度约81 d,覆盖的纬度范围为81.35°S~81.35°N[12-14]。
在对比重跟踪参数前,还要对数据进行一定筛选,由于本文主要验证该重跟踪估计器在开放海域的性能,为避免海冰对高度计测量的影响,选择数据的纬度范围在南北纬50°之间,同时为避免近岸区域和浅水区域回波的影响,选择观测海洋区域的水深小于-1000 m。
以Cycle 035的Pass 171的降轨数据为例,对经数据筛选后得到的覆盖南北纬50°之间时长约28.4 min的海洋回波观测数据进行重跟踪,在CPU为i5-8500的同一计算机上的相同软件条件平台下,利用文献[4,5]的3项卷积模型的半解析模型进行重跟踪,模型的详细参数参考文献[5],重跟踪计算时间约为61.9 min;而本文使用SAMOSA模型设计的重跟踪估计器计算时间约为2.8 min,计算时间缩短为半解析模型重跟踪的4.5%。综合考虑计算效率等因素,本文仅对SAMOSA模型的重跟踪估计器性能进行评估。
数据处理流程如图3所示。由于L1A级产品数据提供的是高度计的正交解调信号的I, Q两通道波形信号,需要先对回波信号进行延迟校正和多视处理,实现进行波形重构,波形重构方法参考文献[9,15]。在完成波形重构后,平均每秒可以生成约20个重构波形,称为20 Hz数据,再利用式(17)构造的加权最小二乘重跟踪估计器从重构波形中提取回波时延对应的距离门、有效波高和回波功率。
图3 数据处理流程图
设定回波模型估计参数的初值,并利用式(17)进行迭代直到拟合的模型与回波波形间的误差逐渐收敛到稳定值后,认为此时模型对应的距离、有效波高和后向散射系数为重跟踪得到的20 Hz的高度计测量数据[7],即平均每秒钟包含20个距离单元、有效波高和回波功率数据。其回波波形和SAMOSA模型的拟合结果如图4所示,其中蓝色图线为L1B级波形数据,红色图线为重跟踪的拟合波形。
图4 重跟踪拟合波形
分别统计式(20)-式(22)得到的标准差随有效波高在0~10 m范围内的散点图分布情况,为便于表示将其转化为散点密度图,最终标准差随有效波高的分布情况如图5-图7的(a), (b)所示,并以0.5 m有效波高为间隔统计有效波高范围内的标准差的中位数作为参数估计精度。最后,比较重跟踪得到的卫星到海面距离、有效波高和归一化后向散射系数的精度,并将重跟踪得到的数据与对应的L2级产品数据进行比较,以评估重跟踪估计器的性能。
图5为1 Hz的卫星到海面距离精度在非加权和加权两种情况下随有效波高变化的分布情况,图5(c)将两种重跟踪方法的精度与Sentinel-3A的L2级数据计算得到的距离精度进行对比,其中蓝色图线为L2级数据产品距离精度,红色图线为最小二乘重跟踪距离精度,黄色方形图线为加权最小二乘距离精度,从中可以看出,最小二乘得到的距离精度与Sentinle-3A的L2级距离精度相近,当Hs=2 m时,最小二乘的距离精度为1.1 cm,加权最小二乘的距离精度为1.0 cm,测距精度提高了9%;当Hs=4m时,最小二乘的距离精度为1.5 cm,加权最小二乘的距离精度为1.3 cm,测距精度提高了13%,距离精度提升统计结果如表1所示。随着有效波高的增加,加权最小二乘重跟踪得到的距离精度整体上优于Sentinle-3A的数据产品精度。
表1 加权最小二乘重跟踪距离精度提高统计
图5 距离精度随有效波高的变化
1 Hz有效波高精度统计结果如图6所示,图6(c)是两种方法重跟踪的有效波高精度与L2级有效波高精度的比较结果,当有Hs<1 m时,三者的精度接近,随着有效波高的增加,加权最小二乘重跟踪的有效波高精度较最小二乘重跟踪的精度有明显提高,其精度提高比统计如表2所示,与L2级有效波高数据相比,整体精度优于L2级数据的有效波高。
表2 加权最小二乘重跟踪有效波高精度提高统计
图6 有效波高精度随有效波高的变化
图7 后向散射系数精度随有效波高的变化
另一方面,对比图5、图6和 图7,可以看出数据集中Hs>8 m的样本数远少于Hs<8 m的样本数,这是由于海洋上风、浪较大的天气并不常见,对于Hs>8m的情况,其参数估计精度需要更长时间的观测数据,才能更精确描述重跟踪参数精度和L2数据集参数精度的大小关系。