劉 東,張 民，沈 明，曹志剛，齊天賜，邱銀國，段洪濤,3

(1:中國科學院南京地理與湖泊研究所，中國科學院流域地理學重點實驗室，南京 210008)(2:中國科學院南京地理與湖泊研究所，湖泊與環(huán)境國家重點實驗室，南京 210008)(3:西北大學城市與環(huán)境學院，西安 710127)

湖泊為人類生產(chǎn)生活提供飲用水、娛樂場所和棲息地等，是人類賴以生存的基礎(chǔ)[1].然而，受人類排污和氣候變化的影響，全球越來越多的湖泊呈現(xiàn)富營養(yǎng)化，一個典型結(jié)果就是藻華暴發(fā)呈增加態(tài)勢[2].藻類大量繁殖后，藻顆粒腐爛分解會消耗溶解氧造成水體缺氧，藍藻等還會釋放藻毒素，因此藻含量是評價湖泊富營養(yǎng)化程度或水質(zhì)好壞的一個關(guān)鍵指標，并通常用葉綠素a濃度(Chl.a)表征藻含量[3-5].關(guān)于富營養(yǎng)湖泊藻含量、藻華時空變異等，目前已有大量研究報道[2,6]，但為合理防治藻華帶來的各種生態(tài)環(huán)境危害，還需提前預測藻華暴發(fā)情況，相關(guān)方面的報道還較少.湖泊藻華暴發(fā)是藻增殖至一定生物量后，在適宜氣象條件下的產(chǎn)物[6-7].對富營養(yǎng)太湖，Shi等[7]報道氣溫和總磷是藻類增殖的主要影響因素，而氣壓和風速是藻華暴發(fā)的主控因子.因此，在已知藻含量和氣象條件時，可進行藻華暴發(fā)概率預測.基于耦合水動力藻含量物理模型，Li等[8]先模擬出藻含量，再結(jié)合氣象因子預測太湖未來3天的藻華概率.基于物理模型的藻華預測機理清晰，但模型輸入?yún)?shù)繁多、構(gòu)建難度極大，目前大多研究是利用不同氣象因子構(gòu)建藻華概率預測的經(jīng)驗模型，選用的氣象因子包括氣溫、氣壓、降雨量和風速等[9-12].這些研究要么顯式限定在藻華暴發(fā)頻繁的夏秋季[12]，要么通過輸入氣溫隱式限定藻華暴發(fā)于高溫季節(jié)[9-11]，不適用于預測不同季節(jié)尤其是冬季藻華，究其原因是缺乏藻含量輸入.

衛(wèi)星遙感具有高時空分辨率優(yōu)勢，越來越多被應用于湖泊藻含量動態(tài)監(jiān)測.針對單個湖泊或湖泊群，已經(jīng)發(fā)展了波段比值、波段差值和分段擬合等多種表層藻含量的遙感算法[2,6-7].然而，富營養(yǎng)湖泊藻顆粒會隨風速等環(huán)境條件在水柱內(nèi)垂向遷移：低風速時會上浮至表層，高風速時在水柱內(nèi)呈垂向均一分布[13].Xue等[14]報道了富營養(yǎng)巢湖水柱藻含量剖面隨風速等變化呈現(xiàn)均一、冪函數(shù)衰減、指數(shù)衰減或高斯分布類型.因此，當水柱藻含量分布不均時，遙感監(jiān)測到的表層或真光層藻含量并不能真實反映湖泊富營養(yǎng)狀況.為更準確監(jiān)測湖泊藻含量，近年來一些研究開始考慮藻類在水柱內(nèi)的不均一分布，并進行水柱藻總量遙感，通常用單位面積水域內(nèi)的Chl.a垂向總含量表征水柱藻總量[14-16]；所采用方法大體分為兩類：一類是先構(gòu)建表層藻含量與水柱藻總量的經(jīng)驗關(guān)系式，再由遙感的表層藻含量估算藻總量，即面向結(jié)果方法[15-16]，該方法已應用于我國滇池和巢湖[15-16]；另一類是先實現(xiàn)水柱不同深度藻含量的函數(shù)化描述，再不同深度積分計算藻總量，即面向過程方法，該方法已應用于我國太湖和巢湖[17].

以富營養(yǎng)巢湖為研究區(qū)，本文首先分析了巢湖不同湖區(qū)藻華暴發(fā)與否日的藻總量和氣象因子差異，然后構(gòu)建了不同湖區(qū)藻華暴發(fā)概率的Logistic預測模型，用于計算不同氣溫、氣壓、風速和藻總量情形下的藻華暴發(fā)概率.綜合考慮藻總量和氣象因子，所構(gòu)建的Logistic模型可用于預測巢湖不同季節(jié)、不同湖區(qū)的藻華暴發(fā)概率，對其藻華防治具有重要現(xiàn)實意義.

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

