基于MODIS 數據的渤海海冰厚度反演算法優化

朱星源，蘇潔,，宋梅，楊茜,，梁韻,

(1.中國海洋大學 海洋與大氣學院，山東 青島 266100；2.中國海洋大學 物理海洋教育部重點實驗室，山東 青島 266100；3.中國科學院西北生態環境資源研究院 冰凍圈科學國家重點實驗室，甘肅 蘭州 730000；4.中國科學院南海海洋研究所熱帶海洋環境國家重點實驗室，廣東 廣州 510301)

1 引言

渤海位于37°～41°N 之間，由遼東灣、渤海灣、萊州灣和中央淺海盆地組成。三面被陸地包圍，屬于半封閉海區，與外界海水的熱量交換較少，在冬季西伯利亞高壓、太平洋副熱帶高壓等因子的影響下，結冰現象顯著[1?4]。入冬之后，隨著負積溫的累積和冷空氣的侵襲，特別是強寒潮的暴發和延續，海冰面積會不斷擴大，厚度也會隨之增加[5?8]。海冰會對海上的交通運輸和生產活動產生影響，甚至會引發自然災害，如1969 年渤海大冰封，持續的冰封不僅導致海上航行受阻，“海二井”石油平臺還直接被海冰推倒[9]；2009–2010 年1 月中下旬的渤海冰情為近40 年同期最重，冰情等級為4 級，對環渤海地區產生了嚴重影響，直接經濟損失達到了63.18 億元[10]；近10 年渤海最重冰情在2012–2013 年冬季，冰情等級為3.5 級，水產養殖受災面積超過2.292×104hm2，直接經濟損失達3.22 億元[11]。為了了解渤海海冰的變化趨勢，研究海冰災害的成因并進行預測，對渤海海冰的準實時監測是非常必要的。

渤海海區范圍較小，微波數據的分辨率較低，不足以進行有效的監測，采用更高分辨率的可見光數據進行海冰反演不失為一種合理有效的方法。渤海海冰密集度可見光數據反演的研究已有不少[12–16]，但作為海上冰情重要指標的海冰厚度，由于反演參數的獲取較為困難，進展相對緩慢。Grenfell[17]提出的海冰厚度與反照率呈現指數關系的計算公式是目前較常用的可見光冰厚反演算法；謝鋒等[18]將此算法應用于AVHRR（Advanced Very High Resolution Radiometer）數據，對遼東灣海域冰厚進行反演；Yuan 等[19]根據不同區域海冰光譜特性將渤海劃分為5 個區域，也采用此算法基于AVHRR 數據分區域進行了渤海海冰厚度的計算；Su 和Wang[20]基于Grenfell[17]和謝鋒等[18]的結論，使用MODIS（Moderate-resolution Imaging Spectroradiometer）數據得出渤海冰厚計算值，但并未進行與實測數據的定性分析驗證；Liu 等[21]則將該算法應用于GOCI（Geostationary Ocean Color Imager）衛星數據，反演渤海海冰厚度。在這類反照率與海冰厚度指數關系模型的應用中發現，反照率不僅與海冰厚度有關，也受到海水中泥沙懸浮物的影響[22]，然而渤海海區受黃河泥沙輸入影響，再加上渤海海區內水動力環境相當復雜，導致泥沙懸浮物濃度在時間和空間尺度上變化劇烈[23]，這也影響到了反照率與海冰厚度指數關系算法的準確性[18–19,24]。除海冰厚度與反照率指數關系的算法外，吳龍濤等[14]簡單地使用MODIS 不同波段的分段線性法計算冰厚，其結果與石油平臺獲得的冰厚相比偏厚；Ning 等[25]通過海冰在MODIS 和TM（Thematic Mapper）衛星不同波段的光譜特性對海冰進行分類，進而大致估算海冰厚度；Yuan 等[26–27]從海冰的光輻射傳輸過程出發，提出了一種光學遙感海冰厚度半經驗模型，該模型考慮到了其中各參數的時空異質性，并利用MODIS 數據多通道的反射率反演海冰厚度，但是存在一定局限性，容易低估高反射率區域的海冰厚度。

