基于遙感技術(shù)的1988—2018年渤海海冰冰情時空變化特征*

孫靜琪，陳笑娟，李 倩，高吉喆，張 化，許映軍

(1.環(huán)境演變與自然災(zāi)害教育部重點實驗室，北京 100875； 2.北京師范大學(xué) 地理科學(xué)學(xué)部，北京 100875; 3.河北省氣象災(zāi)害防御中心，河北 石家莊 050021； 4.應(yīng)急管理部-教育部減災(zāi)與應(yīng)急管理研究院，北京 100875；5.北京師范大學(xué) 國家安全與應(yīng)急管理學(xué)院 ，北京 100875； 6.地表過程與資源生態(tài)國家重點實驗室，北京 100875)

渤海是北半球緯度最低的結(jié)冰海域[1]，渤海及環(huán)渤海地區(qū)是中國的經(jīng)濟中心之一，已發(fā)現(xiàn)勝利油田、遼河油田等豐富的油氣資源，且海上運輸、漁業(yè)養(yǎng)殖等人類經(jīng)濟活動頻繁[2-3]。海冰凍結(jié)和漂移均會對渤海海運、油氣勘探和生產(chǎn)等造成不同程度的影響，形成海冰災(zāi)害。截止目前，海冰已導(dǎo)致了多起石油平臺坍塌、船舶損壞、航道堵塞等災(zāi)害事件發(fā)生，造成了嚴(yán)重的經(jīng)濟損失[4-6]。海冰是全球冰雪圈層的重要組成部分，對氣候變化較為敏感，研究海冰年際變化、日變化差異，可以精確指示局部區(qū)域內(nèi)的氣候變化[7]。海冰淡化因海冰平均鹽度遠(yuǎn)低于海水，有利于降低淡化成本，可以為中國北方地區(qū)解決缺水問題提供新的戰(zhàn)略水資源[8-10]。海冰資源開發(fā)利用的基礎(chǔ)之一就是要較準(zhǔn)確地估算海冰儲量，而海冰面積和厚度時空變化特征研究是海冰資源量估算的基礎(chǔ)[11]。自20世紀(jì)70年代以來，遙感在海冰監(jiān)測中發(fā)揮著重要作用[12-15]。基于主被動遙感數(shù)據(jù)可長期、大范圍、準(zhǔn)實時、動態(tài)地監(jiān)測渤海海冰冰情參數(shù)變化[16-17]。因此，基于多源遙感數(shù)據(jù)研究海冰冰情時空變化特征，對預(yù)防和減輕海冰災(zāi)害、研究環(huán)渤海沿岸地區(qū)氣候變化、估算海冰資源儲量和完善中長期海冰預(yù)報統(tǒng)計模型等均有重要意義。但是，遙感衛(wèi)星由于時空分辨率低、云層遮擋等原因，難獲取海冰冰情指標(biāo)長時序連續(xù)日數(shù)據(jù)，因此，國內(nèi)外的學(xué)者們利用光學(xué)遙感例如MODIS和AVHRR提取的海冰數(shù)據(jù)來建立渤海海冰和氣象因子的統(tǒng)計學(xué)模型。目前，很少有學(xué)者利用GOCI數(shù)據(jù)來研究長時序的渤海海冰，而GOCI的高時間分辨率和高空間分辨率更適合對渤海海冰進(jìn)行大范圍、近實時監(jiān)測。此外，有關(guān)渤海海冰長達(dá)數(shù)十年的面積和厚度的日數(shù)據(jù)集的研究也很少。本文基于AVHRR、MODIS、GOCI三種遙感數(shù)據(jù)解譯出的1988—2018年的渤海海冰面積和厚度數(shù)據(jù)集分析了渤海海冰冰情指標(biāo)的變化特征，以及結(jié)合氣象資料建立統(tǒng)計學(xué)模型，補全了1988—2018年渤海海冰面積日數(shù)據(jù)集。

1 研究區(qū)概況

