基于Sentinel-1双极化数据的北极海域假彩色图像合成方法

于皓，田忠翔，李春花

（国家海洋环境预报中心自然资源部海洋灾害预报技术重点实验室，北京 100081）

1 引言

北极是对全球气候变化响应的敏感地区之一，海冰的快速变化影响着大气和海洋之间的能量平衡和物质平衡。自1979年有卫星观测以来，北极海冰呈快速减少的趋势[1-4]，这也使北极航道的开通成为可能[5]。随着北极航道的开通利用，船舶航行的安全保障日益重要。海冰的存在严重威胁到北极航道通航船只的安全，尤其是无冰级的船舶，因此，实时获取北极航道区域精细化的海冰分布对于航行船舶至关重要。自2013年以来，我国中远海运特运集团夏季在北极东北航道共通航了56艘次，而且有逐年增多的趋势[6]。使用高分辨率合成孔径雷达（Synthetic Aperture Radar，SAR）数据开展海冰的识别和分类[7]，并及时、全面和细致地获取北极航道海域不同类型海冰的分布，具有非常重要的科研意义和现实意义。

目前，天基遥感卫星平台搭载的用于海冰监测的传感器主要有可见光传感器和微波传感器两类。可见光传感器是一种被动传感器，其图像特征丰富且解译难度低，但容易受云和天气等影响，无法实现全天候海冰监测。微波传感器包括被动微波传感器和主动微波传感器两类。用于海冰监测的被动微波传感器主要是辐射计，它可以被动地接收自然状况下地面反射或发射的微波，不能直接获取海冰图像。用于海冰监测的主动微波传感器多为搭载在卫星上的各种高度计和雷达等，SAR 是其中之一。与其他方法相比，SAR 具有不受天气和云的影响、可以全天候监测海冰变化的特点。哨兵一号（Sentinel - 1）卫星是由欧洲委员会（European Commission，EC）和欧洲航空局（European Space Agency，ESA）共同研制而成的地球观测卫星。Sentinel-1 的超宽幅模式（Extra-Wide swath，EW）采用步进的条带扫描方式（Terrain Observation with Progressive scans SAR，TOPSAR）获取5个子条带的图像，然后拼接成一景图像。该产品已经被广泛应用于极地海冰的观测、监测以及海冰分类中[7-17]。在TOPSAR 模式下，通过逐级扫描观测得到的图像会受到热噪声的影响，尤其是在HV 极化图像中[18-19]。HV 极化图像中的热噪声是一种加性噪声。ESA 的热噪声去除算法与前人的处理方法相似，即假设热噪声功率在空间上具有恒定不变性，在原始SAR 图像中减去重构的热噪声场来去除热噪声[20-21]。然而，实际的热噪声功率在方位向和距离向是变化的[18]，因此，使用ESA 提供的热噪声去除算法处理后的HV 极化图像仍然包含较强的噪声，不能达到理想的效果[18]，PARK 等[18-19]通过计算后向散射系数的最优缩放因子和平衡因子，结合局部信噪比和NESD（Noise Equal Standard Deviation）模型，利用滑动窗口对局部后向散射系数进行缩放，并通过噪声功率补偿取得较好的热噪声去除效果。SUN等[20]在PARK 等[19]提出的去噪算法基础上改进了最优缩放因子和平衡因子的求解过程，将计算效率提高了10倍以上，但该方法牺牲了部分扇贝噪声的去除效果。HH 极化图像对入射角更加敏感，入射角的改变会导致后向散射系数发生变化，使得HH 极化图像会沿距离向出现显著的亮度变化。研究证明，入射角变化引起的后向散射系数的变化会在海冰分类过程中引入误差，而对入射角进行归一化后，可以显著提高海冰分类精度[16]。因此，在开展海冰识别和分类任务之前，对HH极化图像进行入射角校正是非常必要的。入射角校正一般有两种思路：一种是人工选取窗口采集HH 极化图像中不同类型海冰的后向散射系数，与入射角建立线性回归方程[12,14,22-24]，但是这种方法获取的校正方程针对性较强，往往具有较强的局地性和时间局限性，不适合批量处理图像；另一种是使用整张HH极化图像建立后向散射系数与入射角的线性回归方程[25]，这种方法适用于对不同时间、不同区域获取的大量HH极化图像进行入射角校正，但会损失一部分校正精度。