冰水分離是反演海冰厚度的重要步驟，近年來人們做了大量的研究。吳奎橋等[28]結合了MODIS 通道1 反射率和通道32 亮溫進行冰水識別；Su 和Wang[20]使用雙通道比值閾值判別法識別海冰；Su 等[29]利用灰度共生矩陣紋理分析方法進行渤海海冰探測；Zhang 等[30]介紹了一種分類與回歸樹的方法檢測海冰；Su 等[31]提取了海冰的溫度特征和紋理特征，然后使用支持向量機進行海冰范圍估計；由于渤海是我國的內海，黃河口附近海區注入大量泥沙、水體較為渾濁，這導致目前常用的方法大多存在懸浮泥沙誤判的問題；通過海冰紋理識別冰區的方法雖然能排除高濃度泥沙海區，但是對于平整薄冰的檢測能力不足。Li和Yang[32]提出了一種基于多種海冰特征的線性分解法用于冰水分離，在提取海冰形狀特征時，創造性地使用了Canny 算子對海冰的邊緣和裂縫進行識別來區分泥沙與海冰[33]，為提取海冰范圍提供了一種新思路，遺憾的是該方法沒有充分完善Canny 算子的提取結果，不同水色海區交界區域會被誤判。另外，通過線性分解模型提取海冰范圍也存在不易提取平整薄冰等問題。

綜上所述，這些研究為獲得渤海海冰厚度可見光遙感數據奠定了一定的基礎，但仍存在一些問題：在冰水分離方面，Canny 算子識別海冰邊緣和裂縫作為一種提取海冰范圍新思路，使用時仍需解決不同水色海區交界區域被誤判和不易提取裂縫較少的平整薄冰區等問題；在冰厚計算方面，泥沙懸浮物提高了海水反照率，影響了反照率與海冰厚度指數關系算法的準確性。為此，本文將針對Canny 算子提取海冰時存在的問題進行改進和完善，構造出一套基于Canny算子的完整的自動化海冰范圍提取方法；而對于海冰厚度計算，為了降低泥沙影響和提高計算結果準確性，本文將基于渤海海區的物理特征，通過試驗確定海冰厚度與反照率指數關系算法中的相關參數，實現對算法的改進，并且將最終的反演結果與平臺實測數據進行對比和分析。

2 數據來源

2.1 MODIS 遙感數據

海冰厚度反演使用的數據為1 000 m 分辨率的MODIS L1B 數據和MODIS 03 地理信息數據，由美國宇航局提供并且已經過幾何校正。中分辨率成像光譜儀MODIS 搭載在TERRA 和AQUA 衛星上，可以覆蓋從可見光到近紅外（0.405～14.385 μm）的36 個通道，最高空間分辨率可達250 m[34]，可以在每日上午、下午各獲得一次渤海海區觀測數據，具有免費獲取、時空分辨率和光譜分辨率高等特點，使其可以有效地進行渤海海冰的監測。

本文主要使用的MODIS 通道和相對應的光譜范圍如表1 所示，其中波段1、3 和4 的輻射率用于獲取真彩圖和灰度圖，波段1、6 的反射率用于云剔除，波段31、32 的亮溫用于冰水判別，波段1、2、3、4、5 和7 的反射率用于計算海冰厚度。數據預處理包括輻射定標和太陽高度角訂正。為了避免天氣原因造成的影響，本文選取的MODIS 數據均在晴空條件下。

表1 MODIS 部分波段光譜范圍Table 1 Spectral range of MODIS partial bands

2.2 渤海海上石油平臺實測數據