渤海是中國唯一的內(nèi)海，三面被陸地環(huán)繞，也是中緯度典型季風(fēng)氣候區(qū)，位于37°～41°N，117°～123°E之間，每年冬季都會結(jié)冰，是北半球緯度最低的有冰海域[18]。渤海東北西南向長約300 n mile，東西向長約187 n mile，面積約為77 284 km2，平均水深18 m，最大水深85 m。海冰的形成和分布是區(qū)域氣象要素和物理海洋要素共同作用的結(jié)果，具有明顯的區(qū)域差異。根據(jù)地理位置不同及冰情差異，渤海可分成遼東灣、渤海灣、萊州灣三個海區(qū)[19]。每年冬季渤海冰冰期長短、分布范圍大小、厚度、表面堆積程度等因素差異較大，年內(nèi)和年際變化明顯，海冰災(zāi)害時有發(fā)生，造成了一定的區(qū)域社會經(jīng)濟等損失。例如2009/2010年冬季，渤海及黃海北部發(fā)生了海冰災(zāi)害，造成遼寧、河北、天津、山東等沿海三省一市受災(zāi)人口6.1萬人，船只損毀7 157艘，港口及碼頭封凍296個，水產(chǎn)養(yǎng)殖受損面積約2.08×105hm2，因災(zāi)直接經(jīng)濟損失63.18億元(災(zāi)情數(shù)據(jù)來源于2010年中國海洋災(zāi)害公報)。因此，渤海海冰冰情監(jiān)測及其時空變化特征研究對于海冰災(zāi)害防災(zāi)減災(zāi)具有重要意義。

圖1 研究區(qū)地理位置(該圖根據(jù)審圖號為GS(2015)732號的標(biāo)準(zhǔn)地圖制作，邊界無修改)

2 數(shù)據(jù)及方法

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

遙感數(shù)據(jù)來源于美國國家海洋和大氣管理局(NOAA)系列氣象衛(wèi)星上搭載的多光譜通道高分辨率輻射掃描儀(Advanced Very High Resolution Radiometer, AVHRR)1988年1月至2000年3月期間累計298幅，美國宇航局大型空間遙感儀器中分辨率成像光譜儀(Moderate-resolution Imaging Spectroradiometer, MODIS)2000年12月至2015年3月期間累計385幅，以及韓國新一代海洋水色遙感器所涉地球靜止海洋彩色成像儀，(Geostationary Ocean Color Imager, GOCI)2015年12月至2018年3月期間累計184幅無云有冰、質(zhì)量較好的影像。

氣象數(shù)據(jù)來源于中國地面國際交換站氣候日數(shù)據(jù)集，可從中國氣象數(shù)據(jù)服務(wù)中心網(wǎng)站(http://data.cma.cn/en)下載，包括渤海沿岸18個氣象站點(營口、熊岳、大連、綏中、瓦房店、興城、黃驊、秦皇島、樂亭、唐山、天津、塘沽、東營、龍口、惠民、濰坊、長島、大洼)1988年1月至2018年3月期間平均氣溫、最低氣溫和最高氣溫的觀測數(shù)據(jù)。

2.2 方法

在海冰面積提取和海冰厚度反演之前，需要先對下載下來的遙感圖像做預(yù)處理，其中包含格式轉(zhuǎn)換、輻射校正、幾何校正和大氣校正等。

2.2.1 海冰面積提取

對于AVHRR圖像，基于分區(qū)閾值法將渤海劃分成遼東灣、渤海灣、萊州灣、渤海中部四個海洋環(huán)境相對較均一的區(qū)域，使各子區(qū)域內(nèi)海水反射率小于海冰反射率，根據(jù)冰-水反射率分區(qū)設(shè)定閾值提取海冰面積信息[20-21]。對于MODIS數(shù)據(jù)和GOCI數(shù)據(jù)采用基于樣本的面向?qū)ο蠓诸惙椒ㄌ崛『１娣e，分割尺度與海冰最小尺寸相關(guān)，設(shè)置為20；合并尺度則表示類相似度，設(shè)置為70。上述步驟之后，通過目視解譯進(jìn)一步校正，去除明顯錯分的海冰多邊形，最終得到海冰覆蓋區(qū)域。

2.2.2 海冰厚度估算

