基于GOCI的渤海海域赤潮信息遙感監測與分析

程 玉，張圣佳，李金宇，朱紅春

(山東科技大學 測繪與空間信息學院，山東 青島 266590)

赤潮又名紅潮，是由海水中的某些浮游植物、原生動物或細菌在適宜的海域條件下，短時間突發性的劇增或聚集引起水體變色的異?，F象[1]。赤潮不僅對漁業資源尤其是海水養殖業造成影響，還會給沿海旅游業造成經濟損失。此外，赤潮還對海洋生物造成巨大威脅，人類可能因意外飲食受赤潮污染的海洋產品造成中毒或死亡[2-3]。2009—2019年間，渤海共發生赤潮92次，主要發生在每年的5—9月份，其中，6—8月份赤潮發生的頻率最高。赤潮已經成為我國三大近海污染之一，防止和減少赤潮造成的危害已成為十分緊迫的任務。

赤潮的形成機理十分復雜，雖然赤潮發生的具體原因還未詳細查明，但普遍認為海域的富營養化是赤潮發生的物質基礎，溫度、鹽度、氮磷含量是發生赤潮的主要條件，氣溫、降水、氣壓等是赤潮發生的誘導因素[4]。目前，尚無有效的手段能阻止赤潮的發生，只能通過監測和預報赤潮盡可能減少赤潮的危害。

傳統的赤潮監測主要依靠船只調查和岸濱人工觀測。常規船只監測赤潮需要對采樣進行生物、化學分析。飛機等飛行物監測的費用高且受天氣限制，數據覆蓋面小，同步性差。衛星遙感技術為赤潮的快速、大范圍連續監測提供了重要手段[5]。常用的赤潮監測方法有葉綠素濃度法、波段比值法、赤潮指數法等[6-9]。姜德娟等[10]基于中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer,MODIS)數據，比較分析了葉綠素濃度異常法、波段比值法、紅波段差值法、karenia brevis指數法等對渤海赤潮遙感監測的有效性，結果表明幾種方法都能較好地監測到赤潮，葉綠素濃度異常法的赤潮監測效果最好。薛存金等[11]針對單一海洋參數對赤潮提取不足的問題，提出一種多海洋參數赤潮監測方法，利用海表面溫度、葉綠素a濃度、懸浮泥沙濃度等海洋參數，根據專家知識設置參數閾值，設計出一套渤海赤潮災害提取規則，成功提取渤海灣的赤潮。李陽東等[12]利用懸浮物濃度和赤潮指數，通過海洋災害公報面積設置不同的閾值組合提取赤潮，以2011年5月29日溫州海域和2014年5月23日舟山海域的赤潮事件為例進行赤潮日變化規律研究。利用海洋參數提取赤潮，需要設置合適的閾值進行圖像分割，從而達到赤潮提取的目的，因此，分割閾值的確定是一個重點問題。但是，目前所用的閾值分割方法大部分基于經驗統計閾值，或根據海洋公報面積，或假彩色圖像通過目視匹配的方法設置閾值，這些方法不僅繁瑣而且還具有較強的主觀性，因此，需要采取自適應閾值分割方法確定最優閾值，減少人工干預。此外，對于赤潮的研究大都集中在赤潮的提取和變化監測方面，對于赤潮演變導致水色參數變化的研究很少。

針對這一問題，本研究探討了赤潮發生前后海洋水色參數的變化特征。通過對葉綠素a和赤潮指數圖像進行多種方法的自適應閾值分割提取，根據效果從中選取合適的赤潮提取方法，分析赤潮的晝變化特征，并對赤潮發展演變造成的影響進行研究。

1 研究區及數據

1.1 研究區概況

渤海位于37°07′～41°00′N，117°35′～122°16′E，是半封閉內海，僅東面與黃海相連，也是中國最大的內海，面積約77 000 km2，海岸線全長約3 800 km，平均水深18 m。渤海處于北溫帶，受北方大陸性氣候影響，夏無酷暑，冬無嚴寒，全年區域溫差不顯著。沿岸周邊河流眾多，其中主要有黃河、海河、遼河和灤河等。河流攜帶周邊三省一市的工農業廢水和生活污水直接流入渤海，陸域活動的影響直接導致渤海水質下降和富營養化，呈現二類水體特征，赤潮災害頻繁發生。

1.2 實驗數據