本文采用的實測數據為渤海海上石油平臺的觀測數據，包括2013–2014 年冬季和2015–2016 年冬季的JZ20-2 平臺（40.500°N,121.352°E），2020–2021 年冬季的JZ20-2 和JZ9-3 平臺（40.664°N,121.462°E）冰厚觀測數據，觀測時間為當地時間早上8 點至下午7 點，每小時進行一次海冰最大厚度和平均厚度觀測；2009–2010 年冬季和2012–2013 年冬季數據引自Zeng等[24]和Karvonen 等[35]的文章，來源于渤海海上石油平臺JX1-1（39.977 5°N,121.058°E）、JZ25-1S（40.243°N,121.025°E）、JZ20-2 和JZ9-3 平臺的海冰厚度現場觀測數據，測量方法為目視觀測，測量時間為當地時間下午6 點。平臺具體位置如圖1 所示。本文在使用實測數據時，將全部實測數據隨機分為訓練數據集和測試數據集，分別占總數據量的63%和37%，訓練數據集用于海冰厚度算法的改進，測試數據集用于驗證和分析算法改進效果。

圖1 獲取實測數據的石油平臺位置Fig.1 Positions of oil platform observing the measured data

3 反演算法及改進

本文首先基于MODIS 數據中云的光譜特征進行云剔除；繼而使用Canny 算子提取海冰裂縫和邊緣，再通過進一步處理實現冰水分離；在計算海冰厚度時對海冰厚度與反照率指數關系模型中參數設置進行了分析和優化，使其更加符合渤海海區的水文特征，其中，將海水反照率參數由固定值變為隨海域實際情況而改變的動態數值，利用渤海海上石油平臺實測數據獲取了更為準確的海冰衰減系數估計值。具體反演流程及改進情況如圖2 所示。

圖2 渤海海冰厚度反演算法流程Fig.2 Flow chart of Bohai Sea ice thickness retrieval algorithm

3.1 云剔除

本文選取渤海冰區晴空下的MODIS 數據進行反演，而水區上空的云在冰水分離階段容易被誤判為冰，還會影響算法中需要的海水反照率的計算，因此，首先要對渤海海區進行云剔除。云和海冰在可見光波段都有較高的反射率，而海冰在近紅外波段反射率顯著下降，因此本文采用的云剔除方法是：構建band 1和 b and 6反 射率r的歸一化指數R1_6[36]，

對于式（1）結果的頻數分布曲線，Bayes 分類準則認為將云峰和水冰峰之間的谷值作為閾值[37]，誤判概率最小，剔除低于閾值的像元。

3.2 冰水判別方法的改進

3.2.1 海冰范圍提取原理

冰區有不規則裂縫和邊緣，而水區表面平滑，這是冰區和水區之間表面特征的重要差異，Canny 邊緣檢測算法可以提取冰區的邊緣以及內部裂縫[32]，其重要優勢是雙閾值檢測，可以有效地減少渤海清、濁水的邊緣誤判。本文利用Canny 算子的優勢，針對引言中分析的算法缺陷，加入一系列處理步驟，包括真彩圖灰度化、高斯模糊、空洞填充、灰度圖二值化和表面溫度（ST）自動化閾值判別法，提出一個可以自動化進行冰水判別并提取渤海海冰范圍的新方案。

以2010 年1 月22 日的MODIS L1B 渤海海區數據為例（圖3），具體分析冰水分離過程。

圖3 2010 年1 月22 日渤海MODIS 真彩圖Fig.3 MODIS true color image of the Bohai Sea on January 22,2010

首先將真彩圖進行灰度化（圖4a），使其可以進行基于單色圖像的Canny 算子裂縫和邊緣提取，提取效果如圖4b 所示。由圖可見，絕大多數的海冰被密集的裂縫覆蓋，但是仍有部分平整薄冰裂縫較為稀疏，如果只是提取裂縫密集的區域，則會遺漏部分平整薄冰，因此需要分別采用不同方式提取裂縫密集冰和平整薄冰；同時，部分裂縫較為密集的區域位于清水濁水交界處或者海冰邊緣外側，不屬于冰區；更重要的一點是圖4b 由線狀的裂縫組成，無法對裂縫密集的區域直接提取，綜上所述，進行海冰提取還需要解決3 個問題：（1）通過線變面將密集的裂縫形成面狀區域，使裂縫密集區作為區域得以直接提取；（2）通過缺失冰區填充得以提取裂縫較為稀疏的薄冰區；（3）去除非冰區的裂縫密集區。