海冰面積提取后，采用太陽短波輻射吸收理論的海冰遙感反演模型來反演厚度，并確定不同分區(qū)遙感反演模型相關(guān)參數(shù)，在此基礎(chǔ)上得到渤海海冰厚度分布。海冰光學(xué)性質(zhì)的理論模式和實驗觀測表明一定的厚度范圍內(nèi)，海冰反照率隨冰厚的增加而增加[22]，海冰厚度與海冰反照率之間呈現(xiàn)良好指數(shù)關(guān)系[23]，尤其對于一年冰或冰厚小于1 m的海冰而言，海冰反照率反演模型具有時間和空間上的優(yōu)勢[24-25]。海冰厚度與反照率之間的關(guān)系式可表達(dá)如下：

α(h)=αmax[1-k·exp(-μαh)]。

(1)

式中：α(h)是與海冰厚度相關(guān)的海冰反照率；αmax為無限大冰厚斜對應(yīng)的反照率；k為與αmax有關(guān)的相關(guān)系數(shù)；μα是關(guān)于反照率的衰減系數(shù)，h是冰厚(m)。

2.2.3 傳感器校正

為了更準(zhǔn)確的獲取海冰面積和厚度時間序列，本文選擇交叉定標(biāo)進(jìn)行傳感器間校正[26]。利用2000—2011年間的MODIS和AVHRR的數(shù)據(jù)集建立線性回歸模型，然后將獲得的模型用于校準(zhǔn)AVHRR結(jié)果，因為MODIS結(jié)果似乎更準(zhǔn)確[27]。選擇2012—2015年間MODIS和GOCI數(shù)據(jù)集建立線性回歸模型，并對GOCI提取海冰面積進(jìn)行傳感器間校正，選擇MODIS數(shù)據(jù)作為標(biāo)定基準(zhǔn)進(jìn)而校正GOCI數(shù)據(jù)，最終將不同分辨率的遙感數(shù)據(jù)統(tǒng)一成MODIS數(shù)據(jù)所對應(yīng)的分辨率。這是由于MODIS數(shù)據(jù)與AVHRR和GOCI數(shù)據(jù)序列均有部分時段重疊。通過傳感器間校準(zhǔn)提高海冰面積和厚度時間序列的一致性和準(zhǔn)確性。

2.2.4 氣象數(shù)據(jù)計算

通常將累積凍冰度日(CFDD)定義為一定時間段內(nèi)每日平均氣溫和海洋冰凍溫度之差的總和(當(dāng)平均地面氣溫低于海洋冰凍溫度)，累積融冰度日(CTDD)定義為一定時間段內(nèi)海洋冰凍溫度和每日平均氣溫之差的總和(當(dāng)平均氣溫高于海洋冰凍溫度)[28]，如式(2)和式(3)所述。

(2)

(3)

式中：Tf是渤海的海洋冰凍溫度，參考標(biāo)準(zhǔn)選取-1.8℃[29]；Ta是日平均氣溫，本文取渤海沿海的18個臺站日均氣溫的平均值；t是冰期天數(shù)，ds和de分別為開始和結(jié)束的日期。

通常將低于0℃的日平均氣溫累加值稱為負(fù)積溫[30]，本文將3 d-t℃積溫定義為連續(xù)的3 d內(nèi)低于t℃的日平均氣溫之和。

2.2.5 海冰指標(biāo)計算

國家海洋局每年以遼東灣、渤海灣和萊州灣三個海區(qū)的海冰結(jié)冰范圍和冰厚作為主要指標(biāo)，把海冰冰情劃分成5個等級：輕冰年(1.0)、常冰年偏輕(2.0)、常冰年(3.0)、常冰年偏重(4.0)、重冰年(5.0)，并在這些等級范圍內(nèi)，細(xì)分出常冰年明顯偏輕(1.5)、常冰年略偏輕(2.5)、常冰年略偏重(3.5)和常冰年明顯偏重(4.5)4個亞等級。1989—2018年間渤海冰情等級來自于國家海洋局每年發(fā)布的海洋公報數(shù)據(jù)；1988—1989年的渤海冰情等級數(shù)據(jù)來自于文獻(xiàn)[31]。

本文將海冰曾結(jié)冰頻率定義為：

(4)

式中：Pij為第i月第j個柵格的曾結(jié)冰頻率，ei總為該柵格在1988—2018年i月的遙感影像總數(shù)，ei結(jié)冰為該柵格在1988—2018年i月的遙感影像中結(jié)冰的天數(shù)。