近些年来，出现了很多基于Sentinel-1影像的自动和半自动海冰分类算法，包括贝叶斯分类器[16]、支持向量机分类器[17]、神经网络[14-15]和随机森林[10]等。然而，单极化图像包含的地物信息非常有限，同时使用HH 和HV 极化图像可以显著地提高自动分类算法的准确度[10,18]。BOULZE 等[15]证明通过叠加HH 和HV 极化图像，利用卷积神经网络（Convolutional Neural Networks，CNN）可以开展多种类型海冰的分类任务，但是该方法没有考虑入射角的影响，而且简单叠加两层灰度图像无法将合成后的产品进行有效可视化，不能很好地分析影响海冰分类精度的因素。ESA 提供了一种RGB 假彩色图像的制作方法（来源：https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-1/sar-ice/.），但由于缺少入射角校正和有效的热噪声去除以及通用的可视化阈值，导致得到的RGB假彩色图像质量较差。因此，本文基于Sentinel-1 双极化数据，在现有的去噪算法基础上，进一步使用图像处理技术减弱残余噪声对图像的影响，尤其是入射角的影响；增强图像信息，优化数据处理流程，避免人工干预引入误差，显著提高了不同类型海冰的特征信息，合成的高质量的RGB假彩色图像为提高海冰分类结果的精度奠定了良好的基础。

2 数据与方法

2.1 Sentinel-1数据

Sentinel-1 包含两个具有相同轨道的极轨卫星Sentinel-1A 和Sentinel-1B。两个卫星先后发射于2014 年 4 月 3 日和 2016 年 4 月 25 日，其中，Sentinel-1B 卫星由于技术性故障已于2021 年12 月23 日停止工作 。Sentinel-1 搭载的 C 波 段 SAR 拥有 HH 和HV 两种极化方式以及4 种成像模式，即条带模式（Stripmap，SM）、EW、宽幅干涉模式（Interferometric Wide swath，IW）和波模式（Wave mode，WV），其中，EW 模式数据的方位向分辨率和距离向分辨率分别为40 m 和25 m。Sentinel-1 卫星不但具有全天时和全天候观测的特点，而且具有及时、可靠、重访周期短和覆盖范围广等方面的优势，其中，Sentienl-1 EW模式数据多适用于极地海冰监测[8]。因此，本文使用的数据为EW 模式下的GRD（Ground Range Detected）产品。在进行RGB 假彩色图像质量检验中，本文使用了2017—2018年西北航道和斯瓦尔巴群岛附近的122景包含不同地物的Sentinel-1影像。

2.2 RGB假彩色图像质量检验方法

为了客观和定量地分析最终获取的RGB 假彩色图像的可靠性，本文选取结构相似性指数（Structural Similarity Index Measure，SSIM）作为判别依据。SSIM 是一种定量度量增强图像和原始图像之间相似度的指数，取值范围为[0，1]，SSIM 值越大意味着两个图像的结构相似度越高[26]。也就是说，当SSIM 值接近1 时，表明增强图像与原始图像具有高度相似性。SSIM可以根据式（1）进行计算：

式中：X、Y是原始图像和重建图像的两个相同尺寸窗口的灰度值；α、β、γ是这些参数的权重；l、s、c是亮度、结构和对比度，分别按照式（2）—（4）来计算：

式中：μX是X的平均值；μY是Y的平均值是Y的方差是X的方差；σXY是X、Y的协方差；c1和c2是稳定分母的两个参数

由于图像的统计特征和失真情况存在空间分布不均的情况，因此需要使用滑动窗口对大尺寸图像进行分块[27]。考虑到窗口形状对分块的影响，本文采用高斯加权计算每一个窗口的均值、方差以及协方差，然后计算对应块结构的SSIM 值，最后将整幅图像的平均结构相似性（Mean SSIM，MSSIM）作为两图像的结构相似性进行度量。前人研究发现[15]，基于Sentinel-1 数据开展多类型海冰分类时，使用大小为50×50 像元的滑动窗口可以得到最佳海冰分类结果，而MSSIM 算法要求滑动窗口尺寸须为奇数，因此本文所选窗口大小为49×49像元。

3 Sentinel-1数据预处理

在进行RGB假彩色图像制作前，为了使得到的RGB 假彩色图像达到最佳效果，需要对Sentinel-1数据进行预处理。数据预处理包括HH 极化数据的入射角校正和HV极化数据的热噪声去除。

3.1 入射角校正

SAR 的侧视成像方式会导致Sentinel-1 HH 极化数据中的后向散射系数随着入射角的变化而变化。一般来说，后向散射系数σ0会随入射角θ的增大而减小。由于入射角变化引起的后向散射系数的变化会在海冰分类算法中引入显著的误差[10]，因此需要对HH 极化数据进行入射角校正。为了尽可能扩大RGB假彩色图像合成方法的适用范围，本文根据北极斯瓦尔巴群岛附近和北极高纬地区16 景Sentinel-1 EW 影像获得的回归模型对所有Sentinel-1影像进行入射角校正[23]。这16景影像均获取自海冰冻结期，且覆盖了多种类型的海冰、开阔水域以及全部变化范围的入射角。入射角校正公式可以表达为：