3.2.2 提取裂縫密集區

第一步，將裂縫密集區作為區域進行整體提取。Li 和Yang[32]通過計算以每個像元為中心的圓形區域內裂縫的密度來給每個像元賦值，從而將線圖變為面圖，但是這種方法最終給出的圖像較為模糊，不易區分裂縫密集區的邊界。本文采用的方法是對圖像進行高斯模糊，高斯模糊的優勢在于增大了距離中心像元較近像元的權重，中心像元的權重最高，因此裂縫密集區邊界更為清晰。對圖4b 進行高斯模糊處理后，得到圖4c，可以看到裂縫密集的區域亮度較高，再對圖4c 進行二值化處理即可提取裂縫密集區域（圖4d）。

3.2.3 提取裂縫稀疏薄冰區

第二步，提取裂縫較為稀疏的薄冰區。由于Canny算子對海冰邊緣更為敏感，導致容易忽視海冰內部半封閉或封閉的平整薄冰，如圖4d 的紅圈區域。為了避免出現這種情況，我們具體做法為進行一次膨脹算法，將半封閉缺失薄冰變為封閉缺失薄冰，二值圖空洞填充后進行一次腐蝕算法，即可完成對內部缺失薄冰的再提取，結果如圖4e 所示。

3.2.4 誤判像元去除

第三步，去除誤判為冰的水像元。由圖4e 可見提取到的區域存在水像元并且邊緣過于光滑，不利于3.3 節改進冰厚計算方法時所需的海冰外緣線的判斷，所以要對誤提取的部分進行修正。圖4f 為基于圖4e 的裂縫區灰度圖，由圖可見水像元的灰度值明顯較低，可以使用最大類間法進行低灰度像元去除，而最大類間法在對像元進行分類時，傾向于在類內方差較大的那一類像元（海冰像元）中選擇閾值，為了避免造成海冰范圍損失，采用約束灰度范圍的改進型最大類間方法[38]。結果如圖4g 所示，消除了部分水像元而且冰區的邊緣更加符合對海冰外緣線的人工視覺判斷。

3.2.5 表面溫度自動化閾值判別法

解決了以上3 個問題后，冰區（圖4g 白色區域）仍包括了一些泥沙濃度較高的濁水或清–濁水交界處的像元，將其去除才能完成海冰范圍提取。海冰和海水在溫度上的差異可以作為一個重要的判據[32,39]，可以有效地修正冰水分離時泥沙懸浮物造成的海冰誤判，使用MODIS 數據31、32 波段亮溫計算出圖4g 中白色區域的表面溫度（ST）如圖4h[40]，由于海冰ST 往往低于海水ST，因此可以通過設置合適的閾值，剔除ST 高于閾值的部分。

關于閾值的取值，本文參考強度比算法的思路，提出了基于ST 頻數比例的提取方法。強度比算法最早用來處理可見光航拍圖像，以每個像元與相鄰像元灰度值的差值大于臨界值作為頻數統計的判斷依據[41]，但是此法難以適用于渤海這種情況較為復雜的大范圍區域。為了實現對閾值進行較為準確自動化選取，本文的算法如下：

以2010 年1 月22 日的數據為例，（1）對整個渤海海區的ST 值域進行頻數統計，記為 ?(k)，k為份數，間隔為0.02 K。（2）對圖4h 的識別為冰區像元的ST 進行頻數統計（圖5a），記為 δ(k)。（3）最后計算二者的比值χ(k)=δ(k)/?(k)，將 χ (k)稱為粗糙度，畫出粗糙度曲線如圖5b 所示。由于海冰在之前的一系列提取流程中被保留下來，表現為低溫區的粗糙度較高而接近1，而隨著ST 升高達到某個值時粗糙度迅速下降，粗糙度迅速下降的區間就是閾值所在的區間，這里，我們取粗糙度為0.4 時對應的ST 作為閾值，所對應的閾值為272.2 K，該閾值隨渤海ST 分布的變化而變化。將ST 小于閾值的像元剔除，最終的冰水分離結果如圖4i 所示，白色部分即為最終提取到的海冰范圍。