冰期指海冰從出現(xiàn)到消融的總天數(shù)(一般為12月至次年的3月)。本文以海冰面積達(dá)到最大值的日期為分界，將冰期分為凍冰期和融冰期，即凍冰期為出現(xiàn)海冰的第1 d到此分界日的總天數(shù)，融冰期為此分界日到海冰消失至無的總天數(shù)。

2.2.6 多元非線性回歸模型構(gòu)建

選取與海冰面積相關(guān)系數(shù)高的四個氣象因子建立多元非線性回歸模型，利用最小二乘法分別在不同發(fā)展階段擬合出各個階段的系數(shù)，補全遙感影像缺失的日期。

3 結(jié)果

3.1 1988-2018渤海海冰時間序列變化特征

3.1.1 年際變化

1988—2018年渤海冰情等級時間序列如圖2所示。整體來看渤海冰情等級有明顯的年際波動，且2000年之后的年際波動更為明顯。劉煜等人[32]的研究表明黃渤海的冰情等級從1950—1980年平均冰情等級的2.9級下降到1980—2010年平均冰情等級的2.3級，降幅為0.6，與全球和東亞地區(qū)氣候變暖的趨勢一致。本文根據(jù)資料[16,31]記錄的冰情等級，計算了1950—1984年和1985—2018年的平均冰情等級，分別為2.91和2.40，下降了0.51。然而，以冰情等級變化來看，1989—2003年的平均冰情等級比2004—2018年的平均冰情等級上升了0.3，表明2004年后的冰情比之前略有加重且波動明顯。1988—2018年間冰情等級最低為2002年的1.0級，屬于輕冰年；最高為2001年的4.0級，屬于常冰年偏重。1988—2018年間并未出現(xiàn)重冰年，而屬于常冰年偏輕(2.0)的年份最多，共有10個年份。

圖2 1988—2018年渤海冰情等級時間序列圖

圖3中可以看到1988—2018年間渤海年平均海冰面積數(shù)值主要在4 000～20 000 km2之間，年際變化明顯。年均海冰面積均值為8 390.96 km2。其中，最大值出現(xiàn)在2010年，為18 895.59 km2，遠(yuǎn)大于平均值，這一年發(fā)生了1988—2018年間最嚴(yán)重的一次海冰災(zāi)害。最小值出現(xiàn)在2015年，為4 533.01 km2。2015年也是冰情等級最低的年份之一，冰情等級為1.0。圖3還可以看到1988—2018年間海冰厚度年平均數(shù)值主要在5～17 cm之間，有明顯的年際變化，年均海冰厚度均值為11.09 cm。其中，年均海冰厚度最大值也出現(xiàn)在2010年，為16.16 cm；最小值出現(xiàn)在2007年，為5.38 cm。2002年冰情等級為1.0，年內(nèi)平均海冰厚度也非常小，為5.66 cm，接近最小值，也是冰情等級最低的年份之一。

圖3 1988—2018年渤海年平均海冰面積和海冰厚度時間序列

3.1.2 年內(nèi)變化

從圖4中可以看出，1988—2018年渤海日平均海冰面積基本呈現(xiàn)為先凍結(jié)增加達(dá)到年最大值，之后再消融減少，整個過程呈單峰型。日平均海冰面積的峰值出現(xiàn)在1月25日，對應(yīng)日平均海冰面積為16 278.74 km2。渤海日平均海冰厚度變化也基本呈現(xiàn)為單峰型，即日平均海冰厚度先增加達(dá)到峰值，再減小。日平均海冰厚度峰值也出現(xiàn)在1月25日，對應(yīng)日平均海冰厚度為16.37 cm。從平均趨勢來看，海冰凍融過程轉(zhuǎn)折點為1月25日，此前海冰整體呈增加的趨勢，此后呈消融趨勢，逐漸消融至無。

圖4 1988—2018年渤海日平均海冰面積和海冰厚度時間序列

3.2 1988-2018渤海海冰空間分布特征

3.2.1 海冰曾結(jié)冰頻率空間分布