式中：σnew为校正后的后向散射系数，单位是dB；θ为入射角，min(θ)是入射角的最小值，单位是°。经过线性回归模型处理后的HH 极化σ0值转为dB 值后近似正态分布[21]，可以更有效地去除极值点。

3.2 热噪声去除

虽然Sentinel-1 HV 极化数据对入射角不敏感，但存在很强的热噪声（见图1a），因此对HV 极化数据的预处理主要是去除热噪声。Sentinel-1 HV极化数据在低后向散射区域的热噪声非常显著。HV 极化数据中的热噪声由两部分组成，一部分为距离向的由条带间功率不平衡造成的噪声阶跃，另一部分为方位向的扇贝效应。由于ESA 提供的热噪声去除算法不能有效地去除热噪声，本文综合PARK等[18-19]提出的热噪声去除算法进行HV 极化数据的热噪声去除。HV 极化数据中的热噪声是一种加性噪声，可以用式（6）描述：

式中：PS是信号强度；PN是噪声强度；G是 SAR 信号的总增益系数。热噪声去除过程可以简单地描述为在给定的信噪联合功率PSN中减去噪声功率G·PN。为了得到准确的噪声功率G·PN，热噪声去除过程包括以下3个步骤：

（1）初步重构噪声场。ESA 将每幅影像的热噪声信息以噪声矢量的格式保存在独立的XML 文档中，读取XML 文档中的噪声矢量信息可以得到PN。通过准确计算每个子条带burst 开始时的零多普勒时间，可以得到每条方位线对应的准确时间，利用天线方位向的增益方向图能为每条方位线找到正确的增益系数G。将PN和G相乘可以得到初步的噪声场（见图1b）。

（2）噪声场条带间功率平衡。由于ESA 提供的噪声向量与真实噪声存在误差，直接使用ESA 提供的噪声矢量进行热噪声去除会在距离向出现信号阶跃现象。为了消除各子条带间的噪声阶跃现象，可以假设去噪过程模型满足线性关系[18]：

式中：s(k)是去除热噪声后的σ0值是未经校正的原始σ0值是利用ESA 提供的热噪声矢量经过双线性插值计算得到的σ0；Kns,n是最佳噪声比例因子是条带间噪声功率平衡因子，n为子条带数，n={1，…，5}。Kns,n可以利用大量HV 极化数据通过加权最小二乘法求解得到:

式中：αn和βn分别是不同子条带线性模型的斜率和截距；i是条带间边界处距离向的像元数；n={2，3，4，5}。由于只有4个条带间边界，所以设置为0。条带间功率平衡后的噪声场见图1c。

（3）局部残余噪声功率补偿。由于XML 文档中记录的噪声矢量信息存在误差，原始影像减去利用上述方法获取的热噪声后，会导致一些像元点的σ0值变为负值。为了消除部分负噪声功率的影响，需要对负噪声功率进行噪声补偿。首先，定义信噪比（Signal-Noise Ratio，SNR）为经过高斯滤波器后的σ0值（s0g）与重构噪声场（NESZ）的比值，即：

然后，进一步计算得到功率补偿后的σ0[19]：

式中：s0o为局部残余噪声功率补偿后的σ0；s0offset为噪声场补偿值，可以取重构噪声场的平均值。最终我们可以得到去除热噪声的HV 极化灰度图像（见图1d）。

图1 HV极化数据热噪声去除过程（图像取自斯瓦尔巴群岛附近，2018年1月25日04:11:15 UTC）Fig.1 Thermal noise removal process of HV polarized data（Image is acquired from the region near Svalbard at 04:11:15（UTC））on 25 January 2018

4 RGB假彩色图像制作

本文以经过预处理后的HH 和HV 极化数据为基础，经过RGB假彩色图像颜色通道填充和两次直方图均衡化，得到最终的RGB假彩色图像（见图2）。

图2 RGB假彩色图像制作流程图Fig.2 Flow chart of synthesizing RGB false color image

4.1 RGB原始图层数据

由于Sentinel-1 EW 模式中的海冰信息主要包含在HH 极化数据中[29]，为了更贴近地物的真实颜色，在进行颜色通道选择时蓝色通道使用HH 极化图像，红色通道则使用HV 极化图像。为了减小颗粒噪声对图像质量的影响，需要对HH和HV极化数据按照式（13）和式（14）进行偏移量处理：