圖4 2010 年1 月22 日渤海冰水分離過程Fig.4 Process diagram of ice-water separation of the Bohai Sea on January 22,2010

圖5 頻數分布（a）和粗糙度分布（b）Fig.5 Frequency distribution (a) and roughness index distribution (b)

3.3 海冰厚度算法及改進

Grenfell[17]和 Allison 等[42]在極地調查中發現，當海冰厚度從2 cm 增加到9 cm 時，海冰反照率由0.11 增加到0.24。Grenfell[17]提出反照率與海冰厚度呈指數關系模型為

式中，h為海冰厚度，單位為m；αmax為無限厚冰反照率，一般取0.7；系數為海水反照率；μα為海冰的反照率衰減系數，謝鋒等[18]的估計值為1.209。

Su 和Wang[20]使用上述模型基于MODIS 數據在渤海海區進行海冰厚度計算。將式（2）中的海水反照率 αsea設 為0.06。計算反照率 α(h)時使用Liang[43]提出的反照率 α經驗公式，在MODIS 的6 個波段之間建立 了函數關系為

計算海冰厚度時一般將海水反照率 αsea設為常數[18,20–21]，但渤海的泥沙懸浮物濃度在時間和空間尺度上變化劇烈[23]，將海水反照率 αsea設為固定數值的做法與實際情況明顯不符，也會影響冰厚計算結果的準確性，因此有必要在不同的時間、地點根據實際情況設置不同海水反照率。而冰下的海水由于被海冰覆蓋，無法進行直接觀測。為了解決這個問題，比較可行的辦法是將與海冰相鄰的無冰海區的海水反照率外推到有冰區。下面仍以2010 年1 月22 日渤海數據為例具體分析外推過程。

圖6a 為當日的渤海海區反照率分布，紅線為基于形態學提取的海冰外緣線[44]。為了排除外緣線以外可能存在的冰水混合像元誤判，經試驗，本文將海冰外緣線進一步外推4 個像元得到圖6a 中的藍線，將藍線所在的寬度為3 個像元的狹長條帶區域定義為最鄰近海冰海水區，根據空間相關性理論，距離近的事物關聯更緊密[45]，理論上這個區域內的海水反照率比較接近海冰區內的海水反照率。將最鄰近海冰海水區像元作為插值節點，冰區像元作為被插值點，使用反距離加權插值法對反照率進行插值，結果如圖6b 所示，由圖可見插值后的海水反照率分布基本符合渤海泥沙分布特點。在本研究中使用上述算法將海水反照率 αsea由固定值0.06 變為隨海域實際情況而改變的范圍在0.05～0.13 之間的動態數值。

圖6 渤海海區反照率（a）和海水反照率插值（b）Fig.6 The albedo of Bohai Sea (a) and the albedo of sea water interpolation (b)

謝鋒等[18]在使用AVHRR 數據反演海冰厚度時，對衰減系數 μα的估算方法為：將式（2）變形為式（4），海冰厚度h使用結（融）冰度日法通過實測氣象數據估算[46]，α(h)由衛星觀測給出，αsea和 αmax分別取0.1 和0.7，然后反推出 μα的值，因此這個方法并沒有采用冰厚實測數據參與計算，并且結（融）冰度日法本身存在誤差，同時他們也指出衛星數據來源不同，定標方式不同，μα的取值也不同，再加上本文已經對算法中的αsea做 了優化，不再為常數，所以對 μα進行重新計算是十分有必要的。

為了得到相對準確的衰減系數 μα的估計值，本文隨機抽取部分冰厚實測數據h作為訓練數據集（具體占比見2.2 節），將對應日期的MODIS 相關波段數據代入式（3）得到整個海區的反照率。再通過最鄰近海冰海水區反照率插值算法，獲得實測數據所在位置冰下海水反照率 αsea。然后將冰厚實測數據h，海水反照率 αsea以 及海冰反照率 α(h)代入式（4），最終反推出 μα。

