基于遙感藻總量和氣象因子的巢湖不同湖區(qū)藻華預(yù)測(cè)*

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

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

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

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

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

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

1.1 研究區(qū)概況

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

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

為獲得歷史時(shí)期的巢湖藻華分布和水柱藻總量，研究使用了搭載在美國(guó)Aqua衛(wèi)星上的中分辨率成像光譜儀(MODerate-resolution imaging spectroradiometer，MODIS/Aqua)掃描獲取的2003-2020年日尺度多光譜影像，數(shù)據(jù)來源于NASA(https://oceancolor.gsfc.nasa.gov/).原始數(shù)據(jù)利用SeaDAS軟件空間重采樣為250 m，并進(jìn)行瑞利校正以去除大氣分子的吸收和散射信號(hào)，從而獲得瑞利校正后的遙感反射率(Rrc).Rrc數(shù)據(jù)共包含7個(gè)波段，波長(zhǎng)范圍分別為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ù)進(jìn)行藻華識(shí)別[19]和水柱藻總量反演[17].

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

圖1 巢湖及不同類型實(shí)測(cè)站點(diǎn)地理區(qū)位

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

2 研究方法

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

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

(1)

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

2.2 水柱藻總量遙感估算

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

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

(2)

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

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

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

(3)

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

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

(4)

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

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

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

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

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

3)湖區(qū)藻華暴發(fā)概率Logistic模型構(gòu)建.統(tǒng)計(jì)分析藻華/非藻華暴發(fā)日對(duì)應(yīng)的藻總量和各氣象因子值，判斷藻華暴發(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)度的檢驗(yàn)標(biāo)準(zhǔn)(越臨近0.5效果越好)，并采用Wald統(tǒng)計(jì)量檢驗(yàn)各自變量對(duì)藻華暴發(fā)的影響程度(值越大相對(duì)影響越大).

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

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).風(fēng)速季節(jié)差異較小，整體呈冬春季高、夏秋季低，但同一月份內(nèi)風(fēng)速變異很大，這反映出巢湖風(fēng)速日變化劇烈(圖3c).總磷濃度在夏秋季高、冬春季低；西巢湖8月總磷濃度((0.2±0.05)mg/L)能達(dá)1月((0.1±0.04)mg/L)的2倍，中巢湖總磷濃度季節(jié)變異弱于西巢湖，東巢湖季節(jié)變異更弱；值得注意的是，總磷濃度高的夏秋季其變異也大.不同環(huán)境因子的季節(jié)變異共同驅(qū)動(dòng)巢湖藻華的季節(jié)變異.

不同環(huán)境因子的湖區(qū)差異不一致.不同湖區(qū)的氣候態(tài)月平均氣溫和水汽壓等差異很??；雖然風(fēng)速湖區(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)境因子空間差異不明顯，但在一些具體日期仍差異明顯，尤其是風(fēng)速.

圖3 不同湖區(qū)氣溫、水汽壓和風(fēng)速氣候態(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)明顯空間差異，這可能與巢湖強(qiáng)水動(dòng)力導(dǎo)致的水藻橫向輸運(yùn)和強(qiáng)烈混合有關(guān).由遙感反演的2003-2019年不同日期的水柱藻總量計(jì)算的不同湖區(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所示，當(dāng)平均藻總量高于100 mg/m2時(shí)，不同湖區(qū)的藻總量表現(xiàn)出明顯差異，這可能是不同湖區(qū)藻華暴發(fā)差異的一個(gè)重要原因.

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

不同湖區(qū)水柱藻總量季節(jié)變異顯著.平均藻總量在冬季低值時(shí)小于10 mg/m2，而在夏季高值時(shí)能達(dá)200 mg/m2(圖5a).由2003-2019年日尺度藻總量可知：高藻總量不會(huì)持續(xù)很久，大多時(shí)間湖區(qū)平均藻總量為60 mg/m2上下.利用Census X-12方法對(duì)2003-2019年的月平均藻總量進(jìn)行時(shí)間序列分解，可得：巢湖不同湖區(qū)藻總量年內(nèi)呈“雙峰、雙谷”的季節(jié)變異特征，主次雙峰出現(xiàn)于8月和3月前后，主次雙谷出現(xiàn)于12月和6月前后(圖5b).對(duì)西巢湖，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)時(shí)間序列分解后的月藻總量；(c)藻華或非藻華日氣溫

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

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

