基于MODIS 數(shù)據(jù)的渤海海冰厚度反演算法優(yōu)化

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

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

1 引言

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

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

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

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

2 數(shù)據(jù)來源

2.1 MODIS 遙感數(shù)據(jù)

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

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

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

2.2 渤海海上石油平臺實測數(shù)據(jù)

本文采用的實測數(shù)據(jù)為渤海海上石油平臺的觀測數(shù)據(jù)，包括2013–2014 年冬季和2015–2016 年冬季的JZ20-2 平臺（40.500°N,121.352°E），2020–2021 年冬季的JZ20-2 和JZ9-3 平臺（40.664°N,121.462°E）冰厚觀測數(shù)據(jù)，觀測時間為當(dāng)?shù)貢r間早上8 點至下午7 點，每小時進(jìn)行一次海冰最大厚度和平均厚度觀測；2009–2010 年冬季和2012–2013 年冬季數(shù)據(jù)引自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 平臺的海冰厚度現(xiàn)場觀測數(shù)據(jù)，測量方法為目視觀測，測量時間為當(dāng)?shù)貢r間下午6 點。平臺具體位置如圖1 所示。本文在使用實測數(shù)據(jù)時，將全部實測數(shù)據(jù)隨機(jī)分為訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集，分別占總數(shù)據(jù)量的63%和37%，訓(xùn)練數(shù)據(jù)集用于海冰厚度算法的改進(jìn)，測試數(shù)據(jù)集用于驗證和分析算法改進(jìn)效果。

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

3 反演算法及改進(jìn)

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

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

3.1 云剔除

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

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

3.2 冰水判別方法的改進(jìn)

3.2.1 海冰范圍提取原理

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

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

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

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

3.2.2 提取裂縫密集區(qū)

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

3.2.3 提取裂縫稀疏薄冰區(qū)

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

3.2.4 誤判像元去除

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

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

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

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

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

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

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

3.3 海冰厚度算法及改進(jìn)

Grenfell[17]和 Allison 等[42]在極地調(diào)查中發(fā)現(xiàn)，當(dāng)海冰厚度從2 cm 增加到9 cm 時，海冰反照率由0.11 增加到0.24。Grenfell[17]提出反照率與海冰厚度呈指數(shù)關(guān)系模型為

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

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

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

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

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

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

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

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

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

3.4 冰厚計算敏感性試驗

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

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

4 反演結(jié)果的對比與驗證

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

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

由于與Yuan 等[26]半經(jīng)驗?zāi)Ｐ头ǖ南嗨菩愿咧淮砼c同類型反演算法的比較情況，無法作為誤差的檢驗標(biāo)準(zhǔn)。因此本文利用渤海海上石油平臺實測數(shù)據(jù)來檢驗算法改進(jìn)前后的結(jié)果；同時為了檢驗本文對反演算法中海水反照率 αsea與衰減系數(shù) μα的改進(jìn)效果，分別改變其中一個參數(shù)而另一個參數(shù)保持不變，將所得結(jié)果與實測數(shù)據(jù)進(jìn)行比較，4 個試驗的結(jié)果見表2。需要指出的是在進(jìn)行誤差分析時使用的數(shù)據(jù)集為測試數(shù)據(jù)集，與3.3 節(jié)反推衰減系數(shù) μα?xí)r使用的訓(xùn)練數(shù)據(jù)集不同。

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

續(xù)表 2

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

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

圖10 2020–2021 年冬季反演結(jié)果與實測數(shù)據(jù)的時間序列比較（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 結(jié)論與討論

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

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

（2）基于渤海的水文特征，對反照率與海冰厚度指數(shù)關(guān)系模型中的參數(shù)重新設(shè)置，實現(xiàn)對算法的改進(jìn)。選取了海冰外緣線外側(cè)的狹長條帶內(nèi)海水像元作為插值節(jié)點，向海冰像元點所在位置進(jìn)行插值，獲得冰覆蓋區(qū)海水反照率 αsea，進(jìn)而將海冰厚度指數(shù)關(guān)系模型中的海水反照率 αsea由固定值變?yōu)殡S海域?qū)嶋H情況而改變的動態(tài)數(shù)值。使用海上石油平臺的冰厚實測數(shù)據(jù)和MODIS 數(shù)據(jù)重新反推衰減系數(shù) μα，為了減小偏差，反推過程中排除了實測厚度6 cm 以下的薄冰，確定其估計值為1.74。

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

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

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