地球靜止海洋水色成像儀(geostationary ocean color imager, GOCI)搭載于韓國通信、海洋和氣象衛星(communication ocean and meteorological satellite, COMS)，空間分辨率為500 m(圖1)；時間分辨率極高，可以做到每天8次，每次間隔1小時的觀測頻率；波段范圍412～865 nm，包括6個可見光波段和2個近紅外波段。在本研究中，根據海洋公報中公布的赤潮災害信息，分別下載了2014年4月28日、5月5日、5月8日、5月12日、5月15日、5月26日、5月28日、6月9日、6月12日、6月13日和6月23日的L1B數據，用于赤潮晝變化和海洋參數變化分析。幾何校正采用自帶的經緯度信息生成地理位置查找表文件(geographic lookup
Table, GLT)，利用GOCI數據處理軟件(GOCI data progressing software, GDPS)內置的二類水體大氣校正算法實現大氣校正，獲得反射率數據，后續進行懸浮物、有色溶解有機物(chromophoric dissolved organic matter, CDOM)等海洋水色參數反演。本研究中葉綠素a濃度、懸浮物、漫衰減系數(Kd(490))和CDOM等水色參數均基于GOCI影像反演計算得到，歸一化熒光高度和海表面溫度數據來自美國國家航空航天局海洋水色網站[13]。

圖1 渤海區域(綠色邊界)GOCI真彩色合成影像圖

2 研究方法

2.1 赤潮指數

Ahn等[14]基于寬視場水色掃描儀(sea-viewing wide field-of-view sensor, SeaWiFS)數據，提出利用赤潮指數(red tide index, RI)來表征光學特性復雜的東北亞海域爆發赤潮的可能性。宋德彬[15]基于GOCI數據分析了赤潮和潔凈水體的反射率曲線，其中赤潮水體在550 nm處的反射率明顯較高；與潔凈水體相比較，在443和490 nm處，赤潮水體反射率與550 nm處的差異顯著增大?；诖瞬町?，采用張賀等[16]改進的赤潮指數來進行渤海赤潮信息的提取，表達式如下：

(1)

其中：RI表示赤潮指數；Rrs(443)、Rrs(490)和Rrs(555)分別表示443、490和555 nm波段的遙感反射率。使用藍綠波段反射率來構建赤潮指數，并通過減去443 nm波段的反射率來減少懸浮泥沙對赤潮信息提取的影響。

2.2 閾值確定方法

2.2.1 最大類間方差法

最大類間方差法(OTSU)要求選擇最佳閾值，使圖像前景和背景間的灰度值差距最大，衡量灰度差異的標準定義為類間方差，即當前景和背景間的類間方差取最大值時，獲得最佳閾值。假設圖像的像素總數為N，L個灰度級[1,2,…,L]，灰度值為i的像素總數為ni，像素i的出現概率

(2)

設圖像的初始閾值為t，若像素的灰度值小于t，則屬于C0類，否則屬于C1類，C0和C1分別代表圖像的背景類和前景類。圖像背景類和前景類像素出現的概率和灰度均值分別用ω0、ω1、μ0、μ1表示，有：

(3)

(4)

(5)

其中，σ2表示類間方差，當σ2取最大值時，獲得最佳閾值。

2.2.2 K-means聚類

K-means聚類算法依據數據中點的歐式距離將給定的數據歸入需要劃分的K類中。假設將數據集X(X中的數據個數為n)分為K個類別，將歐式距離作為相似度指標，起始階段在數據集中隨機選取k個值分別作為每一類的中心，之后計算數據集中每一個點x與各個中心u之間的歐式距離，并且依據歐式距離的大小將每個點歸入與之最近的中心所屬類別中，使各類別的聚類平方和I最小，即

(6)

之后更新每個聚類的中心位置，新的中心位置即為每一類所有點的平均值，不斷重復這一過程，當聚類平方和I達到最小時重復過程終止，此時的分類為K-means聚類的最終結果。

2.2.3 最大熵法

最大熵法通過引入信息論中最大熵原則，統計每個灰度值的出現概率，計算灰度級的熵，保證各類灰度的總熵最大，獲得使前景或背景內部灰度盡可能均勻的最優閾值，或是尋找使分割前后圖像的信息量差異最小的最優閾值。假設某個灰度級出現概率為pi，則灰度熵的計算公式為：

(7)

設灰度值q為初始閾值，該閾值分割的背景類C0和前景類C1的累計概率用P0(q)和P1(q)表示，