巢湖位處長江中下游，湖心經(jīng)緯度分別為117°34′E和31°32′N(圖1)，水面積約770 km2，平均水深約3.0 m，是我國第五大淡水湖[15].巢湖是合肥市的重要水源地，養(yǎng)育著流域約1080萬人[17].環(huán)巢湖入湖河流主要有南淝河、派河、柘皋河、白石天河和杭埠河等，湖水能通過裕溪河排出并流入長江干流；其中，南淝河和派河水污染嚴重，攜帶大量營養(yǎng)鹽入湖，杭埠河入湖流量最大、水質(zhì)較好[18].受河流入湖污染物影響，巢湖水體富營養(yǎng)化程度由西至東逐漸降低，根據(jù)湖體富營養(yǎng)狀態(tài)和生態(tài)環(huán)境特征，通常將巢湖分為西巢湖、中巢湖和東巢湖[15,17].

1.2 衛(wèi)星遙感和站點實測數(shù)據(jù)

為獲得歷史時期的巢湖藻華分布和水柱藻總量，研究使用了搭載在美國Aqua衛(wèi)星上的中分辨率成像光譜儀(MODerate-resolution imaging spectroradiometer，MODIS/Aqua)掃描獲取的2003-2020年日尺度多光譜影像，數(shù)據(jù)來源于NASA(https://oceancolor.gsfc.nasa.gov/).原始數(shù)據(jù)利用SeaDAS軟件空間重采樣為250 m，并進行瑞利校正以去除大氣分子的吸收和散射信號，從而獲得瑞利校正后的遙感反射率(Rrc).Rrc數(shù)據(jù)共包含7個波段，波長范圍分別為459～479 nm(Rrc(469))、545～465 nm(Rrc(555))、620～672 nm(Rrc(645))、841～891 nm(Rrc(859))、1230～1250 nm(Rrc(1240))、1628～1652 nm(Rrc(1640))和2105～2155 nm(Rrc(2130)).為避免內(nèi)陸渾濁水體精確大氣校正而帶來更大誤差，研究直接使用Rrc遙感數(shù)據(jù)進行藻華識別[19]和水柱藻總量反演[17].

為判斷水柱藻剖面類型和構(gòu)建藻華預測模型，研究使用了湖區(qū)附近合肥站和巢湖站(圖1)的氣象日值數(shù)據(jù)，指標包括風速、氣溫、水汽壓、氣壓、降水量和日照時數(shù).數(shù)據(jù)來源于國家氣象科學數(shù)據(jù)中心(http://www.cma.gov.cn/)，時間范圍為2003-2020年.為反映巢湖不同湖區(qū)水質(zhì)，研究使用了湖體20個水質(zhì)站點(圖1)的月尺度總磷數(shù)據(jù)，水樣通常在每月上旬采集，時間跨度為2003-2018年，數(shù)據(jù)來源于安徽省巢湖管理局(http://chglj.hefei.gov.cn/).

圖1 巢湖及不同類型實測站點地理區(qū)位

此外，研究還使用了2013-2017年在巢湖開展的7個航次獲取的106個Chl.a剖面(圖1).站點均勻分布于全湖，但在藻華水域會加測以更多獲取藻華水域的Chl.a剖面；此外，在一些點位還進行了時間序列觀測以獲取水柱Chl.a剖面時間動態(tài)特征；由于采樣船定位不準、隨水漂移等，站點位置以GPS記錄為準，且不同航次點位會出現(xiàn)重合(圖1).對每個站位，采集表層、0.1、0.2、0.4、0.7、1.0、1.5、2.0和3.0 m共9個深度的水樣，然后參照葉綠素測定國家行業(yè)標準(SL 88-2012)用Whatman GF/F(孔徑為0.7 mm)濾膜過濾水樣得到Chl.a樣品[17].實驗室內(nèi)，用90%的丙酮萃取Chl.a，并用Shimadzu UV2600分光光度計測量Chl.a濃度.更多關(guān)于剖面Chl.a的獲取過程，請參考Liu等[17].

2 研究方法

2.1 藻華遙感識別和暴發(fā)日占比計算

對于MODIS/Aqua獲取的日尺度Rrc數(shù)據(jù)(1.2節(jié))，將Rrc(1640)> 0.0215的水體判斷為云覆蓋區(qū)域[20].對無云覆蓋水體，研究計算了浮游藻類指數(shù)(floating algae index,FAI)[21]：

(1)

式中，Rrc(645)、Rrc(859)和Rrc(1240)分別表示中心波長為645、859和1240 nm 3個波段的反射率，“859-645”和“1240-859”為對應中心波長的差值[21].FAI反映的是一個像元覆蓋的250 m×250 m水域的藻華綜合情況，為避免混合像元造成藻華面積高估，參照Duan等[6]使用閾值0.02提取純藻華像元，即：FAI> 0.02的像元對應藻華水體，否則為非藻華水體.根據(jù)2003-2019年日尺度藻華遙感結(jié)果，研究進一步計算了巢湖不同像元、不同月份藻華暴發(fā)日占比(藻華暴發(fā)天數(shù)/總無云觀測天數(shù)×100%).由于云覆蓋存在時空變異，藻華暴發(fā)日占比(%)比藻華暴發(fā)天數(shù)能更好地反映不同月份藻華暴發(fā)情況.

2.2 水柱藻總量遙感估算

富營養(yǎng)湖泊水柱內(nèi)的藻會隨氣象條件等發(fā)生垂向遷移，低風速時會上浮聚集在水表形成藻華，而風速大于一定閾值時會呈均勻混合[13].由于遙感只能觀測到水表或真光層內(nèi)藻含量，為遙感反演水柱藻總量可通過兩種途徑：面向結(jié)果或面向過程方法[15-17].基于實測數(shù)據(jù)，Liu等[17]構(gòu)建了巢湖水柱藻總量遙感的面向過程方法，精度優(yōu)于面向結(jié)果方法，尤其是對水表藻類聚集水域.因此，本文采用Liu等[17]提出的面向過程方法遙感反演巢湖藻總量，主要步驟如下：

1)表層Chl.a濃度遙感反演.基于富營養(yǎng)太湖的表層Chl.a遙感算法[7]，并用巢湖的星地同步數(shù)據(jù)進行重新率定，Liu等[17]給出了巢湖的Chl.a遙感算法(N= 196，R2= 0.59).算法均方根誤差和平均相對誤差分別為25.97 mg/m3和-5.61%，可同時反演藻華/非藻華水域Chl.a[7]，公式如下：

(2)

式中，Ratio是一個自定義的歸一化指數(shù)，用以表征水體藻含量，藻含量越多Ratio值越低.

2)藻類剖面分布函數(shù)定義與參數(shù)化.為減少待定系數(shù)量，Liu等[17]通過分析巢湖106個實測Chl.a剖面，認為可將Chl.a剖面概化為垂向均一和冪函數(shù)衰減兩類(圖2)，如式(3)：