由式（4）可知，衰減系數 μα由比值運算推出，當海冰實測厚度h較小時，分子分母數值均較小，即使是微小的偏差都會導致最終結果產生較大誤差，最終反推出的衰減系數 μα也更容易發散。如圖7 所示，實測厚度較小的薄冰衰減系數 μα較為分散，隨冰厚增大 μα分布越來越集中。實際上，海冰的衰減系數主要由海冰類型，海冰內部鹵水泡、葉綠素等物質的吸收系數及其體積分數決定[17,22]，但是通過衛星無法獲取具體的相關參數，所以在實際應用中一般將衰減系數 μα設為固定值。為了得到衰減系數 μα相對準確的確定值，需要對反推出的 μα進行進一步處理。目前圖7 中所有數據點衰減系數平均值為2.19，為了減小薄冰（實測冰厚小于6 cm）帶來的較大偏差，僅讓厚冰的衰減系數（圖7 中黑點）參與計算，算得其平均值為1.85（圖7中藍色實線），標準差為0.58，然后進行質量控制，僅對處于平均值±標準差以內的數值進行平均（圖7 中兩條紅線之間），最終算得 μα的估計值為1.74（圖7 中藍色虛線）。

圖7 訓練數據集實測冰厚與衰減系數散點圖Fig.7 Scatter plot between attenuation coefficient and measured sea ice thickness of training data sets

3.4 冰厚計算敏感性試驗

本文通過敏感性試驗研究指數關系模型中各參數對反演結果產生的影響。圖8 反映了海冰厚度反演結果對海冰反照率 αice、海水反照率 αsea和衰減系數μα的 敏感性。如圖所示，當海冰反照率 αice與海水反照率 αsea相等時，海冰厚度反演結果為0，隨著海冰反照率 αice的增加，冰厚反演結果呈現指數上升趨勢。在海冰反照率 αice不變的情況下，越低的海水反照率αsea對應越大的冰厚反演結果。那么在泥沙濃度較高的水區（如萊州灣，渤海灣），其實際的海水反照率αsea較高（0.09～0.13），如將其設為以往的固定值0.06，反演的海冰厚度往往偏大，甚至無冰海區，其反演結果也可以達到5 cm 左右。而本文使用的根據實際情況設置海水反照率的做法，會顯著減小高濃度泥沙區海冰厚度的反演結果，對于反照率為0.15 的海冰，計算冰厚時將海水反照率 αsea由0.06 增大到0.1，會使反演結果由8.7 cm 減小到5 cm（μα設為1.74）。另外，本文利用實測數據和MODIS 數據重新估算了渤海海區海冰的反照率衰減系數，取常數為1.74，與之前使用的1.209 相比[18]，冰厚反演結果整體縮小30.5%。

圖8 模型中各參數對反演結果的影響Fig.8 The influence of each parameter in the model on the retrieval result

4 反演結果的對比與驗證

仍以2010 年1 月22 日為例，圖3 為當日渤海海區MODIS 真彩圖。使用本文提出的冰水分離算法提取海冰范圍，然后使用指數關系模型改進前與改進后算法，以及Yuan 等[26]的半經驗模型法進行海冰厚度計算，最終結果如圖9 所示。從冰水分離效果來看，本文所提取海冰范圍與當日真彩圖（圖3）的海冰范圍一致，外緣線也與真彩圖的目測外緣線基本吻合。從海冰厚度反演結果來看，改進后的反演結果（圖9b）與之前相比有所減小，在泥沙懸浮物濃度高的區域尤為明顯，比如渤海灣和萊州灣大面積幾乎透明的薄冰（圖3），改進前算法（圖9a）對這些薄冰的反演結果多處于0.1～0.15 m 之間，而改進后算法結果大多在0～0.05 m 的區間內。再將算法改進前后結果與Yuan等[26]算法結果（圖9c）進行對比分析，空間平均絕對偏差分別為0.045 m 和 0.027 m。

圖9 渤海海冰厚度反演結果對比Fig.9 Comparison of retrieval results of Bohai Sea sea ice thickness

