基于XGBoost的内陆河湖浊度反演与长时序分析

李媛媛, 沈芳, 陈嵩钰, 魏小岛. 基于XGBoost的内陆河湖浊度反演与长时序分析[J]. 环境工程学报, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
引用本文: 李媛媛, 沈芳, 陈嵩钰, 魏小岛. 基于XGBoost的内陆河湖浊度反演与长时序分析[J]. 环境工程学报, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
LI Yuanyuan, SHEN Fang, CHEN Songyu, WEI Xiaodao. Inversion and long-term series analysis of turbidity in inland rivers and lakes based on XGBoost[J]. Chinese Journal of Environmental Engineering, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
Citation: LI Yuanyuan, SHEN Fang, CHEN Songyu, WEI Xiaodao. Inversion and long-term series analysis of turbidity in inland rivers and lakes based on XGBoost[J]. Chinese Journal of Environmental Engineering, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125

基于XGBoost的内陆河湖浊度反演与长时序分析

    作者简介: 李媛媛 (1998—) ,女,硕士研究生,51213904006@stu.ecnu.edu.cn
    通讯作者: 沈芳(1964—),女,博士,教授,fshen@sklec.ecnu.edu.cn
  • 基金项目:
    上海市科委重点项目 (20dz1204700) ;中国长江三峡集团有限公司科研项目资助 (202103552)
  • 中图分类号: X832

