基于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 假彩色圖像，為接下來的海冰分析提供更加有效的支撐。