圖2 巢湖不同站點的實測Chl.a剖面(均一分布剖面用灰色表示，冪函數(shù)衰減剖面用不同顏色表示)

(3)

式中，Chl.a(z)為水深z處的Chl.a濃度；C0為遙感得到的表層Chl.a濃度；n1和n2為待定系數(shù).系數(shù)n2與表層Chl.a濃度存在顯著對數(shù)關(guān)系，故可由遙感的表層Chl.a濃度計算n2；在已知n2和表層(水深為0.01 m)Chl.a濃度時，可計算得到系數(shù)n1；更多關(guān)于n1和n2的計算過程，請參考Liu等[17].

3)藻類剖面分布類型及水柱藻總量遙感.對任一遙感像元，首先用二叉決策樹判斷藻類剖面分布類型：如果FAI> 0.02則為冪衰減分布；如果FAI≤ 0.02 且表層Chl.a≤ 10 mg/m3，則為均一分布；當FAI≤ 0.02且Chl.a> 10 mg/m3時，如果風速≤ 2.75 m/s則為冪衰減分布，否則為均一分布.然后，通過步驟(2)參數(shù)化藻類剖面分布函數(shù)，并結(jié)合水深(depth，水位減湖泊DEM得到)積分得到單位水域面積的水柱藻總量：

(4)

式中, 對不同深度Chl.a濃度積分得到的單位面積水柱藻總量單位為mg/m2.對富營養(yǎng)巢湖，上述方法遙感估算水柱藻總量的平均絕對百分比誤差為31.03%，更多信息請參考Liu等[17].

2.3 湖區(qū)藻華暴發(fā)概率預測模型構(gòu)建

Logistic模型又稱Logistic回歸分析，以事件發(fā)生與否為因變量，以連續(xù)或離散型因子為自變量，將二元因變量與自變量之間的聯(lián)系轉(zhuǎn)換為因變量條件概率同自變量之間的關(guān)系問題，在醫(yī)療衛(wèi)生和災害預測等領(lǐng)域應用較多[10-12].本文以2003-2019年的遙感水柱藻總量和氣象因子日平均值為自變量，以藻華發(fā)生與否為因變量，構(gòu)建了巢湖西、中、東不同湖區(qū)藻華暴發(fā)概率的Logistic預測模型，主要步驟如下：

