• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    1980s以來巢湖藻華物候時空變化遙感分析*

    2023-01-13 07:21:48曹志剛馬金戈齊天賜段洪濤
    湖泊科學(xué) 2023年1期
    關(guān)鍵詞:水華巢湖藍藻

    袁 俊,曹志剛,馬金戈,沈 明,齊天賜,段洪濤

    (1:西北大學(xué)城市與環(huán)境學(xué)院,西安 710127) (2:西北大學(xué),陜西省地表系統(tǒng)與環(huán)境承載力重點實驗室,西安 710127) (3:中國科學(xué)院南京地理與湖泊研究所,中國科學(xué)院流域地理學(xué)重點實驗室,南京 210008)

    湖泊是地表極其重要的水資源,為流域工農(nóng)業(yè)發(fā)展和人類生活提供穩(wěn)定而清潔的水源,是人類活動與發(fā)展的熱點區(qū)域[1]。氣候變化及流域土地利用變化,加速湖泊富營養(yǎng)化[2],有害藻華暴發(fā)加劇[3],全球許多大型湖泊都出現(xiàn)嚴(yán)重的水華,如美國伊利湖,加拿大溫尼伯湖,中國太湖、巢湖和滇池等[4]。藍藻水華發(fā)生導(dǎo)致了一系列的湖泊水生態(tài)環(huán)境問題,顯著削弱了湖泊的生態(tài)服務(wù)功能和價值,嚴(yán)重制約了區(qū)域社會經(jīng)濟可持續(xù)發(fā)展[5-6]。以往研究基于藍藻水華的暴發(fā)頻率分析藻華的時空變化,在此基礎(chǔ)上,進一步了解藍藻生長周期的關(guān)鍵時間點,研究其暴發(fā)狀態(tài)的時序變化,對探究藻華發(fā)生機理、改善水質(zhì)也具有借鑒意義[7]。生物學(xué)領(lǐng)域?qū)⑸镩L期適應(yīng)自然環(huán)境周期性變化,形成與此相適應(yīng)的生長發(fā)育節(jié)律的現(xiàn)象稱之為物候[8],藻類屬于浮游植物,浮游植物物候呈現(xiàn)了在季節(jié)性環(huán)境變化下浮游植物生長繁殖的規(guī)律性變化,能夠有效地反映浮游植物的生長和湖泊的生長條件[9]。不同于生物學(xué)意義上的“物候”變化,遙感領(lǐng)域通常將藍藻水華暴發(fā)面積或葉綠素a等參數(shù)達到一定閾值時的起始和持續(xù)時間等時序變化特征稱之為藍藻水華的物候特征[7,10-11]。理解藍藻水華物候?qū)τ跉夂蜃兓捻憫?yīng)和其在水環(huán)境中的變化機制,有助于準(zhǔn)確掌握湖泊水質(zhì)和生態(tài)環(huán)境變化特征及輔助提升藍藻水華預(yù)測預(yù)警。

    相較于傳統(tǒng)的觀測手段,光學(xué)衛(wèi)星遙感因其速度快、范圍廣、監(jiān)測周期短,已成為湖泊藍藻水華的重要監(jiān)測手段[12]。藍藻水華在水面聚集,引起近紅外波段反射率抬升,形成植被的反射光譜,這是光學(xué)遙感圖像提取藍藻水華的理論基礎(chǔ)。目前,光學(xué)遙感識別藍藻水華的方法分為3類:光譜指數(shù)提取法[13]、監(jiān)督分類法[14]、機器學(xué)習(xí)分類方法[15]。監(jiān)督分類方法需要人工干預(yù),受主觀影響大,機器學(xué)習(xí)方法能實現(xiàn)復(fù)雜地物的提取要求,但是需要大量的訓(xùn)練樣本。相對而言,基于藻華光譜特征構(gòu)建光譜指數(shù),結(jié)合閾值提取,機理性強且提取快速,是衛(wèi)星遙感監(jiān)測湖泊藍藻水華的主流方法。

    藍藻水華暴發(fā)周期長,變化快,以Moderate-Resolution Imaging Spectroradiometer(MODIS)(Terra:1999年至今;Aqua:2002年至今)為代表的中分辨率海洋水色傳感器被廣泛應(yīng)用于水華面積和物候監(jiān)測?;贛ODIS 時間序列數(shù)據(jù),Liniger等[16]發(fā)現(xiàn)了南極洲默茨冰川與B09B冰山2010年撞擊后其水華暴發(fā)開始時間延遲、持續(xù)時間縮短;Rohr等[17]探討了南大洋浮游植物水華物候的變化機制;Wynne等[18]探究了薩吉諾灣和伊利湖近20年的藻華物候,分析了薩吉諾灣和伊利湖藻華物候的影響因子。雖然MODIS等傳感器重訪周期快、波段多,但是衛(wèi)星發(fā)射于1999年后,實際上,太湖、巢湖等一些藍藻水華暴發(fā)頻繁的區(qū)域,2000年之前富營養(yǎng)化已比較嚴(yán)重且藍藻水華已經(jīng)暴發(fā)[19-20],補充這一時間段的藍藻水華時空變化信息有助于提升對湖泊藻華變化規(guī)律的認識。Landsat系列傳感器,包括Landsat 4/5 TM(Thematic Mapper,1982-2011 年)、Landsat 7 ETM+ (Enhanced Thematic Mapper plus,1999-2022 年)、Landsat 8/9 OLI (Operational Land Imager,2013年至今),是持續(xù)時間最長的地球資源衛(wèi)星之一。Landsat可以彌補MODIS在2000年之前觀測的空白,但其時間分辨率低(16 d),觀測藍藻水華頻率確實有所不足[21]。如何綜合利用Landsat系列傳感器和MODIS數(shù)據(jù),考慮不同數(shù)據(jù)源觀測藍藻水華的差異,聯(lián)合監(jiān)測湖泊藍藻水華變化,將使得觀測1980s以來湖泊藍藻水華時空變化成為可能。

    本文的研究目標(biāo)是使用Landsat和MODIS監(jiān)測巢湖1987-2020年間藻華暴發(fā)情況和物候時空變化,解析巢湖藍藻水華近40年的時空規(guī)律。具體內(nèi)容包括:(1)基于光譜指數(shù)和閾值分割方法,研究Terra MODIS和Landsat提取藍藻水華時空變化的一致性和差異;(2)獲取1987-2020年巢湖藍藻水華暴發(fā)面積時空分布,分析1987-2020年巢湖藻華物候特征(暴發(fā)開始時間、暴發(fā)結(jié)束時間、暴發(fā)持續(xù)時間)的變化規(guī)律;(3)基于巢湖總氮(TN)、總磷(TP)濃度和氣溫、風(fēng)速、降雨等環(huán)境因子,分析藻華物候指標(biāo)長時間變化的原因。

    1 數(shù)據(jù)與方法

    1.1 研究區(qū)概況

    巢湖(31°25′28″~31°43′28″N,117°16′54″~117°51′46″E,圖1)位于安徽省中部,屬合肥市管轄,是我國五大淡水湖之一,水域面積約760 km2,東西長55 km、南北寬21 km、湖岸線周長176 km,平均水深2.89 m,多年平均水位為8.03 m,容積20.7×108m3[22]。湖水主要靠地面徑流補給,入湖河流約33條,其中主要出入湖河流有9條,包括南淝河、十五里河、派河、杭埠河、柘皋河、雙橋河、兆河、白石天河、裕溪河[23]。近年來,隨著巢湖流域工業(yè)化以及城市化的發(fā)展,巢湖水體呈現(xiàn)嚴(yán)重富營養(yǎng)狀況,藍藻水華頻繁暴發(fā),對巢湖流域的生態(tài)環(huán)境產(chǎn)生了嚴(yán)重的影響[24]。

    圖1 巢湖及其分區(qū)、氣象站點位置Fig.1 Location of Lake Chaohu and its zoning and meteorological station

    1.2 衛(wèi)星數(shù)據(jù)與預(yù)處理

    1.2.1 Landsat系列衛(wèi)星數(shù)據(jù) 從美國地質(zhì)調(diào)查局(United States Geology Survey, USGS, https://earthexplorer.usgs.gov)網(wǎng)站下載了巢湖上空1987-2020年無云或少云情況下的Landsat L1T級影像共289景,其中Landsat 5 TM 205景、Landsat 7 ETM+ 24景、Landsat 8 OLI 60景(表1)。下載過程中通過USGS提供的Landsat快視圖判斷巢湖上空是否被云覆蓋,對存在云影響的圖像進行去云操作。Landsat 7數(shù)據(jù)在2003年后出現(xiàn)條帶,數(shù)據(jù)質(zhì)量下降, 2003年后未使用此數(shù)據(jù);同時Landsat 5在2011年停止運行,Landsat 8在2013年2月發(fā)射,此期間沒有有效Landsat數(shù)據(jù)。研究使用ACOLITE軟件(https://github.com/acolite)對TM、ETM+和OLI數(shù)據(jù)進行處理得到瑞利散射校正反射率Rrc(Rayleigh-corrected reflectance)[25]。

    1.2.2 MODIS衛(wèi)星數(shù)據(jù) MODIS搭載在Terra和Aqua兩顆衛(wèi)星上,具有較高時間分辨率,一天內(nèi)可以覆蓋全球兩次(部分赤道低緯度地區(qū)除外)[26]。MODIS的陸地波段(前7個波段,光譜范圍為443~2130 nm)具有較高的空間分辨率(250~500 m)且在渾濁湖泊不飽和,因此本研究僅使用這7個波段提取藍藻水華。考慮到Landsat系列衛(wèi)星過境時間為當(dāng)?shù)貢r間上午10∶15左右,因此本文選用了上午過境的MODIS/Terra,以獲取時間相對一致的湖泊觀測。從美國國家宇航局(NASA)獲取了2000-2020年巢湖天氣狀況良好的MODIS/Terra影像1919景,基本涵蓋了2000-2020年全年各月份,特別是水華最為嚴(yán)重的夏秋季節(jié)(表1)。MODIS數(shù)據(jù)利用SeaDAS 7.5.3軟件去除臭氧吸收和分子瑞利散射的影響,獲得Rrc數(shù)據(jù)[27-29]。

    表1 巢湖1987-2020年Landsat/MODIS衛(wèi)星影像數(shù)據(jù)的時間分布*Tab.1 Temporal distribution of Landsat/MODIS data in Lake Chaohu used in this study from 1987 to 2020

    1.2.3 陸地和云掩膜 首先使用Landsat和MODIS影像快視圖進行目視觀察,排除巢湖上空多云覆蓋的圖像;然后,對存在少云覆蓋的圖像采用單波段閾值進行云掩膜處理。云掩膜閾值分別是MODIS_Rrc(1240)≥0.0235[30]和Landsat_Rrc(SWIR)>0.018[31]。采用歸一化水體指數(shù)(normalized difference water index,NDWI)進行水體邊界提取,根據(jù)水體在影像中的分布選擇閾值將水體與其他地物區(qū)分開來[32]。最后,研究將水體向內(nèi)掩膜3個像元,避免陸地臨近效應(yīng)對水體信號的影響。

    1.3 藍藻水華提取算法

    1.3.1 浮游藻類指數(shù)(FAI) 常用的光譜指數(shù)有NDVI、增強型植被指數(shù)法(enhanced vegetation index,EVI)和浮游藻類指數(shù)(floating algae index,F(xiàn)AI)等。NDVI對大氣干擾處理有限,受云層影響較大,不易區(qū)分高渾濁水體,同時當(dāng)藍藻密度較大時,NDVI會出現(xiàn)飽和現(xiàn)象[33]。EVI具有較好的抗大氣干擾能力,可抑制水體背景噪聲,但易將藍藻水華與水生植被混淆[34]。相比于NDVI、EVI等算法,F(xiàn)AI算法對環(huán)境和觀測條件(氣溶膠類型和厚度、太陽/觀測幾何和太陽光)的變化不敏感,并且不易受云層影響,可以有效避免對藍藻水華的誤判,已廣泛應(yīng)用到湖泊藍藻水華的監(jiān)測。

    FAI指數(shù)實質(zhì)上是一種波段減法,即通過構(gòu)建紅光、近紅外和短波紅外的基線來反映藻類在水體表面聚集的光譜變化特征,其計算公式為:

    FAI=Rrc(λNIR)-R′rc(λNIR)

    (1)

    (2)

    其中,MODIS參與計算的波段為λRed=645 nm、λNIR=859 nm、λSWIR=1240 nm;TM參與計算的波段為λRed=660 nm、λNIR=839 nm、λSWIR=1678 nm;ETM+參與計算的波段為λRed=661 nm、λNIR=835 nm、λSWIR=1650 nm;OLI參與計算的波段為λRed=655 nm、λNIR=865 nm、λSWIR=1609 nm[35]。

    1.3.2 藍藻水華閾值確定 使用Hu等[11]在太湖提出的梯度閾值方法提取藍藻水華。對于Landsat系列數(shù)據(jù),使用FAI<-0.01和FAI>0.02排除純凈水體和純藻華像元,生成藻-水混合圖像的FAI梯度直方圖??紤]到藻華-非藻華邊界處FAI值梯度變化大,因此單幅圖像梯度變化最大處所對應(yīng)的FAI值即是單景影像的藻華閾值。最后每景Landsat影像FAI閾值的均值減去其兩倍標(biāo)準(zhǔn)差作為所有Landsat的藻華提取閾值,研究確定的Landsat藻華提取FAI閾值為-0.0089。

    MODIS數(shù)據(jù)由于空間分辨率較低,研究采用藻華像元生長算法[36](algae pixel growing algorithm,APA)進一步提高藍藻水華提取精度。從MODIS圖像確定巢湖的純藻華FAI閾值為0.05,純水體FAI閾值為-0.004,對水和藍藻混合的像元,選用APA精確計算藻華混合像元內(nèi)的藻華面積,估算湖泊全水域藍藻水華實際面積[36]。像元生長算法是通過判定衛(wèi)星影像中的藻華“生長點像元”,采用臨近像元相關(guān)和逐漸擴展的思路,計算出混合像元的藻華蓋度,這里所說的蓋度是指某一混合像元中藻華完全覆蓋的面積占該像元面積的百分比[37],具體計算過程見參考文獻[38]。研究通過統(tǒng)計覆蓋度不為0的藻華像元的個數(shù)CMODIS,計算整個巢湖的藍藻水華面積:

    Bloomarea=0.25×0.25×CMODIS

    (3)

    1.4 藍藻水華物候指標(biāo)

    目前,業(yè)內(nèi)對湖泊水華發(fā)生的定義仍沒有形成統(tǒng)一共識。在以往評估遠洋環(huán)境浮游植物物候的工作中,水華事件開始的一個常見定義是葉綠素a濃度(Chl.a)高于背景中值濃度的5%[10]。考慮巢湖藍藻水華暴發(fā)的季節(jié)性和年際變化周期特征以及《水華遙感與地面監(jiān)測評價技術(shù)規(guī)范(征求意見稿)》中“太湖、巢湖水華狀況判斷暫行辦法”定義藍藻水華物候分析的周期為每年春季(本年3-5月)開始至冬季(12月-次年2月)結(jié)束。暴發(fā)開始時間為每年3月份開始藻華面積第1次為全湖面積的5%時的日期;暴發(fā)結(jié)束時間為次年2月之前藻華面積為全湖的5%時的最后日期;藻華暴發(fā)結(jié)束日期減去暴發(fā)開始日期即為藻華持續(xù)時間。關(guān)于1%、5%和10%閾值識別的巢湖藻華物候指標(biāo)結(jié)果的對比分析見3.1節(jié)。

    Landsat傳感器和MODIS的時空分辨率存在差異。針對二者空間分辨率差異在監(jiān)測藻華分布的影響,本文選取了2000-2020年41景巢湖同步MODIS/Terra和Landsat影像,分析其提取藍藻水華面積的一致性問題(3.1節(jié))。1987-1999年藻華物候特征只能由Landsat系列數(shù)據(jù)獲取,考慮到Landsat的時間分辨率較低,2000年之前的物候參數(shù)使用3年為周期進行統(tǒng)計(1987-1989、1990-1992、1993-1995、1996-1998、1999-2000年),巢湖藻華暴發(fā)開始時間取3年的最小值,暴發(fā)結(jié)束時間和持續(xù)時間取3年的最大值。每周期內(nèi),巢湖每月基本至少2景Landsat影像(表1),雖然2月有影像缺失,但藻華暴發(fā)時間從每年3月開始統(tǒng)計,對結(jié)果影響相對較弱。關(guān)于Landsat影像數(shù)量潛在引起的結(jié)果差異,在3.1.3節(jié)進一步深入討論。

    1.5 環(huán)境驅(qū)動因子

    研究獲取巢湖1987-2020年月度總氮(TN)、總磷(TP)濃度數(shù)據(jù),TN、TP濃度分別由堿性過硫酸鉀消解紫外分光光度法(GB 11894-1989)和鉬酸銨分光光度法(GB 11893-1989)測定[39]。統(tǒng)計計算1987-2020年這兩種水質(zhì)指標(biāo)的年均值,討論分析其對巢湖藻華物候的影響。此外,藍藻水華暴發(fā)受氣象因子(溫度、風(fēng)速、降雨量和日照時長)年際以及季節(jié)性變化影響較大,除了當(dāng)年的氣象因子外,上年年均以及冬季氣象因子也對藻華暴發(fā)有顯著影響[40-41]。研究從中國氣象數(shù)據(jù)中心(http://data.cma.cn)下載了合肥氣象站的日均氣溫、風(fēng)速、降水和日照時長數(shù)據(jù),之后生成年平均及季節(jié)平均結(jié)果,以分析氣象因子與巢湖藻華物候特征的關(guān)系。

    1.6 統(tǒng)計分析指標(biāo)

    本文利用Pearson相關(guān)系數(shù)來分析巢湖藻華物候與環(huán)境因子的相關(guān)性,相關(guān)分析使用線性逐步回歸。采用Theil-Sen median和Mann-Kendall趨勢分析方法,研究巢湖藻華暴發(fā)開始時間、暴發(fā)結(jié)束時間、暴發(fā)持續(xù)時間的長時間變化趨勢,并以P值表示顯著性大小,P<0.05為顯著性水平,P<0.01為極顯著水平。

    2 結(jié)果

    2.1 藍藻水華面積的時序變化

    1987-2020年Landsat和MODIS監(jiān)測結(jié)果 (圖2) 顯示,巢湖藍藻水華規(guī)模逐漸擴大,暴發(fā)面積不斷增加。2000年之前,全巢湖藍藻水華整體暴發(fā)規(guī)模較小,大面積藻華較少,1990年出現(xiàn)峰值507.3 km2,其余年份每年全湖區(qū)藍藻水華暴發(fā)面積峰值都在500 km2以下,且面積峰值集中在夏秋季出現(xiàn)(圖2d)。2000年之后,整個巢湖藍藻水華暴發(fā)面積在2000-2005年呈上升趨勢,2005-2010年較穩(wěn)定,2011年出現(xiàn)最高峰608.4 km2,2012-2014年藻華暴發(fā)面積有所下降,2015年出現(xiàn)次高峰,2016-2020年呈下降趨勢(圖2d)。西巢湖2000年之前藍藻水華暴發(fā)規(guī)模也較小,暴發(fā)峰值出現(xiàn)在1999年,達到216.2 km2,2000年之后,西巢湖藍藻水華面積居高不下 (圖2a)。1987-1999年中巢湖和東巢湖無明顯藻華發(fā)生,其中,中巢湖暴發(fā)峰值出現(xiàn)在1991年(162.7 km2),東巢湖暴發(fā)峰值出現(xiàn)在1990年(204.1 km2)(圖2b,c)。2000年之后,中巢湖藻華面積在2000-2012持續(xù)增加,此后呈現(xiàn)波動式下降的規(guī)律(圖2b);東巢湖藍藻水華面積在2000-2011年逐漸增加,2011年達到暴發(fā)面積次峰值225.5 km2,2012年后暴發(fā)面積整體呈現(xiàn)下降趨勢(圖2c)。

    圖2 巢湖1987-2020年藍藻水華暴發(fā)面積: (a)西巢湖;(b)中巢湖 ;(c)東巢湖 ;(d)全巢湖Fig.2 Area of algal blooms in Lake Chaohu from 1987 to 2020: (a) west lake; (b) central lake; (c) east lake; (d) the whole lake

    過去34年中,全湖以及各湖區(qū)每月藍藻水華平均暴發(fā)面積均呈增加趨勢,其中,全湖增加趨勢最為顯著(slope=0.28,P<0.01),其次是西巢湖(slope=0.14,P<0.01),中巢湖(slope=0.07,P<0.01)以及東巢湖(slope=0.07,P<0.01)增加趨勢較平緩(圖3),分湖區(qū)月平均暴發(fā)面積均呈增加趨勢,可能使全湖整體上月平均面積增加趨勢最為顯著。此外,每年藍藻水華暴發(fā)面積最大值大多集中在6-11月,即夏、秋季節(jié),冬、春季藍藻水華大面積暴發(fā)較少,符合巢湖夏、秋季節(jié)藍藻水華暴發(fā)較為嚴(yán)重,冬、春季節(jié)有所減緩的規(guī)律。

    圖3 巢湖1987-2020年藍藻水華月平均面積長時間變化趨勢: (a)西巢湖;(b)中巢湖;(c)東巢湖; (d)全巢湖Fig.3 Long-term trend of monthly mean area of algal blooms in Lake Chaohu from 1987 to 2020: (a) west lake; (b) central lake ;(c) east lake; (d) the whole lake

    2.2 巢湖藻華物候時空變化趨勢

    本研究以藻華面積大于5%湖泊總面積為藻華判定閾值,結(jié)合1987-2020年Landsat和MODIS數(shù)據(jù)計算出巢湖藍藻水華每年暴發(fā)開始時間、結(jié)束時間和持續(xù)時間 (圖4,表2)??傮w上,巢湖藻華暴發(fā)多開始于4-6月,集中在次年1月和2月結(jié)束。1987-2000年間巢湖藻華暴發(fā)開始時間多在5月,暴發(fā)開始時間最晚的日序數(shù)為154 d,其次是149 d(1993-1995年)、121 d(2000年)和114 d(1987-1989年);2000年之后暴發(fā)開始時間逐步提前至4月甚至3月(圖4,表2)。2000年前,巢湖藻華暴發(fā)結(jié)束時間多集中在1月初,結(jié)束最早的日序數(shù)為362 d(1990-1992年);2000年后,暴發(fā)結(jié)束時間有所延遲,且多集中在2月中下旬結(jié)束,其中,最晚的日序數(shù)為423 d(2016年)(圖4,表2)。1987-2000年,藻華暴發(fā)持續(xù)時間集中在5~7個月,2000年達到10個月;2000年之后,暴發(fā)持續(xù)時間不斷增加,最短為8個月(2020年),最長為12個月(2006年)。

    表2 巢湖1987-2020年藻華物候參數(shù)期及氣象數(shù)據(jù)統(tǒng)計Tab.2 Statistics of algal blooms phenology and meteorological data of Lake Chaohu from 1987 to 2020

    圖4 巢湖1987-2020年Landsat(散點)和MODIS/Terra(折線)藍藻水華暴發(fā)面積以及暴發(fā)開始時間(a)和暴發(fā)結(jié)束時間(b)(括號中標(biāo)注的是年份及藻華開始或結(jié)束的日序,結(jié)束時間出現(xiàn)在來年1-2月,所以,結(jié)束時間的日序數(shù)延續(xù)了前一年的日序數(shù);1987-1999年為3年時間段取物候指標(biāo),年份取中間值)Fig.4 Landsat (scatters) and MODIS/Terra (lines) cyanobacterial blooms area and timing of initiation (a) and timing of termination (b) in Lake Chaohu (Note that the numbers in the parentheses are year and day sequence for the start or end of algal blooms. The day sequence value of the end time is superimposed on the day sequence number of the previous year. 1987-1999 is a three-year time period to take the phenological index, and the year is taken as the middle value)

    總體上,1987-2020年巢湖藻華物候特征的時序變化分為3個階段:①1987-2004年,暴發(fā)開始時間呈顯著提前趨勢(R2=0.66,P<0.01),平均每年提前9.7 d,暴發(fā)結(jié)束時間有所延遲,但不顯著(R2=0.27,P>0.05),平均每年延遲3.9 d,暴發(fā)持續(xù)時間呈顯著增加趨勢(R2=0.67,P<0.01),平均每年增加17.4 d。②2005-2010年,藻華暴發(fā)開始時間呈現(xiàn)顯著延遲趨勢(R2=0.61,P<0.05),平均每年延遲8.9 d,暴發(fā)結(jié)束時間無明顯變化,暴發(fā)持續(xù)時間有所減少,但不顯著(R2=0.15,P>0.05),平均每年減少8.8 d。③2011-2020年,巢湖藻華暴發(fā)開始、結(jié)束和持續(xù)時間變化趨勢均不顯著,其中,暴發(fā)開始時間平均每年提前1.3 d(R2=0.14,P>0.05),暴發(fā)結(jié)束時間平均每年提前5.7 d(R2=0.27,P>0.05),暴發(fā)持續(xù)時間平均每年縮短4.5 d(R2=0.15,P>0.05)(圖5)。

    圖5 巢湖1987-2020年藻華物候變化趨勢: (a)開始時間 ;(b)結(jié)束時間 ;(c)持續(xù)時間(灰色箭頭僅表示為上升或下降趨勢)Fig.5 Long-term trend of algal blooms phenology in Lake Chaohu from 1987 to 2020: (a) the timing of initiation; (b) the timing of termination; (c) the timing of duration (Gray arrows indicate an upward or downward trend)

    2.3 巢湖藻華物候指標(biāo)與環(huán)境因子的關(guān)系

    巢湖藻華物候指標(biāo)的年變化與TN、TP濃度無顯著相關(guān)關(guān)系。1987-2004以及2005-2010年TN、TP濃度均呈下降趨勢,2011-2020年TN、TP濃度略微上升,且從整體上看,1987-2020年TN、TP濃度呈顯著下降趨勢(圖6a)。盡管30多年來,巢湖TN、TP濃度有所下降,但對比《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》 (GB 3838-2002)發(fā)現(xiàn),巢湖1987-2020年平均TN濃度(2.18 mg/L)超過了國家Ⅴ類水體標(biāo)準(zhǔn),TP濃度(0.16 mg/L)位于國家Ⅲ類水體標(biāo)準(zhǔn),巢湖水體富營養(yǎng)化程度一直處于較高水平,滿足藍藻水華生長要求[42],所以,湖體的營養(yǎng)鹽不是影響巢湖藻華物候指標(biāo)的主要因子。

    巢湖藻華不同階段物候指標(biāo)受氣象因子影響不同(圖7)。1987-2004年暴發(fā)開始時間與上年年平均溫度呈顯著負相關(guān)(R2=0.50,P<0.05)(圖7a);暴發(fā)結(jié)束時間和持續(xù)時間與氣象因子均無明顯相關(guān)關(guān)系。研究發(fā)現(xiàn),1987-2004年巢湖藻華暴發(fā)結(jié)束時間和持續(xù)時間都呈延遲趨勢,同時段內(nèi),巢湖年平均溫度顯著上升(R2=0.43,P<0.05)(圖6c),年平均風(fēng)速和冬季平均風(fēng)速均有所下降(圖6b);②2005-2010年,藻華暴發(fā)開始時間與春季降雨量呈顯著正相關(guān)(R2=0.82,P<0.05),與春季平均溫度(R2=0.69,P<0.05)、春季平均日照(R2=0.95,P<0.01)呈顯著負相關(guān)(圖7b~d);暴發(fā)結(jié)束(R2=0.72,P<0.05)和持續(xù)時間(R2=0.93,P<0.01)與年降雨量呈顯著負相關(guān)(圖7e~f);2011-2020年,藻華暴發(fā)開始時間、結(jié)束和持續(xù)時間與氣象因子均無顯著相關(guān)關(guān)系。該時間段內(nèi),藻華暴發(fā)開始時間略微提前,而巢湖上年冬季溫度有所上升(圖6c),暴發(fā)結(jié)束時間提前,持續(xù)時間縮短,相應(yīng)地,冬季平均風(fēng)速不斷增加(圖6b)。

    圖6 1987-2020年巢湖環(huán)境因子變化趨勢:(a)TN和TP濃度; (b)冬季平均風(fēng)速和年平均風(fēng)速;(c)上年冬季平均溫度和年平均溫度Fig.6 Trends of environmental factors in Lake Chaohu from 1987 to 2020: (a) TN and TP; (b) average winter wind speed and annual average wind speed; (c) average temperature in the previous winter and annual average temperature

    圖7 1987-2004年暴發(fā)開始時間與上年年平均溫度的關(guān)系(a); 2005-2010年暴發(fā)開始時間與春季降雨量、平均溫度以及日照時長的關(guān)系(b~d);2005-2010年暴發(fā)結(jié)束和持續(xù)時間與年降雨量的關(guān)系(e~f)Fig.7 Relationship between the start time and the average annual temperature of the previous year from 1987 to 2004 (a); Relationship between the start time and the spring precipitation, temperature and sunshine from 2005 to 2010 (b-d); Relationship between the termination time and duration and the annual precipitation from 2005 to 2010 (e-f)

    3 討論

    3.1 多源衛(wèi)星遙感精度與不確定性

    3.1.1 衛(wèi)星藻華提取的不確定性 基于FAI指數(shù)和閾值分割的湖泊藍藻水華提取有效避免了湖泊水體中泥沙信號和大氣程輻射的影響,提升了精度,但仍有幾個因素會干擾本文藍藻水華衛(wèi)星提取的絕對精度:(1)水陸邊界的臨近效應(yīng):靠近陸地的水體信號強,F(xiàn)AI值也會偏高,導(dǎo)致靠近陸地的像元有時也會被識別為藍藻水華;(2)遙感圖像的云掩膜多使用單波段閾值法,但能夠完全去除云陰影影響的算法尚不成熟,Rrc(1240)>0.023和Rrc(1640)>0.0215會錯誤地將藍藻水華作為“云像元”掩膜去除,過度云掩膜造成有效數(shù)據(jù)丟失[35,43];(3)研究所使用的FAI閾值提取藍藻水華,通常適用于積聚在水體表面的藍藻,即“cyanobacterial scum”,不易監(jiān)測較為分散且濃度較低的藍藻水華,影響了藍藻水華提取精度。

    祁國華等[14]指出不同空間分辨率衛(wèi)星對水華的識別會略有差別,衛(wèi)星影像過境時間存在時間差,水華可能會在幾分鐘內(nèi)發(fā)生微小變化,造成不同衛(wèi)星的提取結(jié)果出現(xiàn)不同程度的差異。2000-2020年間巢湖41景同步MODIS和Landsat圖像提取藻華面積一致性較好(圖8),差異在0.19~39.85 km2之間,其中提取誤差在0~10 km2的占51%,在10~20 km2的占22%,在20~30 km2的占15%,在30~40 km2的占12%。Landsat空間分辨率較高,使得二者提取誤差低值所占比例較大。Landsat和MODIS提取結(jié)果差異最大的在2016年11月14日(MODIS提取面積為79.43 km2,Landsat提取面積為39.58 km2),是由于當(dāng)天MODIS上空存在薄云,造成提取誤差。

    圖8 巢湖MODIS和Landsat同步衛(wèi)星數(shù)據(jù)提取藍藻水華面積一致性對比Fig.8 Comparison of the consistency of algal blooms area from MODIS and Landsat synoptic satellite data in Lake Chaohu

    3.1.2 物候閾值的影響 藻華物候的確定方法多樣,如Soppa等[44]規(guī)定當(dāng)Chl.a濃度值超過中位數(shù)5%時的時間至少持續(xù)兩周時定義為水華開始日期,將Chl.a值下降至閾值以下第1周為水華結(jié)束日期。關(guān)鍵閾值的選取影響物候特征的精度,因此本文對比分析1%閾值和10%閾值的巢湖藻華物候結(jié)果,驗證5%閾值提取物候的精度。1%為閾值即在巢湖藻華面積達到7.6 km2時認為藍藻水華暴發(fā)。MODIS本身空間分辨率較低,該閾值面積僅占有120個像元,提取的物候特征與巢湖藻華暴發(fā)規(guī)律有較大偏差,且年際之間無明顯變化,巢湖整年都在暴發(fā)藻華(圖9)。10%為閾值計算藻華物候,2000年之后藻華暴發(fā)開始時間與5%相差不大,但由于閾值面積達到76 km2,受Landsat時間分辨率和多種不確定性因素的影響,2000年之前有4年(1988、1989、1993和1998年)藻華暴發(fā)開始時間和結(jié)束時間為同一天,且1996年為藻華暴發(fā)的空白期,即僅有的遙感圖像數(shù)據(jù)無法滿足閾值要求,1996年被認定為無藻華暴發(fā)(圖9),因此10%閾值不能準(zhǔn)確刻畫所提取的巢湖藍藻物候特征。

    圖9 10%和1%閾值藻華物候期結(jié)果: (a)暴發(fā)開始時間;(b)暴發(fā)結(jié)束時間Fig.9 10% and 1% threshold algal blooms phenology: (a)the start time; (b)the termination time

    3.1.3 Landsat時間分辨率的影響 Landsat具有較高空間分辨率,但其重訪周期較長(16 d),較少的衛(wèi)星數(shù)據(jù)降低了對巢湖的有效觀測次數(shù),一定程度上影響藻華物候時序變化的精度。為減弱2000年前 Landsat影像稀少對年藻華物候特征的影響,本研究使用3年周期統(tǒng)計巢湖2000年之前的物候參數(shù)。統(tǒng)計周期內(nèi)幾乎每月平均至少2景Landsat遙感影像覆蓋,間接地提升了2000年前Landsat系列數(shù)據(jù)觀測巢湖藻華物候時序變化的時間分辨率。該方法統(tǒng)計的巢湖1987-2020年藻華物候的長時序變化趨勢相對不存在突變(圖5)。為揭示Landsat每3年周期和MODIS數(shù)據(jù)統(tǒng)計藻華物候的時間一致性,我們對比了2000-2020年Landsat和MODIS所提取的藻華暴發(fā)開始時間(圖10)。Landsat每3年周期觀測的藻華暴發(fā)時間雖然比MODIS的結(jié)果略晚一些,但基本上較為接近。同時,二者觀測到的藻華暴發(fā)開始時間的趨勢較為一致,即2000-2010年暴發(fā)開始時間有所延遲,此后藻華暴發(fā)提前。事實上,以往一些研究使用Landsat數(shù)據(jù)對水質(zhì)情況進行分析時也采用了類似的方式[45-47],并揭示出了湖泊水質(zhì)的長時序變化。雖然使用多年周期統(tǒng)計的方式一定程度上彌補了Landsat數(shù)據(jù)不足的缺點,提升了Landsat分析結(jié)果的可信度,但仍需說明的是,2000年前后藻華物候參數(shù)統(tǒng)計的數(shù)據(jù)源不同,且數(shù)據(jù)的時間分辨率存在一定的差異,Landsat低衛(wèi)星數(shù)量統(tǒng)計的年藻華物候特征無法達到與MODIS完全等同的觀測頻率,由此,所得出的藻華物候變化趨勢及其與環(huán)境因子的相關(guān)性不可避免地會存在一定的不確定性。但是Landsat是目前唯一可連續(xù)獲取的2000年前的免費衛(wèi)星數(shù)據(jù),是將藍藻水華監(jiān)測結(jié)果延長15年的最好選擇,同時,探索性聯(lián)合時空分辨率存在差異的多源衛(wèi)星數(shù)據(jù)共同監(jiān)測湖泊藻華物候,可為之后的研究提供一點思路和參考。

    圖10 2000-2020年Landsat藻華暴發(fā)開始時間(每3年)與MODIS的對比Fig.10 Comparison of the start time of algal blooms between Landsat (every 3 years) and MODIS from 2000 to 2020

    3.2 巢湖藍藻水華物候變化的驅(qū)動機制

    巢湖藍藻以微囊藻(Microcystis)和魚腥藻(Anabaena)為主,不同季節(jié)藍藻優(yōu)勢藻藻種不同,5-9月以微囊藻為主,而3-4月和10月后以魚腥藻為主[48],本文“藻華物候”重點在于巢湖藍藻水華暴發(fā)的時間變化,而非生態(tài)上的“物候”的定義,微囊藻、魚腥藻或其他藍藻藻類大量繁殖后在水面聚集均會被識別為藍藻水華,因此,藻種間的交替變化(藻種生長期差異)不會顯著影響巢湖藻華暴發(fā)時間的變化。

    水體中營養(yǎng)鹽濃度升高,藻類吸收營養(yǎng)鹽,迅速繁殖形成有機團聚體,進而在水面形成水華[49]。本研究發(fā)現(xiàn)巢湖TN、TP等營養(yǎng)鹽濃度對于藻華物候變化特征影響較弱,這可能是因為巢湖本身TN、TP等營養(yǎng)鹽濃度一直處于較高水平,滿足藍藻生長要求,因此營養(yǎng)鹽不是改變巢湖藻華暴發(fā)時間變化的主要因子。本研究發(fā)現(xiàn),上年年平均溫度、春季平均溫度越高,日照時長越長,藻華暴發(fā)時間越早。溫度升高、日照時長增加利于藻華生長,使得藻華暴發(fā)時間提前[50]。巢湖1987-2020年上年年平均溫度(R2=0.35,P<0.01)和春季平均溫度(R2=0.40,P<0.01)呈顯著上升趨勢(圖11a),春季平均日照呈不顯著增加趨勢(R2=0.08,P>0.05)(圖11b),同時,1987-2020年巢湖年平均風(fēng)速和冬季平均風(fēng)速均呈下降趨勢,2005年后,平均風(fēng)速都在3 m/s以下(圖6b,表2),溫度升高、日照時長保持穩(wěn)定、風(fēng)速降低利于藻華聚集,形成大面積水華,由此,未來巢湖藍藻水華暴發(fā)時間可能依然呈提前趨勢,暴發(fā)結(jié)束和持續(xù)時間呈延遲趨勢,相應(yīng)地,藍藻水華暴發(fā)規(guī)模變大,周期變長,給巢湖藍藻水華的治理帶來一定困難。

    圖11 巢湖1987-2020上年年平均溫度和春季平均溫度(a)和春季平均日照(b)的變化趨勢 Fig.11 Trends of annual mean temperature in the last year and spring mean temperature (a) and spring mean sunshine (b) from 1987 to 2020 in Lake Chaohu

    研究還發(fā)現(xiàn)不同湖區(qū)影響藍藻水華暴發(fā)的主要因子不同,風(fēng)速對于藍藻水華暴發(fā)也有重要影響。風(fēng)速較低,藍藻容易在水面堆積形成水華,長時間低風(fēng)速天氣,會加劇藍藻水華聚集程度,導(dǎo)致水華面積大幅增加[51],Wang等[52]研究發(fā)現(xiàn)風(fēng)速與滇池藻華發(fā)生呈負相關(guān),并且滇池春季藻華持續(xù)時間明顯短于秋冬季,其認為強風(fēng)是限制藻華聚集的主要因素,滇池春季的高風(fēng)速容易分散藻華使得春季藻華持續(xù)時間較短;Shi等[53]研究指出太湖2003-2017年浮游植物水華暴發(fā)開始時間與春季平均風(fēng)速呈顯著正相關(guān)。

    不同區(qū)域的湖泊或者相同湖泊不同時間受氣象因子的影響程度不同的原因是多方面的。中國的一些淺水湖泊所處地理位置不同,如滇池位于云貴高原,太湖、巢湖位于東部平原,其所處地理位置的氣候條件有明顯差異,氣候差異以及滇池的水質(zhì)特征使得滇池藍藻水華發(fā)生難易與嚴(yán)重程度與太湖、巢湖有所不同[54],溫度和氮、磷營養(yǎng)鹽是導(dǎo)致“三湖”藍藻水華暴發(fā)存在差異的主要因子[55]。對海洋區(qū)域而言,藍藻水華發(fā)生以及浮游植物物候與海洋所處地理位置(高緯度地區(qū)和低緯度地區(qū))、海洋表面溫度、海冰濃度變化、海冰厚度、云層覆蓋、氣象因素也均有一定的相關(guān)性[56-58]。

    3.3 未來巢湖藍藻水華變化的啟示

    藻類水華暴發(fā)是在一定的營養(yǎng)、氣候、水文條件和生態(tài)環(huán)境下形成的藻類過度繁殖和聚集的現(xiàn)象,是水體環(huán)境因子(如 TN、TP、pH值、水動力、溶解氧)和氣象因子(如氣溫、光照、風(fēng)速、降雨量等)綜合作用的結(jié)果,一般認為較高的營養(yǎng)鹽濃度、適宜的溫度、充足的光照以及較低的風(fēng)速條件下,更易暴發(fā)藍藻水華[59]。巢湖流域水體富營養(yǎng)化與合肥市近年來的發(fā)展息息相關(guān),據(jù)安徽省統(tǒng)計年鑒,2000-2020年間,合肥市的人口從2000年的438萬增加到2020年的937萬,而GDP則從324.73億元猛增到10045.72億元。人口的迅速增加、工農(nóng)業(yè)生產(chǎn)的迅速發(fā)展、污水、化肥、農(nóng)藥的大量使用和排放等一系列大規(guī)模的活動嚴(yán)重破壞了巢湖的生態(tài)環(huán)境,巢湖出現(xiàn)了明顯富營養(yǎng)化,制約了該地區(qū)的社會和經(jīng)濟發(fā)展,引起了國家和地方政府的極大關(guān)注,從“六五”至“九五”,圍繞巢湖的富營養(yǎng)化防治,進一步開展了一些管理調(diào)查與對策研究[60]。近年來,國家大力開展“三河三湖”的治理,巢湖的富營養(yǎng)化趨勢得到了明顯的遏制,并且富營養(yǎng)化水平有所下降,但是其營養(yǎng)化水平仍處于高位,足以支持藍藻水華的發(fā)生[23]。巢湖水質(zhì)恢復(fù)是一個流域系統(tǒng)工程,需要科技、資金和全社會的關(guān)心和介入[61]。富營養(yǎng)化及藍藻水華問題的徹底解決,需要流域污染控制與生態(tài)修復(fù)相結(jié)合,不可能一朝一夕完成。通過構(gòu)建藍藻水華預(yù)測體系,預(yù)判水華藍藻發(fā)展?fàn)顩r,采取針對性的應(yīng)急處置措施,可以有效減少藍藻水華帶來的次生災(zāi)害及其對周邊居民的影響。未來,除了研究各種環(huán)境條件下藍藻水華暴發(fā)的規(guī)律,也應(yīng)著重關(guān)注藍藻種類和生理狀態(tài)對于藻華暴發(fā)的影響[62]。

    4 結(jié)論

    1)1987-2020年巢湖藍藻水華暴發(fā)規(guī)模不斷擴大,暴發(fā)面積不斷增加。研究時段內(nèi),全湖區(qū)及其分區(qū)每月藍藻水華暴發(fā)面積均呈增加趨勢,全湖區(qū)增加趨勢最為顯著。此外,每年藍藻水華暴發(fā)面積的最大值大多集中在6-11月,即夏秋季節(jié),冬季和春季巢湖藍藻水華大面積暴發(fā)較少。

    2)在過去34年里,巢湖藍藻水華暴發(fā)時間變化分為3個階段:1987-2004年,巢湖藍藻水華年暴發(fā)開始時間顯著提前,暴發(fā)持續(xù)時間顯著增加;2005-2010年,藻華年暴發(fā)開始時間顯著延遲,但暴發(fā)持續(xù)時間變化不顯著;2011-2020年,巢湖藻華暴發(fā)開始、結(jié)束和持續(xù)時間呈現(xiàn)年際波動,年暴發(fā)開始時間、結(jié)束時間和持續(xù)時間均有所提前,但變化都不顯著。

    3)1987-2020年巢湖藻華物候指標(biāo)年變化與TN、TP年均濃度無顯著相關(guān)性。1987-2004年暴發(fā)開始時間與上年年平均溫度呈顯著負相關(guān),表明上年年平均溫度增加會使得藻華暴發(fā)開始時間提前,2005-2010年,藻華暴發(fā)開始時間與春季平均溫度、春季平均日照呈顯著負相關(guān),與春季降雨量呈顯著正相關(guān),表明春季平均溫度越高、日照越長、降雨量越少,巢湖藻華暴發(fā)時間越早;暴發(fā)結(jié)束和持續(xù)時間與年降雨量呈顯著負相關(guān),表明隨著年降雨量增加,暴發(fā)結(jié)束提前,暴發(fā)持續(xù)時間縮短;2011-2020年,藻華暴發(fā)開始時間、結(jié)束和持續(xù)時間與氣象因子均無顯著相關(guān)關(guān)系。

    猜你喜歡
    水華巢湖藍藻
    藻類水華控制技術(shù)及應(yīng)用
    南美白對蝦養(yǎng)殖池塘藍藻水華處理舉措
    南美白對蝦養(yǎng)殖池塘藍藻水華處理舉措
    巢湖頌歌
    針對八月高溫藍藻爆發(fā)的有效處理方案
    可怕的藍藻
    春季和夏季巢湖浮游生物群落組成及其動態(tài)分析
    巢湖玉卮意蘊長
    大眾考古(2014年7期)2014-06-26 08:00:56
    華能巢湖電廠脫硝系統(tǒng)的改造
    自動化博覽(2014年6期)2014-02-28 22:32:18
    油酸酰胺去除藍藻水華的野外圍隔原位試驗
    母亲3免费完整高清在线观看| www日本在线高清视频| 精品久久久久久久毛片微露脸| 伊人久久大香线蕉亚洲五| 午夜久久久久精精品| 淫秽高清视频在线观看| 激情在线观看视频在线高清| 色视频www国产| 久久精品91无色码中文字幕| 国产蜜桃级精品一区二区三区| 熟妇人妻久久中文字幕3abv| 中文字幕最新亚洲高清| 真人一进一出gif抽搐免费| 伊人久久大香线蕉亚洲五| 少妇的丰满在线观看| www日本黄色视频网| 亚洲一区高清亚洲精品| 男女做爰动态图高潮gif福利片| 久久久久性生活片| 精品国产三级普通话版| 国产激情偷乱视频一区二区| 欧美激情在线99| 亚洲人成网站在线播放欧美日韩| 精品久久久久久久久久久久久| 久久久国产精品麻豆| 亚洲无线观看免费| 国产亚洲欧美98| 丝袜人妻中文字幕| 国产黄色小视频在线观看| 国产午夜精品久久久久久| 久久精品夜夜夜夜夜久久蜜豆| av天堂在线播放| 九九在线视频观看精品| 国产亚洲精品av在线| www.熟女人妻精品国产| 三级毛片av免费| av视频在线观看入口| 黄片小视频在线播放| 国产精品影院久久| 三级男女做爰猛烈吃奶摸视频| 成人欧美大片| 国产激情久久老熟女| 久久午夜亚洲精品久久| 亚洲精品中文字幕一二三四区| АⅤ资源中文在线天堂| 午夜a级毛片| 午夜福利在线在线| 国产精品一区二区精品视频观看| 欧美最黄视频在线播放免费| 久久精品国产99精品国产亚洲性色| 欧美成人一区二区免费高清观看 | 在线免费观看的www视频| 欧美又色又爽又黄视频| 在线十欧美十亚洲十日本专区| 人妻丰满熟妇av一区二区三区| 国模一区二区三区四区视频 | 97超级碰碰碰精品色视频在线观看| 午夜福利视频1000在线观看| av在线天堂中文字幕| www.自偷自拍.com| 国产精品亚洲一级av第二区| 天天一区二区日本电影三级| 99久久精品一区二区三区| 巨乳人妻的诱惑在线观看| 国产亚洲精品久久久久久毛片| 欧美日韩国产亚洲二区| 久久午夜综合久久蜜桃| 成年免费大片在线观看| 给我免费播放毛片高清在线观看| 性色av乱码一区二区三区2| 人妻夜夜爽99麻豆av| 亚洲欧美精品综合一区二区三区| 99精品在免费线老司机午夜| 精品久久蜜臀av无| 久久久国产成人精品二区| 一级作爱视频免费观看| 国产精品久久视频播放| 亚洲狠狠婷婷综合久久图片| 18禁黄网站禁片免费观看直播| 九九在线视频观看精品| 国语自产精品视频在线第100页| 真实男女啪啪啪动态图| 夜夜夜夜夜久久久久| 亚洲精华国产精华精| 亚洲中文av在线| 在线免费观看的www视频| 色在线成人网| 亚洲成人久久性| 成年版毛片免费区| 欧美一级a爱片免费观看看| 亚洲av电影不卡..在线观看| 成年人黄色毛片网站| 精品免费久久久久久久清纯| 哪里可以看免费的av片| 久久久久亚洲av毛片大全| 色哟哟哟哟哟哟| 欧美激情久久久久久爽电影| 国产久久久一区二区三区| 久久草成人影院| 一级作爱视频免费观看| 色精品久久人妻99蜜桃| 亚洲欧美日韩卡通动漫| 国产精品av久久久久免费| 国产又色又爽无遮挡免费看| 我的老师免费观看完整版| 国产午夜福利久久久久久| 日本黄色视频三级网站网址| 日韩中文字幕欧美一区二区| av福利片在线观看| 亚洲一区二区三区不卡视频| 麻豆一二三区av精品| 夜夜躁狠狠躁天天躁| 日韩欧美国产在线观看| 国产久久久一区二区三区| 成年女人毛片免费观看观看9| 熟女人妻精品中文字幕| 精品熟女少妇八av免费久了| 男插女下体视频免费在线播放| 国产成人精品无人区| 在线观看美女被高潮喷水网站 | 欧美日韩瑟瑟在线播放| 香蕉久久夜色| 成人精品一区二区免费| 免费一级毛片在线播放高清视频| 观看免费一级毛片| 色尼玛亚洲综合影院| 成熟少妇高潮喷水视频| 99国产精品一区二区蜜桃av| 老汉色av国产亚洲站长工具| 五月玫瑰六月丁香| 香蕉av资源在线| 精品人妻1区二区| 精品不卡国产一区二区三区| 欧美黄色片欧美黄色片| 18禁美女被吸乳视频| 我的老师免费观看完整版| 精品国产乱子伦一区二区三区| 两性夫妻黄色片| 中文字幕人妻丝袜一区二区| 老司机深夜福利视频在线观看| 亚洲欧美激情综合另类| 天堂动漫精品| 757午夜福利合集在线观看| 午夜精品久久久久久毛片777| 亚洲七黄色美女视频| 香蕉丝袜av| 在线观看美女被高潮喷水网站 | 午夜福利成人在线免费观看| 熟女少妇亚洲综合色aaa.| 欧美极品一区二区三区四区| 一卡2卡三卡四卡精品乱码亚洲| 免费看日本二区| 老司机福利观看| 精品久久久久久成人av| 亚洲第一电影网av| 亚洲人成网站高清观看| 中文在线观看免费www的网站| 国产激情久久老熟女| 欧美国产日韩亚洲一区| 国产精品 欧美亚洲| 欧洲精品卡2卡3卡4卡5卡区| www.熟女人妻精品国产| 亚洲天堂国产精品一区在线| 黄色丝袜av网址大全| 国产69精品久久久久777片 | 十八禁人妻一区二区| 99riav亚洲国产免费| 国产成人精品无人区| 天天添夜夜摸| 日本 av在线| 他把我摸到了高潮在线观看| 亚洲av成人av| 亚洲熟妇熟女久久| 色在线成人网| 在线观看一区二区三区| 97人妻精品一区二区三区麻豆| 亚洲乱码一区二区免费版| 一区福利在线观看| 日韩三级视频一区二区三区| 欧美zozozo另类| 免费在线观看亚洲国产| 国产97色在线日韩免费| 99热这里只有精品一区 | 好男人电影高清在线观看| 国产野战对白在线观看| 国产真人三级小视频在线观看| 美女黄网站色视频| 亚洲精品美女久久久久99蜜臀| 久久久国产欧美日韩av| 日本与韩国留学比较| 亚洲自拍偷在线| 午夜免费成人在线视频| 一a级毛片在线观看| 日韩欧美一区二区三区在线观看| 中文字幕人妻丝袜一区二区| 极品教师在线免费播放| 夜夜爽天天搞| 毛片女人毛片| 女警被强在线播放| 女人被狂操c到高潮| 老司机福利观看| 午夜激情福利司机影院| 99久久综合精品五月天人人| 色尼玛亚洲综合影院| 久久久精品大字幕| 少妇的逼水好多| 久久精品国产综合久久久| 国产伦在线观看视频一区| 亚洲自偷自拍图片 自拍| 成年版毛片免费区| 亚洲精品在线美女| av在线蜜桃| 久久久国产欧美日韩av| 不卡av一区二区三区| 国产成+人综合+亚洲专区| 精品福利观看| 99视频精品全部免费 在线 | 国产不卡一卡二| 欧美黑人欧美精品刺激| 亚洲av电影在线进入| 91av网站免费观看| 午夜免费观看网址| 国产99白浆流出| 18禁观看日本| 91九色精品人成在线观看| 欧美乱码精品一区二区三区| 欧美中文日本在线观看视频| 成人鲁丝片一二三区免费| 欧美乱妇无乱码| 美女高潮的动态| 人妻久久中文字幕网| 性欧美人与动物交配| xxx96com| 国产激情欧美一区二区| 久久热在线av| 亚洲午夜理论影院| 亚洲精品一卡2卡三卡4卡5卡| 国产精品国产高清国产av| 婷婷精品国产亚洲av| 亚洲专区国产一区二区| 九九在线视频观看精品| 亚洲一区二区三区色噜噜| 老熟妇仑乱视频hdxx| h日本视频在线播放| aaaaa片日本免费| 久久香蕉国产精品| 久久欧美精品欧美久久欧美| 91麻豆av在线| 黄色片一级片一级黄色片| 麻豆国产av国片精品| 成年女人毛片免费观看观看9| 国产亚洲精品综合一区在线观看| 亚洲色图av天堂| 人人妻,人人澡人人爽秒播| 美女黄网站色视频| 精品国产乱码久久久久久男人| 麻豆av在线久日| 少妇的丰满在线观看| 一级毛片女人18水好多| 日韩精品中文字幕看吧| 99国产极品粉嫩在线观看| 九九久久精品国产亚洲av麻豆 | av视频在线观看入口| 国产成人一区二区三区免费视频网站| 精品久久久久久久久久久久久| 久久精品亚洲精品国产色婷小说| 三级国产精品欧美在线观看 | 男人和女人高潮做爰伦理| 国产v大片淫在线免费观看| 成人一区二区视频在线观看| 精品福利观看| 神马国产精品三级电影在线观看| 国产精品免费一区二区三区在线| 狠狠狠狠99中文字幕| 久久久久九九精品影院| 亚洲人与动物交配视频| 国产单亲对白刺激| 91老司机精品| 一本综合久久免费| 一本综合久久免费| 国产精品一区二区精品视频观看| 欧美+亚洲+日韩+国产| 久久精品国产亚洲av香蕉五月| 麻豆国产97在线/欧美| 男插女下体视频免费在线播放| 狂野欧美激情性xxxx| 国产探花在线观看一区二区| 国产精品1区2区在线观看.| 日韩免费av在线播放| 国产精品香港三级国产av潘金莲| 国产成年人精品一区二区| 少妇裸体淫交视频免费看高清| 在线观看午夜福利视频| 国产v大片淫在线免费观看| 99热精品在线国产| 亚洲国产欧美网| 12—13女人毛片做爰片一| 亚洲美女黄片视频| 九九久久精品国产亚洲av麻豆 | 伊人久久大香线蕉亚洲五| 日韩成人在线观看一区二区三区| 美女高潮喷水抽搐中文字幕| 日本 欧美在线| 老鸭窝网址在线观看| 亚洲第一电影网av| 天天一区二区日本电影三级| 男女下面进入的视频免费午夜| 色综合站精品国产| 免费大片18禁| 免费观看的影片在线观看| 日本黄色片子视频| 真实男女啪啪啪动态图| 亚洲专区国产一区二区| 亚洲自拍偷在线| 色吧在线观看| 国产黄色小视频在线观看| 国产真实乱freesex| 国产成人精品久久二区二区免费| 我的老师免费观看完整版| 精品乱码久久久久久99久播| 免费人成视频x8x8入口观看| 欧美中文综合在线视频| 国产熟女xx| 91字幕亚洲| 久久久水蜜桃国产精品网| 91久久精品国产一区二区成人 | 他把我摸到了高潮在线观看| 亚洲av日韩精品久久久久久密| 免费观看的影片在线观看| 亚洲av中文字字幕乱码综合| 中文亚洲av片在线观看爽| 色综合亚洲欧美另类图片| 日本 欧美在线| 高潮久久久久久久久久久不卡| 国产亚洲av高清不卡| 国产高清激情床上av| 欧美成狂野欧美在线观看| 日韩欧美国产在线观看| 午夜视频精品福利| 麻豆国产97在线/欧美| 精品一区二区三区视频在线 | 在线观看免费午夜福利视频| av国产免费在线观看| 给我免费播放毛片高清在线观看| 91在线观看av| 男女做爰动态图高潮gif福利片| 亚洲成av人片免费观看| 啦啦啦韩国在线观看视频| 日韩人妻高清精品专区| 精品国产三级普通话版| 黄色视频,在线免费观看| 国产视频内射| 国产成+人综合+亚洲专区| 一级黄色大片毛片| 久久99热这里只有精品18| 日日干狠狠操夜夜爽| 国产亚洲精品久久久com| 欧美3d第一页| 真人做人爱边吃奶动态| 悠悠久久av| 他把我摸到了高潮在线观看| 久久天躁狠狠躁夜夜2o2o| 人妻久久中文字幕网| 亚洲自偷自拍图片 自拍| 国产av在哪里看| 国产精品爽爽va在线观看网站| 亚洲色图 男人天堂 中文字幕| 色尼玛亚洲综合影院| 国产精品一区二区三区四区免费观看 | 色精品久久人妻99蜜桃| 欧美国产日韩亚洲一区| 精品不卡国产一区二区三区| 丰满人妻一区二区三区视频av | 欧美三级亚洲精品| 国产黄a三级三级三级人| 亚洲一区高清亚洲精品| 人妻丰满熟妇av一区二区三区| 亚洲av成人不卡在线观看播放网| 亚洲性夜色夜夜综合| 在线免费观看不下载黄p国产 | 久久九九热精品免费| 国产综合懂色| 欧美乱码精品一区二区三区| 日本与韩国留学比较| 日本a在线网址| 亚洲 欧美一区二区三区| 亚洲美女视频黄频| 夜夜夜夜夜久久久久| 亚洲自偷自拍图片 自拍| 午夜亚洲福利在线播放| 亚洲国产欧美网| 美女黄网站色视频| 美女 人体艺术 gogo| 久久久久亚洲av毛片大全| 免费观看的影片在线观看| 成人国产综合亚洲| 一个人观看的视频www高清免费观看 | 国产成人精品久久二区二区免费| 男女床上黄色一级片免费看| 国产av一区在线观看免费| 午夜免费成人在线视频| 亚洲欧美精品综合久久99| 日本免费一区二区三区高清不卡| 男插女下体视频免费在线播放| 亚洲 欧美一区二区三区| 欧美一区二区国产精品久久精品| 亚洲人成网站高清观看| 在线永久观看黄色视频| 真实男女啪啪啪动态图| 久久精品aⅴ一区二区三区四区| 国产69精品久久久久777片 | 亚洲一区二区三区色噜噜| 观看美女的网站| 亚洲欧美日韩东京热| 啪啪无遮挡十八禁网站| 亚洲av五月六月丁香网| 这个男人来自地球电影免费观看| 国内精品美女久久久久久| 欧美成狂野欧美在线观看| 无限看片的www在线观看| 久久精品国产清高在天天线| 欧美3d第一页| 一区二区三区高清视频在线| 51午夜福利影视在线观看| 成人欧美大片| 亚洲精品国产精品久久久不卡| 真人做人爱边吃奶动态| 欧美绝顶高潮抽搐喷水| 成人一区二区视频在线观看| 久久精品国产亚洲av香蕉五月| 首页视频小说图片口味搜索| 亚洲精品在线观看二区| 成人av在线播放网站| 一级黄色大片毛片| 亚洲成av人片免费观看| 欧美又色又爽又黄视频| 午夜福利成人在线免费观看| 亚洲av美国av| 99国产精品99久久久久| 制服丝袜大香蕉在线| 国产精品国产高清国产av| 这个男人来自地球电影免费观看| 老熟妇仑乱视频hdxx| 两人在一起打扑克的视频| 巨乳人妻的诱惑在线观看| 欧美丝袜亚洲另类 | 久久久成人免费电影| av黄色大香蕉| 亚洲av成人av| 国产伦在线观看视频一区| cao死你这个sao货| 午夜两性在线视频| 欧美色视频一区免费| 男人舔奶头视频| 久久草成人影院| 色综合亚洲欧美另类图片| 精品国产乱子伦一区二区三区| 午夜福利免费观看在线| 日本熟妇午夜| 欧美成人一区二区免费高清观看 | 桃色一区二区三区在线观看| 国内揄拍国产精品人妻在线| 伦理电影免费视频| 丰满人妻一区二区三区视频av | 久久伊人香网站| 亚洲专区中文字幕在线| 夜夜夜夜夜久久久久| 亚洲片人在线观看| 精品午夜福利视频在线观看一区| 久久精品夜夜夜夜夜久久蜜豆| 国产在线精品亚洲第一网站| 成年女人永久免费观看视频| 夜夜躁狠狠躁天天躁| 天天添夜夜摸| 日日摸夜夜添夜夜添小说| 美女黄网站色视频| 高清在线国产一区| 欧美精品啪啪一区二区三区| 狠狠狠狠99中文字幕| 特级一级黄色大片| 岛国在线免费视频观看| 亚洲中文字幕日韩| 美女cb高潮喷水在线观看 | 99国产精品一区二区蜜桃av| 国产乱人伦免费视频| 亚洲国产欧美一区二区综合| 国产爱豆传媒在线观看| 久9热在线精品视频| 特大巨黑吊av在线直播| 精品久久久久久成人av| av天堂中文字幕网| 动漫黄色视频在线观看| 岛国在线免费视频观看| 国产伦一二天堂av在线观看| 国产av在哪里看| 久久人人精品亚洲av| 69av精品久久久久久| 久久伊人香网站| 日韩欧美 国产精品| 久久久水蜜桃国产精品网| 人妻久久中文字幕网| 国产午夜精品论理片| 久久婷婷人人爽人人干人人爱| 特大巨黑吊av在线直播| 亚洲一区二区三区色噜噜| 一个人免费在线观看电影 | 亚洲专区国产一区二区| 偷拍熟女少妇极品色| 欧洲精品卡2卡3卡4卡5卡区| 我的老师免费观看完整版| 91在线精品国自产拍蜜月 | 欧美日韩精品网址| 国产精品久久久久久人妻精品电影| 三级国产精品欧美在线观看 | 精品不卡国产一区二区三区| 三级毛片av免费| 色哟哟哟哟哟哟| 久久久久九九精品影院| 亚洲乱码一区二区免费版| 成人鲁丝片一二三区免费| 少妇丰满av| 国产探花在线观看一区二区| 国产亚洲欧美在线一区二区| 国产 一区 欧美 日韩| 久久久精品欧美日韩精品| 色综合婷婷激情| 日本 欧美在线| 国产高清三级在线| 精品久久久久久,| 99在线视频只有这里精品首页| 国产亚洲精品一区二区www| 久久中文看片网| 青草久久国产| 久久精品91蜜桃| 国产成人精品无人区| 天堂av国产一区二区熟女人妻| 久久久成人免费电影| 熟女少妇亚洲综合色aaa.| 免费电影在线观看免费观看| 国产成人av激情在线播放| 亚洲熟妇熟女久久| 亚洲中文日韩欧美视频| 国模一区二区三区四区视频 | 最近最新免费中文字幕在线| 夜夜看夜夜爽夜夜摸| 久久九九热精品免费| 亚洲国产精品sss在线观看| 久久草成人影院| 亚洲片人在线观看| 色尼玛亚洲综合影院| 深夜精品福利| 国产激情久久老熟女| 此物有八面人人有两片| 女人高潮潮喷娇喘18禁视频| 变态另类丝袜制服| 国产精品 欧美亚洲| 亚洲av成人不卡在线观看播放网| 成人午夜高清在线视频| 99国产精品一区二区三区| 亚洲国产欧美人成| 99热这里只有是精品50| 五月玫瑰六月丁香| 99国产精品99久久久久| 亚洲国产欧洲综合997久久,| 天堂网av新在线| 色综合站精品国产| 淫妇啪啪啪对白视频| 久久久久久久午夜电影| 亚洲一区二区三区不卡视频| 午夜福利免费观看在线| 国产免费男女视频| 国产欧美日韩精品亚洲av| 99国产精品99久久久久| 欧美黄色片欧美黄色片| 亚洲成av人片在线播放无| 免费观看人在逋| 熟妇人妻久久中文字幕3abv| 岛国在线免费视频观看| 国产欧美日韩精品一区二区| 噜噜噜噜噜久久久久久91| 国产精品自产拍在线观看55亚洲| 国产av一区在线观看免费| 色哟哟哟哟哟哟| 国产美女午夜福利| 国产午夜福利久久久久久| 亚洲成a人片在线一区二区| 日韩人妻高清精品专区| 成熟少妇高潮喷水视频| 男女下面进入的视频免费午夜| 黑人欧美特级aaaaaa片| 宅男免费午夜| 亚洲国产欧美网| 亚洲av电影在线进入| 欧美不卡视频在线免费观看| 亚洲精品中文字幕一二三四区| 国产精品一区二区精品视频观看| 黄频高清免费视频| 精品不卡国产一区二区三区| 18禁美女被吸乳视频| 亚洲人成网站在线播放欧美日韩| 亚洲精品久久国产高清桃花| 伊人久久大香线蕉亚洲五| 亚洲人成伊人成综合网2020| 久久这里只有精品中国| 国产成人啪精品午夜网站| 草草在线视频免费看| 亚洲国产精品久久男人天堂| 亚洲熟妇熟女久久| www.熟女人妻精品国产| 综合色av麻豆| 国产精品1区2区在线观看.| 18美女黄网站色大片免费观看| av在线蜜桃| 又大又爽又粗| 黑人巨大精品欧美一区二区mp4|