(8)

背景類和前景類的概率密度函數可分別表示為：

(9)

(10)

背景類和前景類對應的熵分別表示如下：

(11)

(12)

在該閾值下，圖像總熵為：

H(q)=H0(q)+H1(q)。

(13)

計算所有分割閾值下的圖像總熵，將最大熵對應的分割閾值作為最佳分割閾值。

3 赤潮信息提取及晝變化分析

通過比較渤海區域云量，選取2014年5月28日云量最少的1景GOCI影像，即編號為00(8:00)的影像，分別使用OTSU、K-means聚類以及最大熵方法對葉綠素a和RI指數計算結果確定最佳的赤潮分割閾值，利用該閾值進行遙感影像二分類，得到赤潮類型的提取結果并統計其面積。以目視解譯結果為參考結果，采用不同方法進行閾值分割后提取的面積及誤差見表1。

表1 3種閾值分割方法面積提取及誤差分析

3種閾值分割方法對葉綠素a影像確定的閾值不斷增大，OTSU方法確定的閾值為10.5 mg·m-3，分割后得到的面積最接近目視解譯面積，為8 089.25 km2，相對誤差為9.1%；K-means聚類方法確定的閾值為11.7 mg·m-3，赤潮面積為6 619.5 km2，相對誤差25.6%；最大熵法確定的閾值為12.3 mg·m-3，提取的赤潮面積為6 100 km2，相對誤差為31.4%，是3種閾值分割算法中誤差最大的算法。對于RI圖像閾值分割結果，OTSU算法的赤潮提取結果最好，相對誤差為5.8%；K-means聚類算法提取的赤潮面積相對誤差為8.7%；最大熵法的相對誤差結果最大，為9.4%。從3種方法提取面積的誤差分析，OTSU算法對于葉綠素a影像和RI影像中赤潮提取結果最好，K-means聚類方法次之，最大熵法的結果最差。對于葉綠素a影像，3種分割方法提取赤潮面積的平均相對誤差為22%，而RI影像3種方法的平均相對誤差為7.9%，因此，RI影像更適合赤潮提取與分析。綜合來看，使用RI影像并采用OTSU閾值分割算法提取赤潮更有效，其赤潮提取結果見圖2。根據海洋公報的記錄，2014年5月28日， 河北秦皇島灤河口至遼寧綏中海域發生大規模赤潮，實驗中提取的赤潮范圍與赤潮記錄基本一致，可以證明RI影像采用OTSU閾值分割算法提取赤潮準確率較高。赤潮的發生是多因素共同作用的結果，通過查詢歷史氣象資料，2014年5月24日，河北秦皇島附近地區有中雨，最高氣溫18 ℃，降雨導致河流流量增加，攜帶更多營養鹽進入海洋，在5月25日之后，氣溫上升5～8 ℃，最高氣溫達30 ℃，且天氣狀況良好，充足的光照有利于浮游植物的生長，促進赤潮的大范圍暴發。本次赤潮的主要藻類是抑食金球藻，其最適宜的生長溫度為20 ℃[17]，秦皇島附近海域5月26日的海表面溫度在19 ℃左右，基本滿足抑食金球藻的溫度條件，促進了藻類的快速增殖。而未發生赤潮的遼寧錦州附近海域雖然也有降雨，但氣溫上升幅度較小，海表面溫度在14 ℃左右，不利于藻類的生長，因此錦州附近海域并未有赤潮暴發，同時也間接證明提取的赤潮結果準確性較高。

圖2 赤潮提取結果和兩個感興趣區域(黃色方框1和2)

赤潮指數圖像能反映赤潮在一天內的變化，通過閾值分割提取赤潮面積，可大致了解赤潮加劇和消亡的過程。由于太陽高度角對于地表的太陽輻照度影響較大[18]，從而影響地表尤其是水體的表觀輻亮度[19]，進而影響分割閾值的確定，因此對2014年5月28日8景GOCI影像計算得到的RI結果采用OTSU方法分別計算了各自的分割閾值，如表2所示。分別提取每幅影像中的赤潮，選取兩個無云覆蓋的區域，分別統計區域內8景影像中的赤潮面積，其面積變化結果見圖3。

圖3 2014年5月28日感興趣區域的赤潮面積變化

表2 2014年5月28日8景影像的分割閾值