1)湖區(qū)藻華/非藻華暴發(fā)日確定.參照已發(fā)表的藻華判別方法[11]，研究將MODIS/Aqua遙感監(jiān)測到的藻華面積占湖區(qū)總面積比 > 5%的日期確定為藻華暴發(fā)日；同時，為避免云覆蓋影響，將湖區(qū)無云覆蓋面積比 > 70%且藻華面積比 ≤ 5%的日期確定為非藻華暴發(fā)日.理論上，可以用上述確定的所有樣本構(gòu)建模型，但由于遙感判斷藻華存在不確定性，藻華像元少的藻華暴發(fā)日和藻華像元多的非藻華暴發(fā)日可能是由于遙感誤判.為規(guī)避遙感誤判的影響，同時確保足夠樣本量，研究只選取藻華像元多的前30%作為藻華暴發(fā)日樣本，賦值為1；同時，只選取藻華像元少的前30%作為非藻華暴發(fā)日樣本，賦值為0.

2)湖區(qū)藻華/非藻華暴發(fā)日氣象因子匹配.根據(jù)步驟(1)確定的藻華/非藻華暴發(fā)日，提取對應的日平均水柱藻總量、氣溫和水汽壓等氣象因子值.其中，平均水柱藻總量為各湖區(qū)不同像元遙感結(jié)果的平均值；對氣象因子，西巢湖用合肥站的值，東巢湖用巢湖站的值，而中巢湖用兩個站的均值(圖1).

3)湖區(qū)藻華暴發(fā)概率Logistic模型構(gòu)建.統(tǒng)計分析藻華/非藻華暴發(fā)日對應的藻總量和各氣象因子值，判斷藻華暴發(fā)與否日的差異，篩選出湖區(qū)藻華暴發(fā)概率估算的自變量最優(yōu)組合，然后基于SPSS構(gòu)建湖區(qū)藻華暴發(fā)概率Logistic模型，形式如下：

(5)

式中，P為藻華暴發(fā)概率，用0～1區(qū)間值表示；x1、x2…xk為輸入變量；β0、β1、β2…βk為擬合系數(shù).采用Cox & Snell和Nagelkerke偽決定系數(shù)R2為模型擬合優(yōu)度的檢驗標準(越臨近0.5效果越好)，并采用Wald統(tǒng)計量檢驗各自變量對藻華暴發(fā)的影響程度(值越大相對影響越大).

式(5)計算得到的藻華暴發(fā)概率表示湖區(qū)藻華暴發(fā)的可能性，值越大表示藻華暴發(fā)的可能性越大.基于歷史數(shù)據(jù)對不同湖區(qū)構(gòu)建了如式(5)的藻華暴發(fā)概率Logistic模型后，假設不同湖區(qū)藻總量在短期內(nèi)穩(wěn)定(藍藻生長周期是8～10 d[22])，則可以將遙感監(jiān)測到的近幾日平均藻總量和氣象預報資料輸入模型，計算得到未來幾日的藻華暴發(fā)概率預測值.

3 結(jié)果分析

3.1 巢湖環(huán)境特征

巢湖環(huán)境因子存在明顯的季節(jié)變異特征.氣溫和水汽壓夏季(6-8月)高、冬季(12-2月)低，年內(nèi)呈明顯單峰變異特征：7月氣溫為(28.86±2.83)℃，水汽壓為(30.97±2.67)hPa；1月氣溫為(3.10±3.13)℃，水汽壓為5.69±11.88 hPa(圖3a～b).風速季節(jié)差異較小，整體呈冬春季高、夏秋季低，但同一月份內(nèi)風速變異很大，這反映出巢湖風速日變化劇烈(圖3c).總磷濃度在夏秋季高、冬春季低；西巢湖8月總磷濃度((0.2±0.05)mg/L)能達1月((0.1±0.04)mg/L)的2倍，中巢湖總磷濃度季節(jié)變異弱于西巢湖，東巢湖季節(jié)變異更弱；值得注意的是，總磷濃度高的夏秋季其變異也大.不同環(huán)境因子的季節(jié)變異共同驅(qū)動巢湖藻華的季節(jié)變異.

不同環(huán)境因子的湖區(qū)差異不一致.不同湖區(qū)的氣候態(tài)月平均氣溫和水汽壓等差異很小；雖然風速湖區(qū)差異也不明顯，但西部湖區(qū)的值整體略低于東部湖區(qū)(圖3a～c).西巢湖的氣候態(tài)月平均總磷濃度明顯高于中巢湖，中巢湖比東巢湖高；在4月，西巢湖總磷濃度((0.16±0.06)mg/L)是中巢湖((0.1±0.04)mg/L)的1.6倍，是東巢湖((0.08±0.03)mg/L)的2倍.值得注意的是：雖然不同湖區(qū)的氣候態(tài)月平均環(huán)境因子空間差異不明顯，但在一些具體日期仍差異明顯，尤其是風速.