绿色通道图像是将经过偏移量处理后的HH 和HV极化数据按照式（15）进行组合制作而成：

以该方式组合得到的RGB 假彩色图像可以最大限度地还原原始地貌的颜色特征。之后，使用ESA 提供的参考数据范围[22]对数据进行标准化处理，其中红色通道数据最小值为0.02，最大值为0.10；蓝色通道数据最小值为0，最大值为0.32；绿色通道最小值为0，最大值为0.6。最后，对经过标准化处理后的数据按照式（16）进行gamma 值为1.1 的gamma校正，得到最终的原始图层数据：

式中：σ和σγ分别是各颜色通道在 gamma 校正前和校正后的数据。

4.2 两次直方图均衡化

虽然我们在制作原始图层数据时进行了数据标准化，但在将其转化为灰度图像时，如果原始图层数据离散程度过大，得到的灰度图像的对比度会偏低。为了较好地减弱这种影响，本文使用了以下方法：首先，将原始图层数据的σ0转化为dB 值，此时，dB值近似呈正态分布；然后，对每个颜色通道去除数据两端2.5%的极值点，并转化为灰度图像；最后，对每个颜色通道的灰度图像进行全局直方图均衡化（Histogram Equalization，HE）处理，并按照颜色图层拼接后得到第一次HE后的RGB假彩色图像（见图3a—c）。全局HE是一种常规的图像处理方法，具体原理和步骤可以参考PIZER等[28]，本文不再赘述。

全局HE 不仅可以有效地改善数据离散程度过大导致的图像对比度过低现象（见图3），而且可以避免选取有效阈值时的人工干预，使批量制作数据集更加可靠、简单。但受限于有限的不连续灰度级数，在某些包含地物类型较多的HH 极化图像中，灰度分布直方图的变化比较剧烈，甚至出现某些灰度区间不存在像素点的情况，这会造成累积概率密度函数发生剧烈变化，从而导致图像的直方图出现较强的不均匀性，使得图像中的部分信息不能被正确地显示（见图4a和4c）。另外，由于入射角校正不能完全去除入射角对σ0的影响，使得RGB 假彩色图像中蓝色通道图像的亮度不均匀。这不但会导致灰度图像的可识别特征减少[30]，还会使合成的RGB假彩色图像色调不平衡，出现近轨侧偏蓝、远轨侧偏红的现象（见图3c）。为了减弱这种影响，本文进一步对每个颜色通道全局HE 之后的灰度图像进行限制对比度直方图均衡化（Contrast Limited Adaptive Histogram Equalization，CLAHE）[29]，即第二次直方图均衡化。CLAHE 是一种对比度受限情况下的自适应HE 算法，目前多应用于医学影像处理领域[30]。CLAHE首先需要确定两个参数，即直方图裁剪阈值（ClipLimit）和子图像块窗口大小（Block Size）。当图像离散熵函数曲率最大时，CLAHE 处理后的图像包含的地物信息最多[31]。本文将CLAHE 处理后图像离散熵函数曲率最大作为确定ClipLimit 和Block Size 两个参数的依据，即选取离散熵函数曲率最大时对应的两者的值作为最终确定的参数。本文选取的ClipLimit 和Block Size分别为3 和5；然后，使用大小为Block Size 的窗口将全局HE 之后的灰度图像裁剪为若干子图像块，计算每个子图像块的概率密度直方图，将概率密度超过ClipLimit的部分平均分配到各个灰度级上；最后，使用双线性插值加快CLAHE 速度，同时消除相邻子图像块之间不连续的边界。

图3 全局HE之后的HH和HV极化图像和相应的RGB假彩色图像（该图像与图1使用的图像相同）Fig.3 HH and HV polarized images after HE and the corresponding RGB false color images（This image is the same as the one used in Fig.1）

通过CLAHE 可以有效地纠正图像中块状不连续的缺陷，改进灰度分布直方图剧烈变化的情况（见图4），进而减弱RGB假彩色图像的色调不平衡，显著改善近轨侧偏蓝、远轨侧偏红的现象（见图3d—f）。

图4 CLAHE前后的HH极化图像和对应的灰度直方图（该图像与图1使用的图像相同）Fig.4 HH polarized images and the corresponding grayscale histograms before and after CLAHE（This image is the same as the one used in Fig.1）

5 RGB假彩色图像质量检验