從圖5來看，渤海海冰曾結(jié)冰頻率空間分布的特征是遼東灣曾結(jié)冰的范圍最大，其次是渤海灣，最后是萊州灣。1988—2018年冬季渤海曾結(jié)冰的范圍達(dá)到54 544 km2，占渤海總面積的70.58%，其中遼東灣約占結(jié)冰范圍的63%。冬季3個月中，2月份渤海曾結(jié)冰范圍最大，達(dá)到了49 655 km2，其中遼東灣曾結(jié)冰范圍約占68%。12月份渤海曾出現(xiàn)海冰面積約占渤海總面積的39.88%，1月份渤海約60.90%的范圍以及2月份約64.25%的范圍曾出現(xiàn)海冰。

圖5 1988—2018年渤海曾結(jié)冰頻率空間分布

3.2.2 海冰厚度空間分布

從1988—2018年渤海海冰12月、1月、2月和冬季平均厚度空間分布圖(圖6)來看，渤海海冰冰厚空間分布的特征是離岸越近，厚度越大，遼東灣的冰厚最大，其次是渤海灣，萊州灣的冬季平均海冰厚度是三個海區(qū)中最小的。對于遼東灣海區(qū)來說，東岸的冬季平均海冰厚度比西岸的更大，冬季平均海冰厚度的最大值出現(xiàn)在營口市和蓋州市附近，約為0.24 cm。沿岸的海冰厚度相較于其他區(qū)域更厚，且渤海中間區(qū)一般沒有海冰的分布，海冰厚度平均值為0。從時間上來看，三個海區(qū)1月和2月的平均冰厚明顯比12月的平均厚度增大，且海冰范圍增大、海冰外緣線離岸距離增大。

圖6 1988—2018年渤海海冰平均厚度空間分布

3.3 海冰冰情指標(biāo)對氣象因子的響應(yīng)

海冰的生消與該地區(qū)的溫度密切相關(guān)[33]，本文選取一些相關(guān)的氣象因子例如：日平均溫度、日最低溫度、3 d-1.8℃積溫、累積凍冰度日等來分析渤海海冰冰情指標(biāo)與氣象因子的關(guān)系。本節(jié)使用的凍冰期和融冰期的分界點是圖4得出的1月25日，將從出現(xiàn)海冰到次年1月25日算做凍冰期，每年的1月25日到海冰消融算做融冰期。

對比表1看出在海冰面積和氣象因子的相關(guān)性分析中，以12月到次年3月期間海冰面積最大值出現(xiàn)的日期作為分界點，將海冰凍融過程分為凍冰期和融冰期分別與氣象因子做相關(guān)性分析的相關(guān)程度高于直接將冰期與氣象因子做相關(guān)分析的相關(guān)程度。尤其是凍冰期海冰面積與CFDD的相關(guān)性為高度相關(guān)(r=0.821)，明顯高于總冰期海冰面積與CFDD的相關(guān)性(r=0.524)；融冰期海冰面積與CTDD的相關(guān)性(r=-0.521)也明顯高于總冰期海冰面積與CTDD的相關(guān)性(r=-0.132)。在凍冰期，與海冰面積相關(guān)系數(shù)最高的氣象因子是CFDD；在融冰期，與海冰面積相關(guān)系數(shù)最高的氣象因子是3 d-1.8℃積溫。

表1 海冰面積厚度與氣象因子相關(guān)系數(shù)

在海冰厚度和氣象因子的相關(guān)性分析中，將海冰凍融過程分為凍冰期和融冰期分別與氣象因子做相關(guān)性分析的相關(guān)程度也高于直接將冰期與氣象因子做相關(guān)分析的相關(guān)程度。尤其是凍冰期海冰厚度與CFDD的相關(guān)系數(shù)為0.465，明顯高于總冰期海冰面積與CFDD的相關(guān)系數(shù)0.257。在凍冰期，與海冰厚度相關(guān)系數(shù)最高的氣象因子是CFDD；在融冰期，與海冰厚度相關(guān)系數(shù)最高的氣象因子是CTDD。

3.4 海冰面積時間序列數(shù)據(jù)集補全

多元非線性回歸方程：

(5)