圖3 不同湖區(qū)氣溫、水汽壓和風速氣候態(tài)平均值：(a～c)按月份；(d～f)按藻華或非藻華日

3.2 巢湖水柱藻總量遙感結(jié)果

圖4為巢湖氣候態(tài)月平均水柱藻總量.水柱藻總量在巢湖南部近岸淺水區(qū)呈明顯低值，藻總量低的冬季月份甚至低于10 mg/m2.而在西巢湖北部和東巢湖東部近岸水域，呈現(xiàn)出較高藻總量，如2月和4月(圖4).此外的其他水域，氣候態(tài)月平均水柱藻總量并未表現(xiàn)明顯空間差異，這可能與巢湖強水動力導致的水藻橫向輸運和強烈混合有關(guān).由遙感反演的2003-2019年不同日期的水柱藻總量計算的不同湖區(qū)氣候態(tài)月平均值也未表現(xiàn)明顯差異：西巢湖平均藻總量為(66.95±27.6)mg/m2，中巢湖為(63.07±26.84)mg/m2，東巢湖為(67.11±25.38)mg/m2.雖然不同湖區(qū)氣候態(tài)月平均水柱藻總量差異小，但一些日期的不同湖區(qū)水柱藻總量仍表現(xiàn)出明顯差異.如圖5a所示，當平均藻總量高于100 mg/m2時，不同湖區(qū)的藻總量表現(xiàn)出明顯差異，這可能是不同湖區(qū)藻華暴發(fā)差異的一個重要原因.

圖4 巢湖2003-2019年氣候態(tài)月平均水柱藻總量

不同湖區(qū)水柱藻總量季節(jié)變異顯著.平均藻總量在冬季低值時小于10 mg/m2，而在夏季高值時能達200 mg/m2(圖5a).由2003-2019年日尺度藻總量可知：高藻總量不會持續(xù)很久，大多時間湖區(qū)平均藻總量為60 mg/m2上下.利用Census X-12方法對2003-2019年的月平均藻總量進行時間序列分解，可得：巢湖不同湖區(qū)藻總量年內(nèi)呈“雙峰、雙谷”的季節(jié)變異特征，主次雙峰出現(xiàn)于8月和3月前后，主次雙谷出現(xiàn)于12月和6月前后(圖5b).對西巢湖，8月和3月的氣候態(tài)平均藻總量分別為(90.12±44.92)和(78.24±45.8)mg/m2，12月和6月的氣候態(tài)平均藻總量分別為(51.96±21.7)和(61.7±25.7)mg/m2.巢湖不同湖區(qū)藻總量的季節(jié)變異規(guī)律與氣溫、總磷的季節(jié)變化特征存在一定相似性(3.1節(jié)).

圖5 不同湖區(qū)平均值：(a)日藻總量；(b)時間序列分解后的月藻總量；(c)藻華或非藻華日氣溫

3.3 藻華暴發(fā)日占比時空變異

雖然巢湖不同湖區(qū)藻總量差異較小(3.2節(jié))，但不同月份的藻華暴發(fā)日占比存在明顯空間差異.整體而言，巢湖藻華暴發(fā)日占比呈明顯“西高東低”特征，尤其是藻華暴發(fā)嚴重的夏秋季(圖6).由2003-2019年不同月份的藻華暴發(fā)日占比可知：冬春季藻華暴發(fā)少，不同湖區(qū)藻華暴發(fā)日占比基本都<10%，且空間差異不明顯；夏秋季藻華暴發(fā)明顯多，尤其是西巢湖，8-9月份一些水域的藻華暴發(fā)日占比甚至 > 30%；此外，西巢湖夏秋季不同月份藻華暴發(fā)日占比高值區(qū)存在不一致，夏季6-8月高值區(qū)集中在北部水域，而秋季9-11月在西南水域也會出現(xiàn)高值區(qū)(圖6).

不同于水柱藻總量的“雙峰、雙谷”季節(jié)特征(3.2節(jié))，巢湖不同湖區(qū)、月份的藻華暴發(fā)日占比均呈明顯“單峰、單谷”特征，尤其是藻華暴發(fā)嚴重的西巢湖(圖6).由氣候態(tài)月平均結(jié)果可知：不同湖區(qū)藻華暴發(fā)日占比9月份達到峰值，1月達到谷值；在9月和1月，西巢湖的藻華暴發(fā)日占比分別為27.9%±6.62%和1.85%±3.37%，中巢湖分別為13.28%±6.09%和1.58%±1.35%，東巢湖分別為3.47%±2.16%和1.37%±0.75%(圖6).此外，不同月份藻華暴發(fā)日占比變異大(標準差大)，這反映出巢湖不同年份的藻華暴發(fā)情況存在顯著差異，尤其是藻華暴發(fā)嚴重的西巢湖和中巢湖(圖6).