由于與Yuan 等[26]半經驗模型法的相似性高只代表與同類型反演算法的比較情況，無法作為誤差的檢驗標準。因此本文利用渤海海上石油平臺實測數據來檢驗算法改進前后的結果；同時為了檢驗本文對反演算法中海水反照率 αsea與衰減系數 μα的改進效果，分別改變其中一個參數而另一個參數保持不變，將所得結果與實測數據進行比較，4 個試驗的結果見表2。需要指出的是在進行誤差分析時使用的數據集為測試數據集，與3.3 節反推衰減系數 μα時使用的訓練數據集不同。

表2 海冰厚度測試數據集實測數據與反演結果Table 2 Measured data of test data sets and retrieval results of Bohai Sea sea ice thickness

續表 2

如表2 所示，T0（改進前算法厚度）與平均實測厚度之間的平均誤差、平均絕對誤差和均方根誤差分別為6.66 cm、7.05 cm 和8.25 cm；T1（改進后算法厚度）各項誤差分別為0.49 cm、2.74 cm 和3.75 cm，與T0 相比分別降低了93%、61%和67%；其中僅改變海水反照率 αsea的T2 平均誤差、平均絕對誤差和均方根誤差分別降低45%、33%和30%；僅改變衰減系數μα時的T3 誤差降低的幅度比T2 更大，各項誤差分別降低61%、47%和43%。4 個試驗結果與平均實測厚度的相關系數分別為0.434、0.485、0.485 和0.434，其中T0 與T3，T1 與T2 相關系數相同，這是由于衰減系數 μα與海冰厚度反演結果之間呈反比例關系，因此僅衰減系數 μα不同不會影響相關性。我們同時也將算法結果與最大實測海冰厚度進行了比較，結果顯示T0 高于最大厚度，T1 低于最大厚度，且相關性更強（0.480），值得注意的是T2 各項誤差均低于T0 和T1，相關系數與T1 一致。考慮到過境時間，上午星TERRA過境時間為當地時間10:30 前后，下午星AQUA 過境時間為當地時間13:30 前后，此時海冰厚度比較接近一天之中的平均厚度，因此在與實測數據進行比較時，仍以平均實測冰厚為準。為了進行對比，本文還計算了Yuan 等[26]半經驗模型法與實測數據的誤差，結果顯示與平均實測厚度的各項誤差低于T0 高于T1，但與最大實測厚度最為接近。

最后，選取2020–2021 年冬季渤海石油平臺的觀測數據，與反演結果進行比較。由圖10 可見，總體來看JZ20-2 平臺實測冰厚和反演結果相較于JZ9-3 平臺波動起伏較大，這是由于JZ20-2 平臺位置更加靠近海冰邊緣（圖11），極易受到海冰范圍增加和縮小的影響；改進前的算法結果與平均實測數據相比偏大，有時甚至高于最大實測厚度（如2021 年1 月8 日JZ20-2 平臺），而改進后算法的整體厚度為3 種方法中最小，除1 月5 日、13 日低于JZ20-2 平臺平均實測厚度，1 月8 日高于最大實測厚度以外，多介于平均厚度和最大厚度之間。接下來，通過真彩圖和海冰厚度反演結果（圖11）對以上提到的3 個特殊時間節點進行具體分析：由圖11a 可見，1 月5 日JZ20-2 平臺周邊的海冰比較松散，海冰的漂移會相對更加明顯，平臺觀測時間與衛星觀測時間的不一致可能是導致反演有偏差的原因；另外，密集度較小的海區海冰可能偏薄，而本反演算法對厚度小于6 cm 海冰的反演結果尚不是十分準確。1 月13 日真彩圖顯示JZ20-2 平臺周邊基本沒有海冰（圖11d），海冰邊緣在平臺以北，這與冰厚反演結果一致但與平臺實測不同（圖11h），其原因也與平臺觀測與衛星觀測不同步有關。1 月8 日的JZ20-2 平臺3 種算法反演結果一致大于最大實測冰厚，由圖11c 可見海冰表面明顯偏白，可能有積雪覆蓋，同樣情況的還有1 月7 日真彩圖（圖11b），積雪會明顯提高海冰反照率進而影響反演結果，導致JZ9-3 平臺1 月7 日和8 日冰厚反演結果均顯著高于降雪之前的1 月5 日反演結果。但由于缺少石油平臺降雪量數據支持，通過真彩圖對降雪的判斷是一種猜測。