Inversion and long-term series analysis of turbidity in inland rivers and lakes based on XGBoost

    Corresponding author: SHEN Fang, fshen@sklec.ecnu.edu.cn
  • 摘要: 浊度是影响水下光场及营养盐循环的关键要素之一,浊度监测可以为河湖水质的污染防控和预警提供科学依据。以长三角示范区的典型河湖为研究区,使用实测数据构建浊度反演模型,并利用1984—2022年Landsat卫星数据分析了研究区河湖浊度的长时序动态变化。对比传统经验模型、半经验模型和机器学习模型,XGBoost机器学习模型精度最高 (R2为0.68,RMSE为4.78 NTU) 。浊度反演结果表明,近40年河流航道和淀山湖北部非渔场区域浊度上升了10%和12%,而元荡湖和大莲湖浊度下降了19%和27%,并且浊度随着建设用地面积的增加而增大;研究区浊度季节性变化显著,秋冬季平均浊度比春夏季高6 NTU,月平均浊度与月平均降水量负相关,相关系数为-0.61 (p<0.05) ,但与月平均风速没有显著的相关性。基于XGBoost的Landsat长时序浊度反演能够把握研究区浊度的时空变化趋势,明确水污染管理与治理方向,最终助力长三角一体化发展。
  • 硝酸盐和高氯酸盐是地下水中常见的共存污染物。农业氮肥的大量施用是地下水硝酸盐污染的最大来源[1]。在美国,约有22%的农业地下水硝酸盐浓度超过饮用水标准[2];在欧洲,有1/3的地下水硝酸盐超过了饮用水标准,硝酸盐达到了100~150 mg·L−1[3];在中国,有61.2%的地下水井存在着硝酸盐严重超标的问题[4]。当长期饮用硝酸盐超标的地下水时,人体患高铁血红蛋白症和癌症的几率会增加[5]。地下水中高氯酸盐主要来源于火箭助推剂、烟花、染料和油漆等工业产品的生产和加工过程[6-7]。有研究表明,智利阿塔卡马沙漠中地表水中的高氯酸盐质量浓度为744~1 480 μg·L−1[8]。有研究[9]表明,由于将含有高氯酸盐的制造业废水排到地下水中,导致位于美国加利福尼亚州东部的萨克拉门地下水高氯酸盐含量达到了8 mg·L−1。KOSAKA等[10]在2007年对日本Usui河和Tone河调研发现,受污域河水中高氯酸盐质量浓度为0.34~2.38 mg·L−1。据调查,我国部分水体存在高浓度高氯酸盐污染状况,如湖南省衡阳市某地表水高氯酸盐含量高达6.8~54.4 mg·L−1[11]。由于高氯酸根离子在电荷和离子半径上都与碘相似,因此,可能会破坏甲状腺中碘的吸收,从而可能影响甲状腺功能,引起甲状腺疾病[12-14],严重威胁着人类的健康。

    目前处理水体硝酸盐和高氯酸盐复合污染的方法主要包括离子交换、膜分离、化学还原和生物法[15-16]。离子交换和膜分离工艺虽然能够去除水中的硝酸盐和高氯酸盐,但只是将污染物浓缩转移,并没有将污染物无害转化,且操作成本较高;对于化学还原法,多采用贵金属催化剂,价格昂贵且催化剂易失活;生物还原工艺能够实现硝酸盐和高氯酸盐的高效去除,易于推广应用,故逐渐成为本领域的研究热点[17-18]。根据所需电子供体的不同,生物还原法分为异养还原和自养还原法,对于异养还原过程,常用的有机碳源包括葡萄糖、甲醇、乙酸盐等[19-20]。采用乙酸盐作为电子供体,还原去除硝酸盐和高氯酸盐的反应分别如式(1)和式(2)所示。

    7NO3+13CH3COO+H2O3C5H7O2N+2N2+11CO2+20OH (1)
    3ClO4+2CH3COO+2H2OCl+4CO2+10OH (2)

    需要指出的是,虽然异养生物还原反应速率快,但其反应过程会产生碱度,造成出水pH的升高。此外,碳源投加不易控制,投加过少则处理不彻底,投加过多则残余水中容易造成二次污染。生物自养过程利用无机碳作为电子供体,避免了有机源的投加,近年来受到研究人员的关注。其中,微生物利用单质硫作为电子供体,还原硝酸盐和高氯酸盐的反应如式(3)和式(4)所示。对于硫自养生物还原过程来说,虽然避免了有机碳源投加带来的隐患,但反应过程中产生副产物硫酸盐,同时消耗水体碱度,造成出水pH降低。

    1.06NO3+1.11S+0.3CO2+0.785H2O0.06C5H7O2N+0.5N2+1.11SO24+1.16H+ (3)
    3ClO4+4S0+4H2O3Cl+4SO24+8H+ (4)

    基于上述研究背景,为了克服异养和自养还原各自的缺点,本研究建立了异养和硫自养协同作用的一体式生物反应器,其中有机碳源投加量低于理论值,这可有效避免水中有机物的二次污染;硫自养仅承担部分负荷,从而削弱硫酸盐的产生;同时平衡了2种反应过程碱度产生与消耗,稳定出水pH;考察了在不同HRT和不同碳源投加量条件下,混合营养生物工艺对废水中硝酸盐和高氯酸盐的去除效果,以期为该工艺在实际水处理工程中的应用提供参考。

    建立了混合营养一体式微生物反应器(图1)。反应器采用有机玻璃柱制成,柱内径为5 cm,柱内填充硫磺颗粒(10目筛子筛分,燕山石化),填充高度42 cm,反应器有效体积310 mL,外部有水浴保温夹层,温度为(30±2) ℃。反应器顶部采用溢流的方式出水,底部通过蠕动泵(YZ1513,保定兰格)均匀的将进水送进反应器。

    图 1  生物反应器实验装置图
    Figure 1.  Diagram of the experimental bioreactor

    反应器进水采用去离子水和自来水(1∶4)配制而成,加入优级纯高氯酸钠和硝酸钠作为目标污染物。进水组成:ClO4质量浓度(20.12±0.12) mg·L−1NO3-N质量浓度(20.21±0.23) mg·L−1,CH3COO质量浓度60 mg·L−1(第I~IV阶段添加,第Ⅴ阶段无添加),KH2PO4 0.5 mg·L−1(以磷计),缓冲物质NaHCO3 300 mg·L−1。另外,为了给微生物提供足够的营养成分,进水加入1 mL·L−1的微量元素溶液,溶液成分为0.05 g·L−1 ZnCl2·H2O、0.11 g·L−1 NiCl2·6H2O、0.49 g·L−1 MnCl2·4H2O、0.05 g·L−1 H3BO3、2.00 g·L−1 CoCl2·6H2O、0.11 g·L−1 CuCl2·2H2O。其余进水水质参数为 72.26~138.52 mg·L−1 SO24、380.23~497.26 mg·L−1碱度、pH为8.23~8.76。

    取郑州市五龙口污水处理厂厌氧段污泥,将污泥与硫磺颗粒均匀的填充到反应器中,污泥接种量为6.5 g·L−1。驯化启动阶段进水ClO4质量浓度为40 mg·L−1NO3-N质量浓度为40 mg·L−1,CH3COO质量浓度为180 mg·L−1,磷源KH2PO4 为0.5 mg·L−1(以磷计),缓冲物质NaHCO3 500 mg·L−1。反应器在HRT为8 h的情况下连续运行14 d后完成驯化。

    反应器驯化启动完成后,开始正式运行,运行方案如表1所示。实验分为5个阶段,前4个阶段为混合营养工况,调整HRT分别为4、2、1和0.5 h,每个阶段运行10 d。根据反应式(1)和式(2),理论上还原20 mg·L−1NO3-N和ClO4需要164 mg·L−1的CH3COO。在前4个阶段,进水投加CH3COO的质量浓度为60 mg·L−1,以营造低碳源投加的混合营养条件。第Ⅴ阶段,保持HRT为0.5 h,取消有机碳源投加,为单独硫自养工艺,研究反应器的处理效能,并与混合营养工艺作对比。

    表 1  反应器实验运行方案
    Table 1.  Experimental operation scheme of reactor
    运行阶段运行时段/d工艺条件HRT/h进水CH3COO/(mg·L−1)
    1~10混合营养460
    11~20混合营养260
    21~30混合营养160
    31~41混合营养0.560
    42~51单独硫自养0.50
     | Show Table
    DownLoad: CSV

    所有出水水样均经过0.22 μm滤膜过滤,并在4 ℃条件下保存,采用离子色谱仪(ICS-2100,美国赛默飞)测定水样中的ClO4(最低检测限为0.2 mg·L−1)和SO24(最低检测限为0.4 mg·L−1)。用紫外可见分光光度法测NO3,N-(1-萘基)-乙二胺光度法测NO2,利用TOC分析仪(TOCL-CPN CN200,日本岛津)来测定进出水NPOC,使用pH计(PHS-3C,雷磁,中国)测定进出水pH,采用酸碱滴定法测定进出水的碱度。

    混合营养生物反应器启动期运行效果如图2所示。在HRT为8 h的条件下,反应器只需要2 d的适应期就可以将40 mg·L−1NO3-N还原,出水NO3-N质量浓度为(0.25±0.06) mg·L−1,并且驯化期间没有出现NO2的积累,出水浓度始终低于检出限;而对于ClO4,生物反应器则需要8 d的适应期才能够逐渐将40 mg·L−1ClO4还原,出水ClO4浓度低于检出限。产生这种现象的原因可能是,接种污泥取自郑州市五龙口污水处理厂的厌氧段污泥,污泥本身就具有反硝化功能,因此,对NO3具有较强的适应能力,而对于ClO4,需要一定的适应时间。至第5天,反应器出水SO24质量浓度增加量由116.8 mg·L−1升高至326.7 mg·L−1,这表明硫自养菌正在逐渐驯化成熟,至驯化末期,反应器出水SO24质量浓度增加量稳定在(292.4±6.3) mg·L−1

    图 2  驯化阶段ClO4NO3-N、NO2-N 和SO24质量浓度变化
    Figure 2.  Variations of ClO4, NO3-N, NO2-N and SO24 concentrations during acclimation process

    各运行阶段,反应器对ClO4NO3的去除效果如图3所示,在整个运行期间,反应器对ClO4NO3均表现出较好的去除效果。对于ClO4来说,当HRT从4 h缩短至2 h时,经过3 d的适应期,出水去除率趋于稳定,出水ClO4质量浓度低于检出限,当HRT降低到1.0 h和0.5 h时,出水ClO4质量浓度分别降低到(0.34±0.06) mg·L−1和(0.35±0.09) mg·L−1。第Ⅴ阶段,单独硫自养工艺下,反应器出水ClO4质量浓度升高至(0.58±0.12) mg·L−1,对应去除率为96.2%,相较于混合营养工艺,ClO4去除率降低了2%,表明混合营养工艺对ClO4的去除效果优于单独硫自养。LV等[21]用甲烷膜生物反应器还原ClO4NO3时发现:当进水中只有ClO4时,ClO4还原速率为1.7 mmol·(m2·d)−1;而当NO3存在时,ClO4还原速率降为0.64 mmol·(m2·d)−1,当NO3被完全去除后,ClO4还原速率又恢复到1.6 mmol·(m2·d)−1。DENIZ等[22]用自养-异养组合工艺去除ClO4NO3-N时发现:NO3-N在被完全还原之后ClO4才能被完全还原。在本研究中,混合营养工艺条件下,反应器实现了对2种污染物的高效去除,加入60 mg·L−1 CH3COO首先促进了NO3的还原效率,而后提高了ClO4的去除率,使得混合营养生物反应器在0.5 h时对ClO4的去除效果优于硫自养生物反应器。

    图 3  不同操作条件下进出水ClO4NO3-N和NO2-N质量浓度变化
    Figure 3.  Variations of ClO4, NO3-N and NO2-N concentrations in influent and effluent under different operate conditions

    反应器运行期间始终保持对NO3较高的去除效率,几乎不受HRT的影响,出水NO3-N质量浓度小于0.62 mg·L−1,去除率大于96.9%,实验期间不存在NO2-N的积累,出水NO2-N度低于检出限。停止添加碳源,反应器对NO3-N的去除几乎没有受到影响。

    围绕2种污染物异养和硫自养过程,本课题组已经开展了部分研究。2017年,本课题组[23]在硫自养固定床生物反应器还原ClO4NO3复合污染过程的研究中发现,2种污染物的还原同时发生,且硝酸盐的还原速率快于高氯酸盐,硫歧化反应伴随着高氯酸盐的还原,共存的NO3对硫歧化反应有抑制作用。2018年,本课题组[24]采用乙酸钠为碳源还原水中的高氯酸盐,发现以高氯酸盐为电子受体,长期驯化得到的混合菌群,对水中的硝酸盐依然有着较好的反硝化效果,在碳源有限的条件下,2种污染物的还原同时发生,且硝酸盐的还原率高于高氯酸盐。一般来说,生物过程对硝酸盐的去除优先于高氯酸盐。在本研究中,在混合营养条件下(HRT为4、2、1和0.5 h),反应器稳态运行时对硝酸盐的去除率高于对高氯酸盐的去除率。当停止碳源投加时,对高氯酸盐的去除受到一定程度的影响,去除率降低了2%,而对硝酸盐的去除没有受到影响,与DENIZ等报道的硝酸盐的去除优先于高氯酸盐的结果相一致[22]

    需要指出的是,第Ⅳ阶段,HRT仅为0.5 h,本反应器对NO3-N去除负荷达到930.24 g·(m3·d)−1,对 ClO4的去除负荷达到了943.20 g·(m3·d)−1。WAN等[25]建立了硫自养填充床生物反应器,其对ClO4NO3的最大去除负荷分别为705.92 g·(m3·d)−1和721.60 g·(m3·d)−1,去除负荷低于本实验。ZHU等[26]用生物活性炭填充床生物反应器,发现实验期间反应器出水高氯酸盐质量浓度为1.82~2.16 mg·L−1,出水高氯酸盐质量浓度远高于混合营养生物反应器。因此,本研究建立的混合营养生物反应器对ClO4NO3具有较好的去除效果。

    反应器进出水SO24质量浓度变化如图4所示。在第Ⅰ~Ⅳ阶段,混合营养工艺条件下,出水SO24质量浓度随着HRT的缩短而减少。当HRT为4 h时,出水SO24质量浓度增加量为(273±10) mg·L−1,随着HRT降低到2、1和0.5 h,出水SO24质量浓度增加量分别为(217±11)、(147±20)和(129±3) mg·L−1。第Ⅴ阶段,单独硫自养工艺条件下,出水SO24质量浓度增加量为(192±4) mg·L−1,同有混合营养工艺相比,SO24的产量增加63 mg·L−1

    图 4  不同操作条件下进出水SO24浓度变化和SO24浓度增加量
    Figure 4.  Variations of SO24 concentration in influent and effluent and its increase amount under different operate conditions

    WAN等[27]利用氢自养还原水中高氯酸盐,探究了共存的硝酸盐和硫酸盐对高氯酸盐去除的影响。实验结果表明,还原菌群对3种电子受体的还原速率顺序为NO3>ClO4>SO24。此外,有研究[28]表明,无论从热力学或动力学方面分析,硝酸盐还原均优先于硫酸盐。理论上,还原20 mg·L−1NO3-N,需要156 mg·L−1的CH3COO。在本研究中,进水投加CH3COO的质量浓度仅为60 mg·L−1,投加碳源不足,硝酸盐将优先于硫酸盐还原,故可推测,硫酸盐直接被有机碳源还原的可能性较低。

    根据式(3)和式(4),硫自养还原20 mg·L−1NO3-N和ClO4共计产生169 mg·L−1SO24。因此,当HRT为1 h和0.5 h时,投加60 mg·L−1 CH3COO能够有效地减少SO24产生,SO24的产量分别减少了22 mg·L−1和40 mg·L−1,HRT越短,SO24的生成受到的抑制越显著。当HRT为4 h和2 h时,虽然投加了60 mg·L−1 CH3COO,但SO24产生量依然大于硫自养的理论值,由此可推测,在较长的HRT下,发生了硫歧化反应。

    JU等[29]报道,在硫自养还原ClO4过程中,出水SO24质量浓度远大于理论值,据此推测出有硫歧化反应发生。WAN等[30]发现,在硫自养还原ClO4NO3复合污染过程中,硝酸盐的还原速率快于高氯酸盐,硫歧化反应伴随着高氯酸盐的还原,共存的NO3对该反应有抑制作用,HRT越长,该反应进行越剧烈,反应如式(5)所示。

    4S0+4H2O3H2S+SO24+2H+ (5)

    在混合营养条件下,当HRT为4 h和2 h时,SO24质量浓度增加量高于硫自养理论值(104±10) mg·L−1 和(48±10) mg·L−1,证实有硫歧化反应发生。因此,为削弱副产物SO24的产生,在保证NO3ClO4去除率的同时,工程上应当尽可能地采用较短的HRT。

    反应器进出水pH和碱度的变化如图5所示。在运行期间,进水pH为8.51~8.72,出水pH稳定在7.18~7.68。实验结果表明,出水碱度消耗量随着HRT的降低呈现降低的趋势。在混合营养工艺下,当HRT为4 h时,碱度消耗量为(213±13) mg·L−1(以CaCO3计);当HRT缩短至0.5 h时,碱度消耗量降为(105±7) mg·L−1,碱度消耗削减了49%。结合式(3)和式(4),理论上硫自养还原20 mg·L−1ClO4NO3-N的碱度消耗为105 mg·L−1(以CaCO3计)。对比第IV和第Ⅴ阶段,单独硫自养出水pH低于混合营养工艺,碱度消耗量为(150±2) mg·L−1,比混合营养工艺的碱度消耗量增加了45 mg·L−1。这是因为硫自养承担所有负荷,导致碱度消耗增加。在混合营养条件下,异养还原硝酸盐和高氯酸过程产生碱度,硫自养还原硝酸盐和高氯酸盐过程消耗碱度,此外,硫歧化反应也消耗碱度。整个实验运行期间,异养产生的碱度能够抵消一部分自养消耗的碱度,自养-异养在加快去除效率的同时,还能够降低反应物的碱度消耗和硫酸根浓度的产生。在HRT为4 h的条件下,混合营养工艺碱度消耗高于硫自养理论值(108±6) mg·L−1,表明即使在混合营养条件下,较长的HRT也会发生较为剧烈的硫歧化反应,导致碱度消耗较高。

    图 5  不同操作条件下进出水pH和碱度消耗
    Figure 5.  pH and alkalinity consumption under different operate conditions

    有机碳源是否有残余是影响反应器出水水质的关键所在。本反应器在运行期间,进出水NPOC质量浓度变化如图6所示。在混合营养条件下,反应器出水的NPOC始终维持在较低的水平,出水NPOC小于(2.68±0.13) mg·L−1,且不受HRT的影响。完全去除20 mg·L−1NO3-N和ClO4需要164 mg·L−1的CH3COO。本研究中,进水仅加入60 mg·L−1CH3COO,有机碳源的投加远低于理论值,因此,反应器出水中无有机碳源残余。整个实验运行期间,进水COD为(60.56±0.23) mg·L−1,出水COD小于5.68 mg·L−1

    图 6  不同操作条件下进出水NPOC质量浓度变化
    Figure 6.  Variations of NPOC in influent and effluent under different operate conditions

    1)混合营养一体式生物反应器能够同步高效的去除ClO4NO3-N。对于20 mg·L−1的进水ClO4NO3-N,在HRT分别为4 h和2 h时,出水ClO4质量浓度均低于检出限,当HRT为1 h和0.5 h时,出水ClO4 <0.34 mg·L−1;HRT对NO3-N的去除影响较小,出水NO3-N<0.62 mg·L−1

    2)在混合营养条件下,当HRT为0.5 h、温度为(30±2) ℃ 时,反应器中ClO4NO3的去除负荷达到最大,分别为943.20 g·(m3·d)−1和930.24 g·(m3·d)−1

    3)混合营养工艺能够有效降低SO24的产生量和碱度消耗量。当HRT为0.5 h时,与单独硫自养相比,SO24的产生减少了63 mg·L−1,碱度消耗减少了45 mg·L−1

    4)由于异养施加不足量碳源,混合营养工艺能有效避免有机碳源残余水中造成二次污染的问题,出水NPOC低于(2.68±0.13) mg·L−1

  • 图 1  研究区及采样点分布

    Figure 1.  Study area and spatial distribution of sampling sites

    图 2  XGBoost模型结构

    Figure 2.  Model structure of XGBoost

    图 3  浊度的光谱响应

    Figure 3.  Spectral response of turbidity

    图 4  浊度与总悬浮物浓度的相关性

    Figure 4.  Correlation between turbidity and TSM

    图 5  实测等效遥感反射率和大气校正遥感反射率的相关性

    Figure 5.  Correlation between measured equivalent Rrs and atmospherically-corrected Rrs

    图 6  实测浊度与预测浊度散点图

    Figure 6.  Scatter plots of measured turbidity and predicted turbidity

    图 7  Landsat TM/ETM+/OLI反演的1984‒2022年研究区浊度分布

    Figure 7.  Distributions of turbidity in the study area for 1984‒2022 retrieved by Landsat TM/ETM+/OLI

    图 8  Landsat TM/ETM+/OLI反演的1984‒2022年研究区浊度变化率

    Figure 8.  Variation of turbidity in the study area for 1984‒2022 retrieved by Landsat TM/ETM+/OLI

    图 9  目标区域浊度的年际变化

    Figure 9.  Inter-annual changes of turbidity in the target area

    图 10  浊度和建设用地面积的相关性

    Figure 10.  Correlation between turbidity and built-up land area

    图 11  Landsat TM/ETM+/OLI反演的1984‒2022年研究区浊度季节性变化

    Figure 11.  Seasonal variation of turbidity in the study area for 1984‒2022 retrieved by Landsat TM/ETM+/OLI

    图 12  研究区1984‒2022年月平均浊度与降水量、风速的关系

    Figure 12.  Relationship between monthly mean turbidity, precipitation, and wind speed in the study area for 1984‒2022

    表 1  采样位置、时间和数量

    Table 1.  Sampling location, time and number of samples

    编号 位置 时间 采样数量 有效数据量
    1 太湖西山岛以北 2022-10-24 35 6
    2 元荡湖 2022-10-25 35 29
    3 太浦河、汾湖、三白荡 2022-11-23 18 15
    4 淀山湖南部 2022-12-13 20 15
    5 淀山湖、元荡湖 2023-03-26 18 16
    编号 位置 时间 采样数量 有效数据量
    1 太湖西山岛以北 2022-10-24 35 6
    2 元荡湖 2022-10-25 35 29
    3 太浦河、汾湖、三白荡 2022-11-23 18 15
    4 淀山湖南部 2022-12-13 20 15
    5 淀山湖、元荡湖 2023-03-26 18 16
    下载: 导出CSV

    表 2  基于QAA计算OLI的颗粒物后向散射系数bbp(655)

    Table 2.  Calculation of bbp at 655 nm for OLI based on QAA

    步骤 计算
    1 rrs(λ)=Rrs(λ)/(0.52+1.7Rrs(λ))
    2 μ(λ)=g0+(g0)2+4g1×rrs(λ)2g1, 其中g0= 0.089, g1= 0.1245
    如果Rrs(655)<0.0015 sr−1 否则655 nm =λ0
    3 χ=log(rrs(443)+rrs(483)rrs(561)+5rrs(655)/rrs(483)×rrs(655))a(λ0)=a(561)=aw(561)+10h0+h1χ+h2χ2h0=1.146,h1=1.366,h2=0.469 a(λ0)=a(655)=aw(655)+0.39(Rrs(655)Rrs(443)+Rrs(483))1.14
    4 bbp(λ0)=bbp(561)=μ(561)×a(561)1μ(561)bbw(561) bbp(λ0)=bbp(655)=μ(655)×a(655)1μ(655)bbw(655)
    5 η=2.0(11.2exp(0.9rrs(443)/rrs(561)))
    6 bbp(λ)=bbp(λ0)(λ0/λ)η
      注:表中aaw分别为总吸收系数和纯海水吸收系数,m−1bbwbbp分别为纯海水和悬浮颗粒的后向散射系数,m−1rrs(λ)为水表面以下遥感反射率,sr−1Rrs(λ)是水面以上遥感反射率,sr−1
    步骤 计算
    1 rrs(λ)=Rrs(λ)/(0.52+1.7Rrs(λ))
    2 μ(λ)=g0+(g0)2+4g1×rrs(λ)2g1, 其中g0= 0.089, g1= 0.1245
    如果Rrs(655)<0.0015 sr−1 否则655 nm =λ0
    3 χ=log(rrs(443)+rrs(483)rrs(561)+5rrs(655)/rrs(483)×rrs(655))a(λ0)=a(561)=aw(561)+10h0+h1χ+h2χ2h0=1.146,h1=1.366,h2=0.469 a(λ0)=a(655)=aw(655)+0.39(Rrs(655)Rrs(443)+Rrs(483))1.14
    4 bbp(λ0)=bbp(561)=μ(561)×a(561)1μ(561)bbw(561) bbp(λ0)=bbp(655)=μ(655)×a(655)1μ(655)bbw(655)
    5 η=2.0(11.2exp(0.9rrs(443)/rrs(561)))
    6 bbp(λ)=bbp(λ0)(λ0/λ)η
      注:表中aaw分别为总吸收系数和纯海水吸收系数,m−1bbwbbp分别为纯海水和悬浮颗粒的后向散射系数,m−1rrs(λ)为水表面以下遥感反射率,sr−1Rrs(λ)是水面以上遥感反射率,sr−1
    下载: 导出CSV

    表 3  遥感反射率光谱特征组合与浊度的相关性

    Table 3.  Correlation between combinations of Rrs spectral features and turbidity

    波段组合 相关系数r 波段组合 相关系数r
    B2 0.49** B3 0.50**
    B4 0.71** B5 0.65**
    B4/B2 0.78** B4/B3 0.88**
    B4/B5 −0.23* (B4-B2)/(B4+B2) 0.76**
    (B4-B3)/(B4+B3) 0.87** (B4-B5)/(B4+B5) −0.38**
      注:*表示相关性在0.05水平上显著;**表示相关性在0.01水平上显著。
    波段组合 相关系数r 波段组合 相关系数r
    B2 0.49** B3 0.50**
    B4 0.71** B5 0.65**
    B4/B2 0.78** B4/B3 0.88**
    B4/B5 −0.23* (B4-B2)/(B4+B2) 0.76**
    (B4-B3)/(B4+B3) 0.87** (B4-B5)/(B4+B5) −0.38**
      注:*表示相关性在0.05水平上显著;**表示相关性在0.01水平上显著。
    下载: 导出CSV

    表 4  浊度反演模型参数化及比较

    Table 4.  Parameterization and comparison of turbidity inversion models

    序号 模型 表达式 留一交叉验证
    R2 RMSE (NTU) MAPE (%)
    1 二次多项式 y=136.38x1297.36x1+19.55 0.78 5.25 31.43
    2 指数函数 y=5e2.72x115.16 0.78 5.25 30.35
    3 幂函数 y=62.7x12.993.28 0.78 5.26 30.93
    4 SESB-based y=117.14ρw(λ)/(1ρw(λ)/0.1686)+9.63 0.39 8.72 57.91
    5 QAA-based y=74.14x0.962 0.54 7.60 48.83
    6 XGBoost y=fXGBoost(x1) 0.83 4.58 23.79
      注:x1表示OLI传感器B4/B3x2表示bbp(655),m−1y表示实测浊度,NTU;fXGBoost表示XGBoost模型函数。
    序号 模型 表达式 留一交叉验证
    R2 RMSE (NTU) MAPE (%)
    1 二次多项式 y=136.38x1297.36x1+19.55 0.78 5.25 31.43
    2 指数函数 y=5e2.72x115.16 0.78 5.25 30.35
    3 幂函数 y=62.7x12.993.28 0.78 5.26 30.93
    4 SESB-based y=117.14ρw(λ)/(1ρw(λ)/0.1686)+9.63 0.39 8.72 57.91
    5 QAA-based y=74.14x0.962 0.54 7.60 48.83
    6 XGBoost y=fXGBoost(x1) 0.83 4.58 23.79
      注:x1表示OLI传感器B4/B3x2表示bbp(655),m−1y表示实测浊度,NTU;fXGBoost表示XGBoost模型函数。
    下载: 导出CSV
  • [1] PETUS C, CHUST G, GOHIN F, et al. Estimating turbidity and total suspended matter in the Adour River plume (South Bay of Biscay) using MODIS 250-m imagery[J]. Continental Shelf Research, 2010, 30(5): 379-392. doi: 10.1016/j.csr.2009.12.007
    [2] NECHAD B, RUDDICK K G, PARK Y. Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[J]. Remote Sensing of Environment, 2010, 114(4): 854-866. doi: 10.1016/j.rse.2009.11.022
    [3] DING W H, ZHAO J X, QIN B Q, et al. Exploring and quantifying the relationship between instantaneous wind speed and turbidity in a large shallow lake: case study of Lake Taihu in China[J]. Environmental Science and Pollution Research, 2021, 28(13): 16616-16632. doi: 10.1007/s11356-020-11544-y
    [4] LIN X N, WU M, SHAO X X, et al. Water turbidity dynamics using random forest in the Yangtze River Delta Region, China[J]. Science of the Total Environment, 2023, 903: 166511. doi: 10.1016/j.scitotenv.2023.166511
    [5] HOU X J, FENG L, DUAN H T, et al. Fifteen-year monitoring of the turbidity dynamics in large lakes and reservoirs in the middle and lower basin of the Yangtze River, China[J]. Remote Sensing of Environment, 2017, 190: 107-121. doi: 10.1016/j.rse.2016.12.006
    [6] YIN Z Y, LI J S, LIU Y, et al. Water clarity changes in Lake Taihu over 36 years based on Landsat TM and OLI observations[J]. International Journal of Applied Earth Observation and Geoinformation, 2021, 102: 102457. doi: 10.1016/j.jag.2021.102457
    [7] DOXARAN D, FROIDEFOND J, CASTAING P, et al. Dynamics of the turbidity maximum zone in a macrotidal estuary (the Gironde, France): Observations from field and MODIS satellite data[J]. Estuarine, Coastal and Shelf Science, 2009, 81(3): 321-332. doi: 10.1016/j.ecss.2008.11.013
    [8] DOGLIOTTI A I, RUDDICK K G, NECHAD B, et al. A single algorithm to retrieve turbidity from remotely-sensed data in all coastal and estuarine waters[J]. Remote Sensing of Environment, 2015, 156: 157-168. doi: 10.1016/j.rse.2014.09.020
    [9] MA Y, SONG K S, WEN Z D, et al. Remote sensing of turbidity for lakes in northeast China using Sentinel-2 images with machine learning algorithms[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 9132-9146. doi: 10.1109/JSTARS.2021.3109292
    [10] VANHELLEMONT Q, RUDDICK K. Acolite for Sentinel-2: Aquatic applications of MSI imagery[C]// Proceedings of the 2016 ESA Living Planet Symposium. Prague: 2016: 9-13.
    [11] VANHELLEMONT Q. Sensitivity analysis of the dark spectrum fitting atmospheric correction for metre- and decametre-scale satellite imagery using autonomous hyperspectral radiometry[J]. Optics Express, 2020, 28(20): 29948. doi: 10.1364/OE.397456
    [12] VANHELLEMONT Q. Adaptation of the dark spectrum fitting atmospheric correction for aquatic applications of the Landsat and Sentinel-2 archives[J]. Remote Sensing of Environment, 2019, 225: 175-192. doi: 10.1016/j.rse.2019.03.010
    [13] VANHELLEMONT Q, RUDDICK K. Atmospheric correction of metre-scale optical satellite data for inland and coastal water applications[J]. Remote Sensing of Environment, 2018, 216: 586-597. doi: 10.1016/j.rse.2018.07.015
    [14] TANG R G, SHEN F, PAN Y Q, et al. Multi-source high-resolution satellite products in Yangtze Estuary: cross-comparisons and impacts of signal-to-noise ratio and spatial resolution[J]. Optics Express, 2019, 27(5): 6426. doi: 10.1364/OE.27.006426
    [15] LI P, KE Y H, BAI J H, et al. Spatiotemporal dynamics of suspended particulate matter in the Yellow River Estuary, China during the past two decades based on time-series Landsat and Sentinel-2 data[J]. Marine Pollution Bulletin, 2019, 149: 110518. doi: 10.1016/j.marpolbul.2019.110518
    [16] LUO W, SHEN F, HE Q, et al. Changes in suspended sediments in the Yangtze River Estuary from 1984 to 2020: Responses to basin and estuarine engineering constructions[J]. Science of the Total Environment, 2022, 805: 150381. doi: 10.1016/j.scitotenv.2021.150381
    [17] VANHELLEMONT Q. Daily metre-scale mapping of water turbidity using CubeSat imagery[J]. Optics Express, 2019, 27(20): A1372-A1399. doi: 10.1364/OE.27.0A1372
    [18] 唐军武, 田国良, 汪小勇, 等. 水体光谱测量与分析Ⅰ: 水面以上测量法[J]. 遥感学报, 2004(1): 37-44. doi: 10.11834/jrs.20040106
    [19] 国家环境保护总局. 地表水和污水监测技术规范: HJ/T 91-2002[S]. 北京: 中国环境出版社, 2002.
    [20] 国家环境保护局. 水质 悬浮物的测定 重量法: GB 11901-89[S]. 1990.
    [21] LUO W, LI R H, SHEN F, et al. HY-1C/D CZI image atmospheric correction and quantifying suspended particulate matter[J]. Remote Sensing, 2023, 15(2): 386. doi: 10.3390/rs15020386
    [22] SHI K, ZHANG Y L, ZHU G W, et al. Long-term remote monitoring of total suspended matter concentration in Lake Taihu using 250 m MODIS-Aqua data[J]. Remote Sensing of Environment, 2015, 164: 43. doi: 10.1016/j.rse.2015.02.029
    [23] ALLAM M, YAWAR ALI KHAN M, MENG Q. Retrieval of turbidity on a spatio-temporal scale using Landsat 8 SR: A case study of the Ramganga River in the Ganges Basin, India[J]. Applied Sciences, 2020, 10(11): 3702. doi: 10.3390/app10113702
    [24] KRATZER S, KYRYLIUK D, EDMAN M, et al. Synergy of satellite, in situ and modelled data for addressing the scarcity of water quality information for eutrophication assessment and monitoring of Swedish Coastal Waters[J]. Remote Sensing, 2019, 11(17): 2051. doi: 10.3390/rs11172051
    [25] LEE Z P, CARDER K L, ARNONE R A. Deriving inherent optical properties from water color: a multiband quasi-analytical algorithm for optically deep waters[J]. Applied Optics, 2002, 41(27): 5755-5772. doi: 10.1364/AO.41.005755
    [26] BALASUBRAMANIAN S V, PAHLEVAN N, SMITH B, et al. Robust algorithm for estimating total suspended solids (TSS) in inland and nearshore coastal waters[J]. Remote Sensing of Environment, 2020, 246: 111768. doi: 10.1016/j.rse.2020.111768
    [27] LIU H Z, LI Q Q, BAI Y, et al. Improving satellite retrieval of oceanic particulate organic carbon concentrations using machine learning methods[J]. Remote Sensing of Environment, 2021, 256: 112316. doi: 10.1016/j.rse.2021.112316
    [28] GHATKAR J G, SINGH R K, SHANMUGAM P. Classification of algal bloom species from remote sensing data using an extreme gradient boosted decision tree model[J]. International Journal of Remote Sensing, 2019, 40(24): 9412-9438. doi: 10.1080/01431161.2019.1633696
    [29] CAO Z G, MA R H, DUAN H T, et al. A machine learning approach to estimate chlorophyll-a from Landsat-8 measurements in inland lakes[J]. Remote Sensing of Environment, 2020, 248: 111974. doi: 10.1016/j.rse.2020.111974
    [30] CHEN T Q, GUESTRIN C. XGBoost: a scalable tree boosting system[C]// Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. New York: Association for Computing Machinery, 2016: 785-794.
    [31] NECHAD B, RUDDICK K G, NEUKERMANS G. Calibration and validation of a generic multisensor algorithm for mapping of turbidity in coastal waters[C]// Remote Sensing of the Ocean, Sea Ice, and Large Water Regions 2009. SPIE, 2009: 7473: 74730H.
    [32] CASTAGNA A, AMADEI MARTÍNEZ L, BOGORAD M, et al. Optical and biogeochemical properties of diverse Belgian inland and coastal waters[J]. Earth System Science Data, 2022, 14(6): 2697-2719. doi: 10.5194/essd-14-2697-2022
    [33] WONG T T. Performance evaluation of classification algorithms by k-fold and leave-one-out cross validation[J]. Pattern Recognition, 2015, 48(9): 2839-2846. doi: 10.1016/j.patcog.2015.03.009
    [34] ZHENG Z B, LI Y M, GUO Y L, et al. Landsat-based long-term monitoring of total suspended matter concentration pattern change in the wet season for Dongting Lake, China[J]. Remote Sensing, 2015, 7(10): 13975-13999. doi: 10.3390/rs71013975
  • 加载中
    Created with Highcharts 5.0.7访问量Chart context menu近一年内文章摘要浏览量、全文浏览量、PDF下载量统计信息摘要浏览量全文浏览量PDF下载量2024-052024-062024-072024-082024-092024-102024-112024-122025-012025-022025-032025-040Highcharts.com
    Created with Highcharts 5.0.7Chart context menu访问类别分布DOWNLOAD: 2.4 %DOWNLOAD: 2.4 %HTML全文: 85.6 %HTML全文: 85.6 %摘要: 12.1 %摘要: 12.1 %DOWNLOADHTML全文摘要Highcharts.com
    Created with Highcharts 5.0.7Chart context menu访问地区分布其他: 93.7 %其他: 93.7 %XX: 4.9 %XX: 4.9 %东莞: 0.1 %东莞: 0.1 %北京: 0.5 %北京: 0.5 %吉林: 0.1 %吉林: 0.1 %开封: 0.1 %开封: 0.1 %新乡: 0.1 %新乡: 0.1 %无锡: 0.1 %无锡: 0.1 %杭州: 0.1 %杭州: 0.1 %沈阳: 0.1 %沈阳: 0.1 %济南: 0.1 %济南: 0.1 %福州: 0.1 %福州: 0.1 %贵阳: 0.1 %贵阳: 0.1 %郑州: 0.3 %郑州: 0.3 %阳泉: 0.1 %阳泉: 0.1 %其他XX东莞北京吉林开封新乡无锡杭州沈阳济南福州贵阳郑州阳泉Highcharts.com
图( 12) 表( 4)
计量
  • 文章访问数:  1976
  • HTML全文浏览数:  1976
  • PDF下载数:  110
  • 施引文献:  0
出版历程
  • 收稿日期:  2023-09-28
  • 录用日期:  2024-01-11
  • 刊出日期:  2024-02-26
李媛媛, 沈芳, 陈嵩钰, 魏小岛. 基于XGBoost的内陆河湖浊度反演与长时序分析[J]. 环境工程学报, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
引用本文: 李媛媛, 沈芳, 陈嵩钰, 魏小岛. 基于XGBoost的内陆河湖浊度反演与长时序分析[J]. 环境工程学报, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
LI Yuanyuan, SHEN Fang, CHEN Songyu, WEI Xiaodao. Inversion and long-term series analysis of turbidity in inland rivers and lakes based on XGBoost[J]. Chinese Journal of Environmental Engineering, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125
Citation: LI Yuanyuan, SHEN Fang, CHEN Songyu, WEI Xiaodao. Inversion and long-term series analysis of turbidity in inland rivers and lakes based on XGBoost[J]. Chinese Journal of Environmental Engineering, 2024, 18(2): 398-408. doi: 10.12030/j.cjee.202309125

基于XGBoost的内陆河湖浊度反演与长时序分析

    通讯作者: 沈芳(1964—),女,博士,教授,fshen@sklec.ecnu.edu.cn
    作者简介: 李媛媛 (1998—) ,女,硕士研究生,51213904006@stu.ecnu.edu.cn
  • 1. 华东师范大学河口海岸学国家重点实验室,上海 200241
  • 2. 上海勘测设计研究院有限公司,上海 200050
  • 3. 中国长江三峡集团有限公司长江生态环境工程研究中心 (上海) ,上海 200050
基金项目:
上海市科委重点项目 (20dz1204700) ;中国长江三峡集团有限公司科研项目资助 (202103552)

摘要: 浊度是影响水下光场及营养盐循环的关键要素之一,浊度监测可以为河湖水质的污染防控和预警提供科学依据。以长三角示范区的典型河湖为研究区,使用实测数据构建浊度反演模型,并利用1984—2022年Landsat卫星数据分析了研究区河湖浊度的长时序动态变化。对比传统经验模型、半经验模型和机器学习模型,XGBoost机器学习模型精度最高 (R2为0.68,RMSE为4.78 NTU) 。浊度反演结果表明,近40年河流航道和淀山湖北部非渔场区域浊度上升了10%和12%,而元荡湖和大莲湖浊度下降了19%和27%,并且浊度随着建设用地面积的增加而增大;研究区浊度季节性变化显著,秋冬季平均浊度比春夏季高6 NTU,月平均浊度与月平均降水量负相关,相关系数为-0.61 (p<0.05) ,但与月平均风速没有显著的相关性。基于XGBoost的Landsat长时序浊度反演能够把握研究区浊度的时空变化趋势,明确水污染管理与治理方向,最终助力长三角一体化发展。

English Abstract

  • 浊度 (turbidity,T,NTU) 是水环境监测的重要参数,反映了光在水体传输中的受阻程度,浊度的变化能够表征水中悬浮颗粒物浓度的变化,对水中污染物质迁移具有重要指示作用[1-2],浊度还会通过影响水下光的传播改变水生生态系统的初级生产力和生态平衡[3],因此浊度对于水资源管理、水生态保护等具有重要意义。长三角地区水环境状况复杂多变,前人对长三角地区的水体浊度[4]、总悬浮物浓度 (total suspended matter concentration,TSM,mg∙L−1) [5]、透明度[6]等开展了众多研究,但研究多为大型河湖水体,中小型湖泊和河流研究较少,长三角示范区是引领长三角地区高质量一体化发展的重要区域,浊度反演研究可为示范区绿色发展提供重要依据。

    如何利用卫星数据准确估算水体浊度是国内外学者研究的重点,常见的有经验模型、半经验模型和机器学习模型。HOU等[5]运用MODIS 555 nm和645 nm波段表面反射率比值构建了总悬浮物浓度经验模型用于长江中下游流域大型湖泊和水库的浊度动态监测。DOXARAN等[7]基于MODIS 858 nm和645 nm波段表面反射率比值构建浊度反演经验模型,用于研究法国吉伦特河口最大浑浊带浊度动态变化。DOGLIOTTI等[8]从固有光学特性出发基于MODIS 645 nm波段遥感反射率 (remote sensing reflectance,Rrs,sr−1) 发展了适用于沿海和河口水域中低浊度的半经验反演模型。MA等[9]运用梯度增强决策树 (gradient boosting decision trees,GBDT) 估算了中国东北地区湖泊浊度。前人研究的浊度反演模型大多应用于开阔的近海河口水域,受人类活动影响较小,研究对象多为光学敏感性较高的高度浑浊水体。而本研究水体主要为中低浊度水体,且研究区内湖荡密布,航道众多,生产生活密集,水动力变化复杂,针对此类水体浊度的长时序监测研究尚不多见。

    本研究以长三角示范区内典型河湖为研究区,运用实测数据构建并对比了经验模型、半经验模型和机器学习模型三类浊度反演模型,随后基于最佳的XGBoost (eXtreme Gradient Boosting) 模型运用Landsat卫星数据反演研究区浊度,分析浊度的长期时空动态变化及其驱动因素。

    • 长三角示范区包括上海市青浦区、江苏省吴江区、浙江省嘉善县,研究区位于长三角示范区中部,包括淀山湖、元荡湖、汾湖、大莲湖、白蚬湖和太浦河等,研究区内河湖水系密布,相互连通,白蚬湖-淀山湖航道、汾湖-太浦河航道为研究区主要航道,研究区位置和采样点分布见图1。以9×9卫星影像像素选择4个目标区域,依次编号为:1-淀山湖入湖口、2-淀山湖出湖口、3-太浦河上游、4-太浦河下游。

    • Landsat TM/ETM+/OLI传感器获取了1982年以来时间连续的地表遥感影像,能够满足近40年研究区长时序浊度监测的应用需求。Landsat TM/ETM+波段范围为450~2 350 nm,分别有7个和8个波段,Landsat OLI波段范围为430~2 290 nm,共9个波段。本研究使用来自美国地质调查局 (United States Geological Survey, USGS,https://earthexplorer.usgs.gov/) 30 m空间分辨率的Landsat Level-1数据。剔除1984‒2022年研究区受云量、条带、耀光影响的数据,最终筛选得到323幅高质量影像。

      VANHELLEMONT和RUDDICK[10]提出的暗光谱拟合大气校正算法 (dark spectrum fitting,DSF) 能够对高空间分辨率影像进行准确的大气校正处理[11-13],经验证该方法得到的遥感反射率精度较高[14],能够应用于悬浮物浓度[15-16]、浊度[17]的反演研究。因此,本研究采用该方法对Landsat数据进行大气校正,处理工具为ACOLITE Python (Version 20210114.0) 。

    • 1) 现场实测数据。使用地物光谱仪 (ASD FieldSpec4 HR NG,北京理加联合科技有限公司) 获取水体高光谱数据,根据水面以上测量方法[18]计算水体高光谱遥感反射率,光谱分辨率为1 nm。使用多参数水质分析仪 (HYDROLAB HL7,美国HACH公司;YSI-EMS-KorEXO,美国Xylem公司) 现场测量水样获取水质数据,包括电导率、温度、pH、溶解氧、氨氮和浊度等。野外实验严格遵循《地表水和污水监测技术规范》[19]的要求。将采集的水样带回实验室,采用称重法[20]测量水体总悬浮物浓度。2022‒2023年进行了5次现场数据采集工作,共采集了126组实测高光谱遥感反射率和水质数据,剔除异常数据后得到实测浊度和高光谱遥感反射率匹配的有效数据81组 (见表1) 。

      2) 等效遥感反射率。为匹配Landsat大气校正后的多光谱遥感反射率和现场实测高光谱遥感反射率,根据Landsat不同传感器的光谱响应函数 (spectral response function,SRF) ,将实测高光谱遥感反射率转化为传感器等效波长处的遥感反射率[21],计算公式如式(1)所示。

      式中:Bi表示卫星传感器第i个波段;Rrs(Bi)是传感器第i个波段的等效遥感反射率,sr−1λ1λ2分别是第i个波段波长的最小值和最大值,nm;Rrs(λ)是波长λ处实测高光谱遥感反射率,sr−1S(λ)为传感器在λ处的光谱响应函数。

    • 研究区覆盖急水港桥、淀峰、太浦河桥和汾湖大桥4个国控水质监测断面,本研究搜集了以上4个站点2021年夏季到2022年秋季的水质监测数据,用于浊度反演模型的真实性检验,数据来源于国家地表水水质自动监测发布系统 (https://szzdjc.cnemc.cn:8070/GJZ/Business/Publish/Main.html) 。

      本研究收集了1984—2022年研究区的降水量 (mm) 和风速 (m∙s−1) 数据,用于分析降水量、风速与浊度季节性变化的关系。该气象数据来自欧洲中期天气预报中心 (European Centre for Medium-Range Weather Forecasts,ECMWF) 的ERA5-Land数据集(https://cds.climate.copernicus.eu/),ERA5-Land数据集收集了1950年至今的陆地月平均气象数据 ,空间分辨率为0.1°×0.1°,时间分辨率为1个月 (30 d) 。

    • 目前发展的浊度反演方法以经验模型为主,通过简单的数学函数,拟合浊度和遥感反射率之间的关系;半经验模型从辐射传输理论出发,计算过程较为复杂;机器学习可以很好地描述客观变量和特征变量之间的非线性关系,能够更加灵活地处理具有高维特征的数据集。下文分别构建经验模型、半经验模型和机器学习模型,对比分析三类模型精度,以获取最佳浊度反演模型。

    • 以多项式、指数、幂函数模型为例的经验回归模型由于构造简单,被广泛运用于浊度和总悬浮物浓度的遥感反演[5, 22]。常用单波段Bi[4, 23]、波段比(Bi/Bj)[23-24]或波段组合(Bi-Bj)/(Bi+Bj)[4, 9]构建浊度反演模型 (对于OLI传感器,BiBj表示OLI 2、3、4或5波段遥感反射率) 。利用相关系数选择与浊度相关性最高的光谱特征组合,分别使用多项式、指数、幂函数模型对浊度和光谱特征组合进行拟合。

    • 常见的半经验浊度反演模型有2个,一个是DOGLIOTTI等[8]从NECHAD等[2]的反射率模型发展的单波段浊度反演模型 (Semi-Empirical Single Band turbidity retrieval algorithm,SESB) ,另一个是基于LEE等[25]准分析算法 (Quasi-Analytical Algorithm,QAA) 的浊度反演模型。

      参考NECHAD等[2]的模型构建方法,对SESB模型引入系数BλT以提高模型的准确性和预测能力,模型如式(2)和式(3)所示。

      式中:T表示浊度,NTU;ρw(λ)表示波长λ处水体表面反射率,sr−1,与离水辐亮度 (water-leaving radiance,Lw(λ),W/(m2∙sr∙nm)) 和水面以上下行辐照度 (downwelling irradiance just above the water surface,E0+d,W/(m2∙nm)) 有关;AλTBλT是波长依赖的两个系数,由实测数据回归获得。Cλ是NECHAD等[2]通过“标准”固有光学特性校准的系数,当λ为655 nm时,Cλ为0.1686。

      BALASUBRAMANIAN等[26]已证明中低浊度水体的悬浮颗粒物浓度与665 nm处的颗粒物后向散射系数 (particulate backscattering coefficient,bbp,m−1) 存在幂函数关系。通过LEE等[25]的QAA算法获取OLI的bbp(655),对bbp(655)和浊度进行幂函数拟合,构建基于QAA的浊度反演算法。本研究所用的QAA版本为QAA_v6,bbp(655)计算步骤见表2

    • 机器学习模型具有良好的泛化能力和鲁棒性,受噪声的影响较小,以XGBoost为代表的机器学习模型被广泛应用于水质参数反演[27-28],对于小样本量数据适用性良好[29]。XGBoost是CHEN和GUESTRIN[30]提出的基于GBDT的一种新实现,通过加法策略将一组弱学习器组合成一个强学习器。XGBoost的弱学习器是决策树,决策树的目标函数值必须为最小值。模型训练时首先使用整个数据集拟合一个学习器,然后添加第二个学习器来拟合前一个学习器的残差,重复训练过程,直到满足停止标准,最终预测结果为每个学习器的预测结果之和,XGBoost模型如图2所示。

      对于由n个样本和m个特征组成的给定训练集,D=(xi,yi),(|D|=nxiRm, yiR),XGBoost模型可视为由t个回归树组成的加法模型,XGBoost模型的预测值如式(4)所示,目标函数如式(5)所示。

      式中:t是树的数量,yiˆyi分别是测量值和预测值,F是CART回归树的函数空间,fk是第k个独立树所代表的函数。式(5)中引入fk最小化目标,ˆyi(k1)表示第k-1次迭代时对第i个实例的预测,l(yi,ˆyi)度量预测值ˆyi和测量值yi之间的差异。Ω(fk)为正则化项,即每棵树的复杂度之和,它可以控制模型的复杂性,防止过拟合,γλ分别为模型的惩罚系数和L2正则项系数,M为叶节点数,ω为叶得分。

      模型的输入包括10个单波段或波段组合:OLI的B2B3B4B5波段,波段比B4/B2B4/B3B4/B5,以及波段组合(B4B2)/(B4+B2)、(B4B3)/(B4+B3)、(B4B5)/(B4+B5)。输出为浊度T。确定模型结构的主要参数包括学习器数量n_estimators和最大树深max_depth,通过10倍交叉验证的网格搜索方法对参数进行调整。为避免过拟合,将subsample和colsample_bytree两个参数设置为0.7。本研究中使用XGBoost的Python库来实现该算法 (XGBoost Version: 1.7.4;Python 3.9.12) 。

    • 为评估浊度反演模型精度,运用决定系数 (R2) 、均方根误差 (Root Mean Square Error,RMSE) 、平均绝对百分比误差 (Mean Absolute Percentage Error,MAPE) 作为评价指标。各指标公式如式(6)~式(8)所示。

      式中:Xmi为实测值,Xei为模型预测值,n为样本数量。

    • 研究区水体属于中低浊度水体,实测浊度为0~50 NTU,中位数为25 NTU。图3为实测浊度的光谱响应,浊度较高的水体高光谱遥感反射率也较高,在可见光波段,560~580 nm和700~720 nm处分别存在明显反射峰,反射谷则落在670~680 nm,在近红外波段,810~820 nm处存在一个较小反射峰,说明红、绿波段是浊度的敏感波段,可能是由悬浮颗粒物散射造成的[5]图4显示浊度和总悬浮物浓度线性相关,相关系数r为0.69 (p<0.01) ,这与已有的研究结果一致[7, 31-32]

    • 表3显示了Landsat OLI遥感反射率不同波段组合与浊度的相关性。单波段中,B4与浊度相关性最高,相关系数为0.71;波段组合中,B4/B3和与浊度相关性最高,相关系数为0.88,证明OLI红绿波段遥感反射率比值对浊度较为敏感。

      利用三类方法 (第2.1至2.3节所述) 构建浊度反演模型。由于实测数据量有限,在构建模型时使用留一交叉验证[33] (Leave One Out Cross Validation,LOOCV) 对模型进行评估,此方法适用于样本量较少的数据集,在完成交叉验证后,将所有数据用于模型构建。表4列出了三类浊度反演模型的交叉验证结果。基于B4/B3的XGBoost浊度反演模型精度最优,R2为0.83,RMSE为4.58 NTU;基于B4/B3波段组合的二次多项式、指数、幂函数回归模型的R2均为0.78,RMSE为5.25~5.26 NTU;基于QAA的半经验模型R2为0.54,RMSE为7.60 NTU,基于SESB的半经验模型R2为0.39,RMSE为8.72 NTU;经验模型和半经验模型效果均低于XGBoost机器学习模型精度,说明XGBoost模型对于处理小样本量数据具有一定优势。

    • 获取采样点同步卫星数据遥感反射率时,采样时间与卫星过境时间相差1.5 h以内,以采样点所在像元为中心使用3×3窗口进行卫星遥感反射率平均。图5显示OLI传感器实测等效遥感反射率与Acolite大气校正的遥感反射率相关系数达0.94 (p<0.01) 。

      将浊度反演模型应用到卫星数据时 (图6) ,浊度实测值和预测值大部分处于1∶1线附近,浊度大于30 NTU的点预测偏差较大。基于XGBoost的浊度反演模型效果最好,R2为0.68,RMSE为4.78 NTU;其次是二次多项式、指数、幂函数经验模型,R2为0.5~0.51、RMSE为5.93~6.03 NTU;半经验模型效果相对较差,基于SESB的半经验模型R2为0.13,RMSE为7.91 NTU,基于QAA的半经验模型R2为0.44,RMSE为6.37 NTU。验证结果表明XGBoost机器学习方法能够一定程度上提高模型精度。

    • 基于XGBoost模型,利用1984—2022年Landsat TM/ETM+/OLI卫星数据反演研究区浊度,分析浊度长期时空动态变化,结合自然因素和人类活动分析浊度变化的驱动因素。

      1) 浊度的年际时空变化。研究区浊度的长时序时空变化分析以3年为周期求平均值[16],结果见图7。1993—1995年、2014—2019年汾湖-太浦河航道浊度较高,可能与1991年底开启的太浦河疏浚工程和2014年底动工的金泽水库开发工程有关。1996年汾湖一分为二,内汾湖浊度低于外汾湖。2005年淀山湖北部渔场开发,该区域在2005‒2007年浊度显著上升。图8显示了研究区近40年浊度整体变化趋势,研究区河流航道、淀山湖北部非渔场区域浊度升高,分别上升了10%和12%,可能与水上货运量增多和区域建设有关。元荡湖、大莲湖浊度下降,分别下降了19%和27%,可能与河湖湿地生态保护修复有关。

      定量分析4个目标区域 (图1) 近40年浊度变化,结果如图9所示,目标区域近40年浊度在15~40 NTU之间,2017年以前,淀山湖入湖口和出湖口浊度差异较小,2017年以后,入湖口浊度大于出湖口,浊度平均相差10 NTU,太浦河上下游浊度差值保持在0~4 NTU之间,浊度差异较小。土地利用分类结果显示,太浦河沿岸建设用地面积近40年增长了1.5倍,太浦河上游浊度和建设用地面积相关系数为0.71 (p<0.01) (图10) ,说明建设用地的增加在一定程度上会影响研究区浊度变化。

      2) 浊度的季节性变化。研究浊度的季节性变化时将近40年浊度反演结果按季节求平均值,如图11所示。研究区近40年浊度季节性变化显著,秋冬季平均浊度为27 NTU,春夏季平均浊度为21 NTU,秋冬季平均浊度比春夏季高6 NTU。白蚬湖-淀山湖航道和汾湖-太浦河航道冬季浊度明显高于其他季节。

      降水和风力是调节水体浊度的重要因素[22, 34],1984—2022年月平均浊度和降水量、风速之间的关系如图12所示。近40年研究区1~3月、9~12月月平均浊度较高 (>25 NTU) ,而4~8月浊度较低,月平均浊度秋冬季高于春夏季,这与HOU等[5]的研究结果一致。近40年月平均降水量在0~10 mm之间,6月降水量最高为9 mm,浊度与降水量整体呈负相关关系,相关系数为−0.61 (p<0.05) 。1984‒2022年月平均风速在1~2 m·s−1之间,7月风速最大为1.8 m·s−1,浊度与风速之间无显著相关性。

    • 1) 与经验、半经验模型相比,XGBoost在模型构建和真实性检验时具有较高的精度 (模型构建时R2为0.83,RMSE为4.58 NTU;真实性检验时R2为0.68,RMSE为4.78 NTU) ,适用于浊度反演研究。

      2) 近40年研究区浊度时空变化明显。人类活动较为频繁的河流航道和淀山湖北部非渔场区域浊度上升,而位于湿地生态保护修复区的元荡湖、大莲湖浊度下降,建设用地面积增加一定程度上使得研究区浊度上升。

      3) 近40年研究区浊度季节性变化明显。秋冬季浊度高于春夏季,月平均浊度与月平均降水量相关性较为明显,相关系数为-0.61 (p<0.05) ,而与月平均风速无明显相关性。

    参考文献 (34)

返回顶部

目录

/

返回文章
返回