圖6 2003-2019年氣候態(tài)月平均藻華暴發(fā)日占比

3.4 不同湖區(qū)藻華暴發(fā)概率預測模型

為確定巢湖不同湖區(qū)藻華暴發(fā)概率Logistic預測模型的自變量，研究首先分析了藻華暴發(fā)與否日不同因子的差異.由藻華/非藻華暴發(fā)日與環(huán)境因子的歷史匹配數(shù)據(jù)可得(2.3節(jié))，不同湖區(qū)藻華/非藻華暴發(fā)日的平均氣溫、平均水汽壓和平均藻總量存在明顯差異，尤其是藻總量(圖3d～f，5c)：對西巢湖，藻華暴發(fā)日的平均氣溫、水汽壓和藻總量分別為22.75℃、21.14 hPa和85.73 mg/m2，而非藻華暴發(fā)日的值分別為7.33℃、6.63 hPa和61.65 mg/m2；對中巢湖，藻華暴發(fā)日的平均氣溫、水汽壓和藻總量分別為21.18℃、20.41 hPa和82.1 mg/m2，而非藻華暴發(fā)日的值分別為10.27℃、8.76 hPa和54.59 mg/m2；對東巢湖，藻華暴發(fā)日的平均氣溫、水汽壓和藻總量分別為16.74℃、16.47 hPa和81.37 mg/m2，而非藻華暴發(fā)日的值分別為12.82℃、10.72 hPa和56.62 mg/m2.因此，本研究以日平均氣溫、水汽壓和水柱藻總量為自變量，通過隨機選取2/3歷史匹配樣本(西巢湖N=220, 中巢湖N=336, 東巢湖N=322)構(gòu)建了不同湖區(qū)的藻華暴發(fā)概率Logistic預測模型，如式(6)：

(6)

式中，P1、P2和P3分別為西、中和東巢湖的藻華暴發(fā)概率，乘以100%則取值范圍為0～100%，值越大表示湖區(qū)藻華暴發(fā)概率越大.利用模型訓練之外的剩余1/3匹配樣本檢驗模型精度：如果式(6)得到的值 > 50%則預測為藻華暴發(fā)，反之為不暴發(fā)；對西、中和東巢湖，式(6)預測藻華暴發(fā)與否的精度分別為90%(N=110)、85%(N=168)和89.5%(N=161).

模型擬合優(yōu)度檢驗結(jié)果表明：式(6)很好地擬合了建模數(shù)據(jù)，且水汽壓和水柱藻總量對不同湖區(qū)藻華暴發(fā)均有顯著影響(P<0.05).對西巢湖，Cox & Snell和NagelkerkeR2分別為0.55和0.74，日平均氣溫、水汽壓和水柱藻總量3個自變量影響均顯著(P<0.05)；Wald值分別為6.84、9.02和6.59，即水汽壓影響最大.對中巢湖，Cox & Snell和NagelkerkeR2分別為0.32和0.49，水汽壓和藻總量影響顯著(P<0.05)；Wald值分別為22.69和39.03，即藻總量影響更大.對東巢湖，Cox & Snell和NagelkerkeR2分別為0.23和0.45，3個自變量影響均顯著(P<0.05)；Wald值分別為7.1、6.24和53.01，即藻總量影響明顯大.

3.5 模型應用案例分析

雖然巢湖藻華暴發(fā)主要發(fā)生在夏秋季，但冬春季也時有藻華暴發(fā)(圖6).為檢驗式(6)對不同季節(jié)藻華暴發(fā)的預測能力，本文對比分析了2020年冬春季和夏秋季巢湖藻華遙感結(jié)果和預測的藻華暴發(fā)概率.基于MODIS/AquaRrc數(shù)據(jù)，研究得到了2020年冬春和夏秋季少云覆蓋的4天的藻華和水柱藻總量分布，如圖7所示.巢湖在這4天內(nèi)均暴發(fā)了藻華，但不同日期藻華空間分布差異大：冬春季的2月9日和12日，藻華主要分布在中、東巢湖，總面積分別為419.00和424.69 km2；夏秋季的6月17日和7月25日，藻華主要分布在西巢湖，總面積分別為97.31和175.56 km2(圖7).不同日期遙感反演的全湖平均藻總量分別為89.29、87.43、67.50和108.64 mg/m2，整體表現(xiàn)為高藻華面積、高藻總量，且藻華暴發(fā)水域水柱藻總量通常也較高(圖7).