圖10 2020–2021 年冬季反演結果與實測數據的時間序列比較（JZ20-2 平臺和JZ9-3 平臺）Fig.10 Time series comparison between measured data and retrieval results in 2020–2021 winter (JZ9-3 and JZ20-2)

5 結論與討論

本文利用MODIS 遙感數據實現了渤海海冰范圍的自動化提取，并由渤海的實際情況出發，針對海冰厚度指數關系模型中海水反照率 αsea與 衰減系數 μα參數的設置進行了反演算法的改進，并將改進前后算法結果及前人算法結果與實測數據進行檢驗和比對，主要結論如下：

（1）基于Canny 邊緣檢測算子，加入真彩圖灰度化、高斯模糊、空洞填充、灰度圖二值化和ST 自動化閾值判別法等步驟，提出了一種渤海海冰的自動化提取算法。該算法在一定程度上克服了Canny 邊緣檢測算子對平整薄冰區檢測能力不足的弊端，可以有效提取高濃度泥沙海區的薄冰，并且通過應用最大類間法和粗糙度法實現了冰水分離過程中的閾值自動化提取。

（2）基于渤海的水文特征，對反照率與海冰厚度指數關系模型中的參數重新設置，實現對算法的改進。選取了海冰外緣線外側的狹長條帶內海水像元作為插值節點，向海冰像元點所在位置進行插值，獲得冰覆蓋區海水反照率 αsea，進而將海冰厚度指數關系模型中的海水反照率 αsea由固定值變為隨海域實際情況而改變的動態數值。使用海上石油平臺的冰厚實測數據和MODIS 數據重新反推衰減系數 μα，為了減小偏差，反推過程中排除了實測厚度6 cm 以下的薄冰，確定其估計值為1.74。

（3）利用海上石油平臺實測數據中的測試數據集進行驗證，結果顯示：改進后的海冰厚度指數關系模型的反演結果與實測數據間誤差明顯縮小，平均絕對誤差由7.05 cm 縮小到2.74 cm，相關系數由0.434 提高到0.485；但反演結果較最大實測厚度偏薄，多介于平均厚度與最大厚度之間。

需要說明的是，計算冰下海水反照率 αsea時，由于本文使用的插值算法為反距離加權插值法，理論上最近的插值節點與被插值點之間距離太遠會增大插值結果與真實值之間的差異，最終冰厚反演結果也會受到影響，尤其是重冰期遠離海水的海冰，插值誤差造成的影響更加明顯，因此還需要進一步考慮如何提高αsea取值的合理性。本文利用海冰平均厚度實測數據對衰減系數 μα進行了重新反推計算，所得衰減系數μα為1.74，本文也嘗試采用了最大實測厚度進行反推，所得衰減系數 μα為1.18，代入算法中使得反演結果高于前者47.5%，但由于篇幅所限，未在正文中展開討論；計算時排除了6 cm 以下的薄冰，因此薄冰反演結果的準確性有所下降，并且最終結果不是通過理論推導所得，而是經驗性的估計值，實際上海冰的衰減系數并不是固定不變的，而是上下波動的，如何通過衛星判斷不同的海冰類型，進而獲得更加準確的衰減系數 μα值，是未來繼續優化反演算法的關鍵之一。

海冰表面覆蓋積雪會顯著提高海冰反照率，導致冰厚反演結果偏大，而且積雪厚度越大，對結果造成的影響越大。在未來的研究中，需要先區分冰上積雪和厚冰，再進一步對冰厚反演算法進行修正。同時也建議平臺及岸站觀測部門在監測海冰實況時，增加對降雪情況的記錄。