區域1和區域2的赤潮面積均呈先增大后減小的變化趨勢。區域1的赤潮面積較大，一天之內的面積變化比較劇烈。赤潮面積從8:00左右的2 411 km2開始增加，10:00左右達到一天之中的最大值，面積為2 565.25 km2，之后赤潮面積開始減少，15:00達到最小值，面積為1 258.25 km2，減少了51%。區域2中赤潮面積從8:00的715.25 km2增加到10:00的772 km2，15:00左右赤潮面積降為384.25 km2，降低了約50.2%。從面積變化上可以看出，赤潮在一天內的變化十分明顯，10:00左右赤潮最為嚴重，分布面積最大，這可能與浮游植物的垂直遷移有很大關系。當光照充足時，浮游植物會上浮到海水表面吸收光照，進行光合作用，產生營養物質，當光照減弱后，浮游植物會下潛至海底吸收營養。但10:00之后光照仍然充足，赤潮面積反而開始下降。根據氣泡調節浮力機制[20]，日出之后，浮游植物上浮至海面進行光合作用，由于小氣泡的作用細胞浮力增加，浮游植物量在10:00的時候形成上層高峰，中午強日照下的光合產物使得細胞內部壓強增大導致氣泡破裂，細胞失去浮力造成浮游植物下沉，因此赤潮面積在10:00之后開始下降。根據項斯端等[21]藍藻型富營養化湖泊藻量的晝夜變化規律，藻量水平在一天中會出現兩個峰值：第一次峰值出現在上午10:00，中午12:00至14:00藻量大幅下降并出現低谷；下午藻量水平上升，18:00至20:00之間出現第二次高峰，之后藻量又開始下降。但由于GOCI衛星在15:00之后不再采集影像，所以此規律還需進一步研究證明。

4 赤潮發生過程的水色參數變化分析

赤潮往往是由某種藻類生物異常繁殖導致，只有當各種理化條件基本滿足某種赤潮藻生理、生態需求時，該種藻類才可能形成赤潮。赤潮的發生是多因素共同作用的結果，與赤潮發生區域的環境密切相關。赤潮生物的生長發育和繁殖都要從環境中吸收營養物質和能量，同時各個階段又受周圍環境因素的制約。赤潮會對環境因素造成影響，使得水色參數在赤潮前后發生變化，因此，赤潮與環境之間相互影響。

根據海洋災害公報，2014年5月15日—8月7日，河北秦皇島灤河口至遼寧綏中海域發生大規模赤潮，但從2014年6月12日的葉綠素a影像和RI影像中并未提取到大范圍赤潮，因此，可以認為此次大范圍赤潮在6月12日左右就已基本消亡。本研究選取一些赤潮發生區域云量較少的影像對赤潮導致的海洋水色參數變化進行分析，按照公報顯示的赤潮記錄，4月28日、5月5日、5月8日和5月12日赤潮還未爆發，5月15日、5月26日和5月28日為赤潮快速生長繁殖時期，6月9日、6月12日、6月13日、6月23日為大規模赤潮消亡時期。由于云層的覆蓋，選取這些時段也只能提取到部分赤潮區域的連續信息。

4.1 漫衰減系數Kd(490)

水體漫衰減系數是描述輻照度隨海水深度增加而衰減的重要海洋光學參數，可以反映水體的環境狀況。如圖4(a)所示，赤潮發生前Kd(490)在0.2 m-1左右，赤潮發生后Kd(490)開始增大，赤潮發生第一天達到0.29 m-1，之后數值迅速增大，5月26日達到最大值0.68 m-1，直至6月12日赤潮基本消亡，Kd(490)值減小到0.09 m-1，此后Kd(490)值基本穩定在0.1 m-1左右。在清澈水體中，水體下行輻射隨海水深度增加而衰減的趨勢減小且出現穩定的狀態，因此漫衰減系數較??；當赤潮發生后，下行輻照度衰減程度較大，導致漫衰減系數增大。

4.2 懸浮物濃度

懸浮物濃度對浮游植物種群生長影響是多種機制疊加的結果，較低濃度的懸浮物有利于赤潮藻類的生長，但高濃度的懸浮物具有抑制作用[22]。由圖4(b)可見，5月上旬，懸浮物濃度達2 g·m-3,伴隨著赤潮的生長繁殖，濃度不斷升高，5月28日上升到3.6 g·m-3，之后濃度下降至1～1.5 g·m-3波動。從赤潮發生前后懸浮物的濃度來看，懸浮物濃度總體較低，赤潮藻類的生長并沒有受到懸浮物的影響，因此，隨著赤潮的生長，懸浮物濃度升高，赤潮消亡后，濃度隨之下降。