圖7 2020年不同日期遙感藻華與藻總量

將日平均氣溫、水汽壓和藻總量代入式(6)，本文進一步計算了上述不同日期、不同湖區(qū)的藻華暴發(fā)概率，如表1所示.冬春季2月9日和12日的藻華暴發(fā)概率由西湖區(qū)向東湖區(qū)遞增，對藻華嚴重的東巢湖計算的暴發(fā)概率分別為86.6%和85.0%，而藻華少的西巢湖分別為7.7%和19.5%；夏秋季6月17日和7月25日的藻華暴發(fā)概率由西湖區(qū)向東湖區(qū)遞減，對藻華嚴重的西巢湖計算的暴發(fā)概率分別為95.3%和98.2%(表1).因此，不同湖區(qū)藻華暴發(fā)概率計算結(jié)果和遙感監(jiān)測結(jié)果整體吻合，在中巢湖存在一定誤差(圖7，表1).中巢湖藻華暴發(fā)概率預測誤差偏大可能有兩方面原因：一是模型輸入的氣象因子采用的是合肥站和巢湖站的平均值，與實際值可能存在一定偏差；二是模型輸入沒有考慮風速，而風速也是巢湖藻華暴發(fā)與否的一個重要影響因子[14,17].

表1 2020年不同日期藻華暴發(fā)概率模型計算結(jié)果

4 討論

4.1 巢湖藻華暴發(fā)概率預測

本文所構(gòu)建的公式(6)可應用于巢湖不同湖區(qū)次日(t+1)藻華暴發(fā)概率預測.巢湖藻華暴發(fā)的主控藻種是藍藻，其生長周期為8～10 d[22]，即可以假設短期內(nèi)湖區(qū)藻總量不變，然后將今日(t)的遙感藻總量作為t+1日的藻總量輸入模型.值的注意的是，由于衛(wèi)星遙感易受云雨覆蓋影響，經(jīng)常會導致t日的藻總量沒有遙感結(jié)果；此種情況下，為預測t+1日的藻華暴發(fā)概率，可以使用前幾日的湖區(qū)藻總量遙感平均值為輸入.對于t+1日的氣溫和水汽壓，可以使用歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecasts，ECMWF)發(fā)布的模式預測數(shù)據(jù)集(時間分辨率為3 h，空間分辨率為0.125度)，數(shù)據(jù)集由中央氣象局統(tǒng)一分發(fā)至各省市氣象局，巢湖區(qū)域數(shù)據(jù)可以從安徽省氣象局獲取.

根據(jù)富營養(yǎng)狀態(tài)，本文只將巢湖粗略分為3個湖區(qū)，而實際應用中會關(guān)注更小區(qū)域的藻華暴發(fā)情況，如自來水廠取水口、臨近市鎮(zhèn)的近岸水域等.為實現(xiàn)這些重點區(qū)域的藻華暴發(fā)概率預測，可以對巢湖進行更精細化劃分，然后基于歷史氣象因子和遙感藻總量，根據(jù)2.3節(jié)方法，分別構(gòu)建各重點水域藻華暴發(fā)概率Logistic預測模型.此外，當假設某一水域藻總量不變，也可預測未來幾日的藻華暴發(fā)概率，但不同區(qū)域水藻交換導致的藻顆粒水平遷移會使模型預測精度隨區(qū)域變小和預測時間變長而降低.

4.2 巢湖藻華預測精度影響因素

水柱藻總量是巢湖藻華暴發(fā)的關(guān)鍵影響因素，尤其是對東巢湖，藻總量遙感精度會直接影響藻華預測結(jié)果(3.4節(jié)).本文以近幾日的遙感藻總量輸入模型計算t+1日的藻華概率，一個重要依據(jù)是巢湖水動力強，進而可以假設三大湖區(qū)的日平均藻總量一致、短期內(nèi)湖區(qū)藻總量也不變(3.3節(jié)).由圖8可知，雖然三大湖區(qū)平均藻總量在大部分時間基本一致，但仍在一些日期存在明顯差異；如果這些日期不同湖區(qū)水藻交換明顯，則會導致輸入模型的藻總量與實際情況不符，進而影響藻華預測結(jié)果.因此，為提高藻華預測精度，需要保證模型輸入藻總量的真實性：一方面可通過算法優(yōu)化提高藻總量遙感精度[17]；另一方面可使用靜止衛(wèi)星數(shù)據(jù)以縮短藻總量遙感與藻華預測的時間差.

圖8 2003-2019年巢湖不同湖區(qū)遙感日平均藻總量對比

