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

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

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

1 研究區(qū)概況

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

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

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ǎng)站(http://data.cma.cn/en)下載，包括渤海沿岸18個氣象站點(營口、熊岳、大連、綏中、瓦房店、興城、黃驊、秦皇島、樂亭、唐山、天津、塘沽、東營、龍口、惠民、濰坊、長島、大洼)1988年1月至2018年3月期間平均氣溫、最低氣溫和最高氣溫的觀測數(shù)據(jù)。

2.2 方法

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

2.2.1 海冰面積提取

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

2.2.2 海冰厚度估算

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

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

(1)

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

2.2.3 傳感器校正

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

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

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

(2)

(3)

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

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

2.2.5 海冰指標計算

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

本文將海冰曾結冰頻率定義為：

(4)

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

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

2.2.6 多元非線性回歸模型構建

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

3 結果

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，遠大于平均值，這一年發(fā)生了1988—2018年間最嚴重的一次海冰災害。最小值出現(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)為先凍結增加達到年最大值，之后再消融減少，整個過程呈單峰型。日平均海冰面積的峰值出現(xiàn)在1月25日，對應日平均海冰面積為16 278.74 km2。渤海日平均海冰厚度變化也基本呈現(xiàn)為單峰型，即日平均海冰厚度先增加達到峰值，再減小。日平均海冰厚度峰值也出現(xiàn)在1月25日，對應日平均海冰厚度為16.37 cm。從平均趨勢來看，海冰凍融過程轉(zhuǎn)折點為1月25日，此前海冰整體呈增加的趨勢，此后呈消融趨勢，逐漸消融至無。

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

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

3.2.1 海冰曾結冰頻率空間分布

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

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

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 海冰冰情指標對氣象因子的響應

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

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

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

在海冰厚度和氣象因子的相關性分析中，將海冰凍融過程分為凍冰期和融冰期分別與氣象因子做相關性分析的相關程度也高于直接將冰期與氣象因子做相關分析的相關程度。尤其是凍冰期海冰厚度與CFDD的相關系數(shù)為0.465，明顯高于總冰期海冰面積與CFDD的相關系數(shù)0.257。在凍冰期，與海冰厚度相關系數(shù)最高的氣象因子是CFDD；在融冰期，與海冰厚度相關系數(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日前后的幾天，真實值比預測值偏大的較為明顯。2016年和2018年的冰情等級分別是3.0和2.5，比用來回歸擬合的26個年份的平均冰情等級2.4大，屬于偏重的冰情年份，這或許導致了預測值在此時間段出現(xiàn)低估的現(xiàn)象。2015年1月21日之后的幾天，真實值比預測值偏小的較為明顯，而2015年的冰情等級為1.0，屬于輕冰年，這或許導致了預測值在此時間段出現(xiàn)高估的現(xiàn)象。由于輕冰年和重冰年的海冰面積差異較大，此回歸模型在預測重冰年的海冰面積時，可能會一定程度上的低估，在預測輕冰年的時候可能出現(xiàn)一定程度上的高估。如果遙感提取樣本數(shù)足夠多，可以按輕冰年、常冰年和重冰年分別建立回歸模型，這會使得預測的結果更接近真實值。目前，遙感提取海冰面積與回歸模型計算的海冰面積整體趨勢一致，且回歸模型具有較好的擬合優(yōu)度，具有一定的補全和預報價值。

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

4 結論與討論

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

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

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

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

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

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