不同于水柱藻總量的“雙峰、雙谷”季節(jié)特征(3.2節(jié))，巢湖不同湖區(qū)、月份的藻華暴發(fā)日占比均呈明顯“單峰、單谷”特征，尤其是藻華暴發(fā)嚴(yán)重的西巢湖(圖6).由氣候態(tài)月平均結(jié)果可知：不同湖區(qū)藻華暴發(fā)日占比9月份達(dá)到峰值，1月達(dá)到谷值；在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ā)日占比變異大(標(biāo)準(zhǔn)差大)，這反映出巢湖不同年份的藻華暴發(fā)情況存在顯著差異，尤其是藻華暴發(fā)嚴(yán)重的西巢湖和中巢湖(圖6).

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

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

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

(6)

式中，P1、P2和P3分別為西、中和東巢湖的藻華暴發(fā)概率，乘以100%則取值范圍為0～100%，值越大表示湖區(qū)藻華暴發(fā)概率越大.利用模型訓(xùn)練之外的剩余1/3匹配樣本檢驗(yàn)?zāi)Ｐ途龋喝绻?6)得到的值 > 50%則預(yù)測(cè)為藻華暴發(fā)，反之為不暴發(fā)；對(duì)西、中和東巢湖，式(6)預(yù)測(cè)藻華暴發(fā)與否的精度分別為90%(N=110)、85%(N=168)和89.5%(N=161).

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

3.5 模型應(yīng)用案例分析

雖然巢湖藻華暴發(fā)主要發(fā)生在夏秋季，但冬春季也時(shí)有藻華暴發(fā)(圖6).為檢驗(yàn)式(6)對(duì)不同季節(jié)藻華暴發(fā)的預(yù)測(cè)能力，本文對(duì)比分析了2020年冬春季和夏秋季巢湖藻華遙感結(jié)果和預(yù)測(cè)的藻華暴發(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)，本文進(jìn)一步計(jì)算了上述不同日期、不同湖區(qū)的藻華暴發(fā)概率，如表1所示.冬春季2月9日和12日的藻華暴發(fā)概率由西湖區(qū)向東湖區(qū)遞增，對(duì)藻華嚴(yán)重的東巢湖計(jì)算的暴發(fā)概率分別為86.6%和85.0%，而藻華少的西巢湖分別為7.7%和19.5%；夏秋季6月17日和7月25日的藻華暴發(fā)概率由西湖區(qū)向東湖區(qū)遞減，對(duì)藻華嚴(yán)重的西巢湖計(jì)算的暴發(fā)概率分別為95.3%和98.2%(表1).因此，不同湖區(qū)藻華暴發(fā)概率計(jì)算結(jié)果和遙感監(jiān)測(cè)結(jié)果整體吻合，在中巢湖存在一定誤差(圖7，表1).中巢湖藻華暴發(fā)概率預(yù)測(cè)誤差偏大可能有兩方面原因：一是模型輸入的氣象因子采用的是合肥站和巢湖站的平均值，與實(shí)際值可能存在一定偏差；二是模型輸入沒有考慮風(fēng)速，而風(fēng)速也是巢湖藻華暴發(fā)與否的一個(gè)重要影響因子[14,17].

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

4 討論

4.1 巢湖藻華暴發(fā)概率預(yù)測(cè)

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

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

4.2 巢湖藻華預(yù)測(cè)精度影響因素

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

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

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

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

4.3 對(duì)其他湖泊藻華預(yù)測(cè)的啟示

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

5 結(jié)論

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

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