營養(yǎng)鹽、日照時長和氣溫等會影響藻類增殖[6-7]，風速、水汽壓和降雨等會影響藻類遷移[7,14]，因此這些環(huán)境因子也會影響藻華暴發(fā).對富營養(yǎng)化湖泊，Webster等[13]報道存在風速閾值，低于該閾值時藻類會大量上浮至表層形成藻華，而高于該閾值時水柱內(nèi)藻類會趨于均勻混合.陳莉瓊等[12]報道日照時長、前三日氣壓差、風速和降雨量是影響富營養(yǎng)滇池藻華形成的關(guān)鍵氣象因素，可用于預測滇池藻華暴發(fā)概率.然而，營養(yǎng)鹽缺乏日尺度資料，風速和降雨在一天內(nèi)不同小時會劇烈變動，藻華暴發(fā)與否日照時長沒明顯差異，因此本文的藻華概率預測模型只輸入了氣溫和水汽壓.不同湖區(qū)模型通過調(diào)整氣溫和水汽壓系數(shù)，充分利用兩者在藻華/非藻華暴發(fā)日差異，以更好預測藻華暴發(fā)概率.此外，由于缺乏湖區(qū)氣象數(shù)據(jù)，使用的是附近合肥站和巢湖站的氣象資料，但不同位置的氣象因子可能存在明顯時空差異.因此，為進一步提高巢湖藻華暴發(fā)預測精度，一方面需使用湖體高頻氣象資料；另一方面，模型構(gòu)建時需同時考慮影響藻類生長和分布的風速和日照時長等.

本文是以水柱積分Chl.a表征藻總量，進而構(gòu)建Logistic模型預測巢湖不同季節(jié)的藻華暴發(fā)概率.但不同季節(jié)藻華暴發(fā)時藻總量存在差異(表1)、優(yōu)勢藻種會不一致[23]，藻類不同生長階段Chl.a含量也會不一致，這些都會影響用Chl.a表征藻總量的準確度，進而影響藻華預測模型的通用性.為提高巢湖不同季節(jié)的藻華暴發(fā)概率預測精度，有必要考慮季節(jié)差異，分季節(jié)構(gòu)建藻華預測模型.

4.3 對其他湖泊藻華預測的啟示

藻華廣泛發(fā)生于全球湖泊[2]，為合理防治藻華危害，需對不同區(qū)位、不同類型湖泊藻華提前預測.本文研究的是富營養(yǎng)巢湖的藻華預測，藍藻是主控性藻種，在藻總量、水汽壓和氣溫滿足一定條件后會上浮至表層形成藻華.因此，對本文構(gòu)建的藻華暴發(fā)概率預測模型進行參數(shù)重新率定，可推廣應用于全球其他藍藻主控的富營養(yǎng)湖泊，如太湖.而對非藍藻主控的富營養(yǎng)湖泊，由于藻華暴發(fā)的控制機制不同，藻華暴發(fā)預測的敏感因子也不同，不能簡單通過參數(shù)率定將本文構(gòu)建的模型應用于其藻華暴發(fā)預測.但湖泊藻華暴發(fā)都是藻類增殖至一定生物量后，在適宜氣象條件下的產(chǎn)物[6-7].因此，對非藍藻主控的其他富營養(yǎng)湖泊，仍可參照2.3節(jié)方法，以藻總量和氣象因子為輸入構(gòu)建相應藻華暴發(fā)概率Logistic預測模型.

5 結(jié)論

以我國長江中下游典型富營養(yǎng)巢湖為研究區(qū)，本文首先分析了不同湖區(qū)水柱藻總量、藻華暴發(fā)日占比及藻華/非藻華暴發(fā)日各氣象因子差異，然后構(gòu)建了巢湖不同湖區(qū)藻華暴發(fā)概率的Logistic預測模型，模型輸入?yún)?shù)包括遙感的水柱藻總量、日平均氣溫和日平均水汽壓.巢湖不同湖區(qū)藻華暴發(fā)存在顯著時空變異特征，模型預測結(jié)果與遙感歷史監(jiān)測結(jié)果基本一致；模型也能預測2020年冬春季和夏秋季藻華暴發(fā)，基本滿足藻華預測的應用需求.雖然，基于遙感獲取的前幾日水柱藻總量和預測的日平均氣溫及水汽壓可以實現(xiàn)巢湖不同湖區(qū)藻華暴發(fā)概率預測，但藻華暴發(fā)和氣象因子的高時空動態(tài)導致的多源數(shù)據(jù)時空匹配難會降低模型預測精度，今后還需進一步結(jié)合湖體高時空分辨的多源氣象因子進行模型優(yōu)化和藻華預測.

致謝：野外實驗數(shù)據(jù)的獲取得到了中國科學院南京地理與湖泊研究所湖泊環(huán)境遙感團隊各位老師和研究生的支持，在此表示衷心感謝.