圖4 各水色參數時間序列圖

4.3 CDOM

CDOM為水體內一類重要的光吸收物質，其濃度和組成能改變水下光場，對水體生態系統有明顯影響。從圖4(c)分析，CDOM變化也遵循先增大后減小的趨勢，且在赤潮最嚴重時達到最大，赤潮消亡，CDOM也隨之減小。但整體變化較小，基本在0.12～0.24 m-1變化。浮游植物的大量生長聚集會伴隨一些植物自身降解，為海水中CDOM提供了充足來源，促使CDOM數值上升。

4.4 海表面溫度

赤潮爆發需要適宜溫度，同時，大量的浮游植物聚集也會對海表面溫度造成影響。但圖4(d)的海表面溫度變化情況并沒有顯示出海表面溫度隨著赤潮的生長消亡發生變化。從4月28日開始，海表面溫度一直上升，原因是隨著季節變化和氣溫的不斷升高海表面溫度也不斷上升，沒有體現出受赤潮影響導致的變化。不過，此次赤潮的主要種類是抑食金球藻，適宜的生長溫度為15～25 ℃，而5月26日的海表面溫度達15 ℃，促進了赤潮的大規模爆發。

4.5 歸一化熒光高度

赤潮爆發時，浮游植物大量聚集，葉綠素a濃度異常增高，其水體光譜存在一個熒光峰，可以用歸一化熒光高度來定量表示(圖4(e))。5月8日前歸一化熒光高度值約為0.35，赤潮爆發后浮游生物量增加，歸一化熒光高度值上升，6月9日歸一化熒光高度值和赤潮爆發前基本持平，之后一直維持在0.37左右。歸一化熒光高度情況基本反映葉綠素a濃度隨赤潮發生的變化情況，盡管葉綠素a濃度與歸一化熒光高度不成正比，但赤潮的生長消亡也會導致葉綠素a濃度發生變化。

綜合各種水色參數變化，赤潮的爆發會導致Kd(490)、懸浮物濃度、CDOM和歸一化熒光高度出現不同程度的升高，盡管本研究中的海表面溫度并沒有隨之變化，但多項研究表明，赤潮的發生會導致海表面溫度的輕微上升[23-24]。Kd(490)、懸浮物和CDOM這3種水色要素在赤潮發生后的數值低于赤潮發生前，這可能是由于赤潮的生長發育消耗了水體中的各種物質，從而達到了凈化水體的作用，導致數值較赤潮發生前低。

5 結論

赤潮是由多種因素作用的結果，同時，赤潮也會引起多種海洋參數發生變化，為赤潮的預測預警提供了依據。

1) 在閾值分割方法上，OTSU方法在葉綠素a影像和赤潮指數影像上提取的赤潮面積精度最高，K-means聚類方法次之，最大熵法的精度最差。在赤潮監測方法上，葉綠素a濃度圖像3種分割方法和赤潮指數圖像的平均相對誤差分別為22%和7.9%，赤潮指數圖像的提取結果更準確。因此，選擇OTSU方法結合赤潮指數圖像進行赤潮提取。

2) 以區域1和區域2的赤潮面積為例研究了赤潮的晝變化情況，赤潮面積存在明顯的晝變化特征，兩個區域的赤潮面積在一天的變化均先增大后減少，在10:00面積達到最大值，這可能與浮游植物的垂直遷移有關。

3)Kd(490)、懸浮物濃度、CDOM、歸一化熒光高度等水色參數均隨赤潮的發生而增高，赤潮消亡后數值也隨之下降。海表面溫度在赤潮發生周期一直上升，季節變化導致赤潮引起的溫度上升不明顯，海洋水色參數的變化可以為赤潮的預警預報提供依據。

由于云層遮擋等原因導致赤潮監測的時間序列不完整，而且GOCI數據的空間分辨率不夠精細，所以赤潮提取面積會有較大誤差，下一步可以考察更高分辨率的數據結合更多的海洋環境參數建立赤潮預警預報模型。同時對于赤潮發生后海洋水色參數的數值低于赤潮發生前的數值的原因還需進一步研究。