式中：Tmin為日最低氣溫，T為3 d-1.8℃積溫，CFDD為累積凍冰度日，CTDD為累積融冰度日。凍冰期系數(shù)：a1=-386.164,a2=-132.346,a3=-3.254,a4=135.72,a5=-70.41,a6=2.932,a7=1.072,a8=-0.685,a9=-3.824,a10=0.124,a11=-0.003,a12=-0.004,a13=-683.058。融冰期系數(shù)：a1=-207.436,a2=-114.502,a3=-58.967,a4=99.492,a5=109.589,a6=13.155,a7=-0.685a8=-1.387,a9=5.137,a10=0.32,a11=-0.002,a12=0.005,a13=1 541.771。

凍冰期的多元非線性回歸方程的R2=0.803；融冰期的多元非線性回歸方程的R2=0.689。

從圖7中可以看出，2016年2月2日之后的幾天和2018年2月20日前后的幾天，真實值比預(yù)測值偏大的較為明顯。2016年和2018年的冰情等級分別是3.0和2.5，比用來回歸擬合的26個年份的平均冰情等級2.4大，屬于偏重的冰情年份，這或許導(dǎo)致了預(yù)測值在此時間段出現(xiàn)低估的現(xiàn)象。2015年1月21日之后的幾天，真實值比預(yù)測值偏小的較為明顯，而2015年的冰情等級為1.0，屬于輕冰年，這或許導(dǎo)致了預(yù)測值在此時間段出現(xiàn)高估的現(xiàn)象。由于輕冰年和重冰年的海冰面積差異較大，此回歸模型在預(yù)測重冰年的海冰面積時，可能會一定程度上的低估，在預(yù)測輕冰年的時候可能出現(xiàn)一定程度上的高估。如果遙感提取樣本數(shù)足夠多，可以按輕冰年、常冰年和重冰年分別建立回歸模型，這會使得預(yù)測的結(jié)果更接近真實值。目前，遙感提取海冰面積與回歸模型計算的海冰面積整體趨勢一致，且回歸模型具有較好的擬合優(yōu)度，具有一定的補全和預(yù)報價值。

圖7 2014年1月—2018年3月海冰面積遙感提取值與回歸計算值對比

4 結(jié)論與討論

(1)1988—2018年中，2004年后冰情比之前略有加重且波動明顯；年平均海冰面積和厚度最大值均出現(xiàn)在2010年，年平均海冰面積最小值出現(xiàn)在2015年，而年平均海冰厚度最小值出現(xiàn)在2007年。

(2)渤海海冰年內(nèi)變化特征為日平均海冰面積和海冰厚度呈現(xiàn)為先增加達(dá)到最大值之后再減少，整個過程呈單峰型。日平均海冰面積和日平均海冰厚度的峰值均出現(xiàn)在1月25日。

(3)渤海海冰主要分布在遼東灣，1988—2018年間曾出現(xiàn)過海冰的范圍占渤海面積的70.58%，達(dá)到54 544 km2。渤海海冰冰厚空間分布的特征是離岸越近，厚度越大，渤海中間區(qū)域不曾結(jié)冰。且遼東灣的冬季平均海冰厚度最大，其次是渤海灣，最小的是萊州灣。

(4)在凍冰期，與海冰面積和海冰厚度相關(guān)系數(shù)最高的氣象因子均是CFDD；在融冰期，與海冰面積相關(guān)系數(shù)最高的氣象因子是3 d-1.8℃積溫，與海冰厚度相關(guān)系數(shù)最高的氣象因子是CTDD。

(5)凍冰期的多元非線性回歸方程的R2=0.803，融冰期的R2=0.689。將回歸計算的海冰面積與遙感提取的海冰面積做對比，其整體的趨勢一致，有一定的補全和預(yù)報價值。

本文建立的多元非線性回歸模型在逐日海冰擬合上有一定的預(yù)測能力，但僅能反映整個海區(qū)的平均值，沒有考慮渤海不同海域的差異，沒有空間分布的預(yù)測；由于樣本數(shù)量的原因，沒有將其按冰情等級分別建立模型，導(dǎo)致該模型在較重冰年會出現(xiàn)低估而在較輕冰年會出現(xiàn)高估；且僅用了氣象相關(guān)的因子，沒有考慮其他影響因素。