RGB 假彩色图像的制作主要包括数据预处理、全局HE和CLAHE 3个步骤。数据预处理部分使用的入射角校正方法比较简单，不能完全去除入射角对HH 极化图像的影响[25]，因此，为了进一步改善HH 极化图像质量，本文引入图像处理中的CLAHE方法。然而，CLAHE 在将图像进行小窗口化的局部HE 过程中，不可避免地会突出小窗口中包含的纹理特征。由于本文应用的热噪声去除算法不能有效去除所有HV 极化图像中的热噪声，导致原本不太明显的残余热噪声会在CLAHE 后出现不同程度的增强，造成RGB假彩色图像质量变差。为了定量评估最终获取的RGB假彩色图像的质量，本文将全局HE 得到的RGB 假彩色图像作为原始图像，将CLAHE 得到的RGB 假彩色图像作为增强图像，将计算了122幅RGB假彩色图像的MSSIM值。

图5 给出了不同MSSIM 值的RGB 假彩色图像。从图中可以清楚地看到在MSSIM 值小于0.7的图像中存在一定的残余热噪声，且残余热噪声在CLAHE 后被放大。红色不连续热噪声条纹会造成整幅图像偏红，导致原有地物色调发生变化，地物可辨识度显著降低。当MSSIM 值大于0.7 时，RGB假彩色图像质量显著提高，改善了近轨侧偏蓝、远轨侧偏红的情况，增强了地物的可辨识度。

图5 不同MSSIM值的RGB假彩色图像Fig.5 RGB false color images with different MSSIM values

图 6 给出了 122 幅 RGB 假彩色图像 MSSIM 值的频率分布直方图，其中，MSSIM 平均值为0.82，中位数为0.84。经统计，约有87.70%（107 幅）图像的MSSIM 值大于0.7。也就是说，上述处理方法对大多数Sentinel-1 影像是适用的。图7 给出了本文获取的不同类型海冰的RGB假彩色图像，将此图像与ESA 的方法获得的图像进行对比。从图中可以看到，本文获取的RGB 假彩色图像对多年冰（见图7a）、一年冰（见图7b）以及新冰和开阔水域（见图7c）的特征都表现出不同程度地增强，而且从中能够清晰地分辨出水道区域生成的新冰（见图7b）。ESA 的方法获取的RGB 假彩色图像更容易受入射角和热噪声的影响，不同地物信息特征差异较小，容易错判海冰/海水类型。由此可见，本文提出的RGB 假彩色图像合成方法可以明显提高地物信息特征，同时，相较于单纯地将HH 和HV 极化图像叠加在一起进行海冰识别和分类，本文获取的RGB假彩色图像引入了色调变化，可以将增强后的地物信息可视化，增加了Sentinel-1 影像的可读性，为进一步开展海冰识别和分类奠定良好的基础。

图6 MSSIM值概率分布直方图Fig.6 Probability histogram of MSSIM

图7 本文方法与ESA方法[26]合成RGB假彩色图像的对比图Fig.7 Comparison between RGB false color images obtained using the method proposed in this paper and the method provided by ESA[26]

6 结论与展望

基于前人的研究结果，本文提出了一种新的Sentinel-1双极化数据假彩色图像合成方法，该方法包括数据预处理和RGB假彩色图像制作两个部分。在数据预处理阶段，本文使用线性回归模型对HH极化数据进行入射角校正，并使用一种更高效的热噪声去除算法进行HV 极化数据热噪声的去除。在RGB 假彩色图像制作阶段，以经过数据预处理后的HH 和 HV 极化数据为基础，通过混合 HH 和 HV 极化数据得到新的颜色图层，并使用全局直方图均衡化和限制对比度直方图均衡化进行图像增强，获得最终的RGB 假彩色图像。两次直方图均衡化可以避免选取可视化数据范围时的人工干预，同时降低入射角对RGB 假彩色图像质量的影响。图像质量定量评估结果表明，本文提出的假彩色图像合成方法能够有效去除热噪声和入射角对图像质量影响，引入的色调变化增强了图像携带的地物特征信息，能够将增强后的图像可视化，可为海冰信息的提取奠定良好的基础。

在不同影像中，入射角差异和热噪声带来的影响有所差别，会导致限制对比度直方图均衡化在减弱入射角影响的同时，放大某些影像中的残余热噪声，进而造成一部分RGB假彩色图像质量较差。因此，对不同地物、不同时间的Sentinel-1 影像分别进行入射角校正，使用更加有效的热噪声去除算法减小HH和HV极化图像中的残余噪声，或者使用暗角校正等图像处理方法进一步提高HH 极化图像的质量，将有助于我们获取质量更好的RGB 假彩色图像，为接下来的海冰分析提供更加有效的支撑。