EN中文

地球信息科学学报 >

2023 , Vol. 25 >Issue 8: 1682 - 1698

DOI: https://doi.org/10.12082/dqxxkx.2023.221010

遥感科学与应用技术

无人机遥感和GWR结合的水华短时预测方法

  • 张寒博 , 1 ,
  • 李彤 2 ,
  • 李晓芳 2 ,
  • 邓滢 2 ,
  • 邓应彬 1 ,
  • 荆文龙 , 1, * ,
  • 胡义强 1 ,
  • 李勇 1 ,
  • 杨骥 1
展开
  • 1.广东省科学院广州地理研究所，广州 510070
  • 2.广东省生态环境监测中心，广州 510308
*荆文龙（1988— ），男，河南开封人，博士，副研究员，主要研究方向为水文水资源遥感。E-mail:

张寒博（1996— ），男，河南开封人，硕士，主要研究方向为水环境遥感。E-mail:

收稿日期: 2022-12-27

  修回日期: 2023-03-20

  网络出版日期: 2023-07-14

基金资助

国家重点研发计划项目(2022YFF0711602)

国家自然科学基金项目(41976190)

国家自然科学基金项目(41976189)

广东省科学院发展专项资金项目(2022GDASZH-2022010202)

广东省科学院发展专项资金项目(2022GDASZH-2022020402-01)

广东省科技计划项目(2021B1212100006)

收起

A Short-term Prediction Method of Water Bloom based on UAV Remote Sensing and GWR

  • ZHANG Hanbo , 1 ,
  • LI Tong 2 ,
  • LI Xiaofang 2 ,
  • DENG Ying 2 ,
  • DENG Yingbin 1 ,
  • JING Wenlong , 1, * ,
  • HU Yiqiang 1 ,
  • LI Yong 1 ,
  • YANG Ji 1
Expand
  • 1. Guangzhou Institute of Geography, Guangdong Academy of Sciences, Guangzhou 510070, China
  • 2. Ecological Environment Monitoring Center of Guangdong Province，Guangzhou 510308, China
*JING Wenlong, E-mail:

Received date: 2022-12-27

  Revised date: 2023-03-20

  Online published: 2023-07-14

Supported by

National Key R&D Program of China(2022YFF0711602)

National Natural Science Foundation of China(41976190)

National Natural Science Foundation of China(41976189)

The GDAS' Project of Science and Technology Development(2022GDASZH-2022010202)

The GDAS' Project of Science and Technology Development(2022GDASZH-2022020402-01)

The Science and Technology Program of Guangdong(2021B1212100006)

Fold

摘要

目前水华已成为国内外一个重要的环境问题，单一的光学遥感方法难以实现精细化的水华预测。针对上述问题，本文以无人机多光谱影像、水质、水温及气象数据为数据源，首先通过像素匹配法（Matching Pixel-by-Pixel，MPP）反演水质参数及归一化植被指数（Normalized Difference Vegetation Index, NDVI）阈值法提取水华信息，然后基于地理加权回归（Geographically Weighted Regression, GWR）建立涵盖面积信息和位置信息的水华短时预测模型，并探讨了预测窗口尺寸对预测结果的影响。结果表明：① 本文方法可以实现精细化的水华短时预测，在应用中对水华面积预测精度达到96.19%，在水华空间分布的预测上对水华和非水华的总体分类精度均大于0.97，生产者精度、Kappa系数均大于0.5； ② MPP反演的总氮、总磷和溶解氧浓度与实测数据都有着较高的相关性，其决定系数R2分别达到0.89、0.85、0.89； ③ 不同尺寸的预测窗口直接影响预测结果精度，相比8×8、12×12、14×14，选用10×10预测窗口得到的总体分类精度、生产者精度、Kappa系数最高，分别为0.98、0.77、0.77。该模型可为短时水华预测提供借鉴，有利于水华监测、预测工作的进一步深入与完善。

关键词： 水华; 预测模型; 无人机遥感; MPP; GWR; 水质参数

本文引用格式

张寒博 , 李彤 , 李晓芳 , 邓滢 , 邓应彬 , 荆文龙 , 胡义强 , 李勇 , 杨骥 . 无人机遥感和GWR结合的水华短时预测方法[J]. 地球信息科学学报, 2023 , 25(8) : 1682 -1698 . DOI: 10.12082/dqxxkx.2023.221010

Abstract

In recent years, the conflict between rapid growth of economy and protection of water resources has become increasingly prominent. Water bloom has also become an important environmental issue, both domestically and internationally. Efficient use of remote sensing to predict water bloom outbreaks is of great significance to promote the management and protection of the lake and reservoir water environment. Using multispectral images of Unmanned Aerial Vehicle and measured water quality parameters as data sources, the water quality parameters were inversed by the Matching Pixel-by-pixel (MPP) algorithm, and the water bloom information was extracted through the Normalized Difference Vegetation Index (NDVI) threshold method. Then, a short-term prediction model of water bloom based on geographical weighted regression was proposed. The characteristic of this model was that it can accurately estimate the area and spatial location of the water bloom within a short period. We also discussed the influence of window size on the prediction results. The results show that: (1) based on the proposed model, the prediction accuracy of bloom area reached 96.19%. For the spatial distribution of water bloom, the overall classification accuracy of water bloom and non-water bloom pixels was both greater than 0.969, and the producer accuracy and Kappa coefficient were both higher than 0.5; (2) the total nitrogen, total phosphorus, and dissolved oxygen concentrations inverted using the MPP algorithm showed high correlation with the measurements (R2 was 0.888 6, 0.854 6, and 0.896 9, respectively). The water bloom information extracted by NDVI threshold method had a high consistency with the true color orthophoto-based visual results in terms of bloom outbreak density and spatial distribution; (3) the prediction window size was closely related to the spatial resolution of data. The prediction window of different sizes directly affected the accuracy of water bloom prediction results. Compared with the window size of 8×8, 12×12, and 14×14, the highest accuracy (0.98), precision (0.77), and Kappa (0.77) were obtained with a window size of 10×10. Hence, the 10×10 prediction window was the most suitable for the water bloom prediction in this study area. The model developed here can effectively predict water bloom in lakes and reservoirs in advance, which provides a reference for improving short-term water bloom prediction and warning.

Key words： water bloom; prediction model; UAV remote sensing; MPP; GWR; water quality parameters

1 引言

近年来，经济社会的快速发展和水资源保护间的矛盾日益突出，水源地内农场养殖、居民生活排污、工业废水排放等行为时有发生，水源地富营养化程度日趋严重[1]。富营养化引起的水华现象不仅影响环境景观，还对湖库系统的水质健康和生态平衡产生了严重危害[2⇓-4]。提前对水华爆发强度及爆发位置进行预测，对防止水华爆发造成的灾害、制定合理的防治措施具有重要意义。
研究表明氮、磷等营养盐浓度，光照、水温等影响藻类的生长速度[5⇓-7]，强风、强降雨、洪水等极端气候会导致水华的垂直结构发生改变[8-9]，使本存在于水体中的蓝藻群体上浮，聚集到水体表面形成可见的水华。
目前预警预报水华爆发的主要方法是通过监测站对水质指标进行实时监测[10]，但这种监测结果只能通过站点插值的方式反映水华的变化和情势，不能掌握水华的爆发程度和空间分布。相较于传统的地面监测，遥感技术为水华的监测预测提供了新的机遇，可以实现对水华的长期、快速和动态的监测，并结合地面实测数据实现水华的短期预测[11⇓-13]。
卫星遥感具有覆盖范围广、成本低等特点，很多产品也相继应用于水华的监测中，如MODIS系列[14-15]、Landsat系列[16]、高分系列[17]等。但水华爆发迅速，卫星遥感的空间分辨率和时效性难以满足水华预测的需求，并且南方地区多云多雨的气候条件也限制了卫星遥感的应用。相比卫星遥感，无人机遥感机动灵活、时效性强、分辨率高，能够有效开展水华的实时、快速监测，如Zhao等[18]基于无人机多光谱数据，通过构建归一化植被指数（Normalized Difference Vegetation Index, NDVI），采用阈值分割法实现了对滇池水华的监测；丁铭等[19]通过无人机高光谱数据的光谱差异实现了水华的提取。
目前对水华进行预测的研究普遍是基于卫星遥感影像的水华面积预测或通过水华爆发概率间接预测水华空间分布，如张虎军等[14]结合实测数据和卫星遥感数据构建多元线性回归预测模型，实现了太湖年度水华的面积预测；吴娟等[15]通过太湖气象水文水质数据与MODIS数据，采用支持向量机、长短记忆神经网络、极端梯度提升树模型成功对太湖水华面积进行了预测；祁国华等[16]通过构建气象数据和遥感数据间的非线性关系实现对巢湖水华发生概率的预测；刘东等[17]构建了一种基于遥感藻总量和气象因子的不同湖区水华暴发概率预测方法，模拟了巢湖的水华分布。基于无人机遥感的研究则还处在水华的提取和监测阶段，还未发现精细化的水华预测方法。地理加权回归模型（Geographically Weighted Regression，GWR）是一种空间分析技术，广泛应用于地理学及空间模式分析的相关学科，它考虑到了空间对象的局部效应，通过建立水质、水文、气象数据与水华关系，为实现水华预测提供了可能[20]。
本文从水华生长环境因素和移动规律入手，以无人机多光谱影像和实测水质参数为数据源，先通过像素匹配法（Matching Pixel-by-pixel，MPP）反演水质参数及NDVI阈值法提取水华信息，后基于GWR实现了水华面积和位置的预测，提出了一种涵盖面积信息和位置信息的水华短时预测模型，在广东省鹤地水库开展了实验和应用，可为其他湖库水华的遥感预测提供参考。

2 研究方法

针对目前缺乏精细化的水华预测方法，本文建立了一种水华短时预测模型，其设计思路为：① 水华面积预测：通过MPP算法进行水质参数反演，构建预测窗口，结合GWR实现水华面积预测；② 水华位置预测：通过气象数据与水华面积的建立位移函数，实现预测结果可视化。模型构建流程如图1所示。
图1 研究技术路线

Fig. 1 Work flow of the study

2.1 水质参数模型构建及水华信息提取

无人机影像的空间分辨率一般为厘米级，单个像元对应的水面区域难以准确匹配采样点，因此，本文选择MPP算法进行多光谱反射率与水质参数的拟合，采用3×3的像素窗口进行模型的构建。
MPP是一种通过遍历采样点对应像素窗口中所有的像素值来寻找最优回归模型的反演算法[21]。其原理如图2所示，对m个采样点对应的i×j像素窗口中共（i×j）m组数据进行4种函数拟合，从而找出最优模型（以R2为依据）。
图2 MPP算法原理

Fig. 2 MPP algorithm principle

为避免陆地及植被等因素对水质参数遥感定量反演结果的影响，通过无人机多光谱影像构建 NDVI(> 0.3，< 0.8)和归一化水体指数（Normalized Difference Water Index，NDWI > 0.2）对植被及水体信息进行精确提取，在此基础上结合阈值法和人工目视判读实现对水华信息的提取，其计算公式如下：
N D W I = b 2 - b 5 b 2 + b 5
N D V I = b 5 - b 3 b 5 + b 3
式中： b 2、 b 3、 b 5分别代表无人机多光谱数据的绿、红、近红外波段的反射率，本文水华信息提取的阈值通过参考潘鑫[22]的阈值结合目视判读最终确定为-0.081~0.264。

2.2 水华面积与位置预测

水华面积的预测是通过引入预测窗口建立不同时相（相同时间步长）的水华面积与水质参数间的函数关系来实现。其中，时间步长的确定需要在模拟效率和预测精度间进行权衡，预测窗口本质上是大小为n×n的局部像元矩阵，预测窗口大小需根据遥感影像的空间分辨率进行调整。如 图3所示，预测窗口是以遍历方式进行移动，水华的面积为预测窗口内水华像元数目与其空间分辨率的乘积，各水质参数为预测窗口内所有值的算术平均数，公式如下：
S = n ' × C 2
α = ∑ i = 1 n 2   α i n 2
式中： n '代表预测窗口内水华的像元个数；C代表遥感影像的空间分辨率/m。本文中C为0.09， α可代表TN、TP、pH、DO、水温（Q）。
图3 预测窗口与面积预测原理

Fig. 3 Prediction window and area prediction principle

GWR（图4）是一种用于量化空间异质性的局部线性参数估计方法[23]，通过引入数据的位置信息来计算各参数在局部回归方程里的权重。其基本公式如下：
y i = β 0 u i , v i + ∑ t = 1 n β t u i , v i x i t + ε u i , v i                                       ( i = 1,2 , ⋯ , m )
式中： y i为第 i个基本单元的水华面积，共有 m个点； x i t为第 t个参数的第 i个基本单元的值； u i , v i表示第 i个基本单元的经纬度； β 0 u i , v i为第 i个基本单元的常数值； β t ( u i , v i )为第 t个参数对第 i个基本单元的线性回归系数， ε u i , v i为模型在第 i个基本单元的残差值。
图4 GWR原理

Fig. 4 GWR Principle

本文使用的GWR选择交叉验证（Cross-validation，CV）方法来确定带宽，选用高斯函数法计算空间权重矩阵，构建出一个单调递减的连续性函数式，公式如下：
C V = 1 n ∑ i = 1 n y i - y ^ ≠ i b 2
W i t = e x p - D i t b 2
W i t → β t u i , v i x i t
式中： y ^ ≠ i b表示在回归参数估计时只根据回归点周围的数据点进行回归计算，当 C V值到达最小时，对应的 b为最佳带宽； W i t为权重，是以 e为底的指数函数，带宽越大，权重随距离增加衰减的越慢，带宽越小，权重随距离增加衰减的越快； D i t为采样点与带宽范围内其余点之间的距离。权重与回归系数的具体转换公式见文献[24]。
基于GWR，以预测窗口为基本单元，时间步长为3.5 h，结合 T 1时刻水华的面积、总氮、总磷、pH（无量纲）、溶解氧、水温、经纬度8个参数来实现 T 2时刻水华面积的预测，具体步骤如下：
（1）将 T 2时刻提取的水华面积为因变量，各参数为自变量，建立GWR训练模型。从回归模型中得到常数项，自变量系数及残差结果。
（2）将训练模型结果栅格化，并自变量系数、常数项重采样为0.09 m，对残差结果通过反距离权重法插值得栅格残差。
（3）按照GWR原理将训练模型结果进行回带建立预测模型，将待预测数据的各自变量参数带入模型并与残差相加，得到最终的水华面积预测结果。
S ' T 2 = f S ' T 1 , T P ' , T N ' ,   p H ' , D O ' ,   Q ' , L o n ' , L a t ' + ε
S T 2 = f S T 1 ,   T P ,   T N ,   p H ,   D O ,   Q ,   L o n ,   L a t + ε
式中：有上标“ '”表示训练模型数据集，没有上标“ '”表示预测模型数据集； S T 1 、 S T 2分别为T1、T2时刻水华的面积； T P为水域总磷含量/（mg/L）； T N为总氮含量/（mg/L）； p H为酸碱度； D O为溶解氧含量/（mg/L）； Q为水温/ ℃; L o n 、 L a t为经纬度； ε为残差; f (   )为非线性函数。
考虑到晴朗天气下水华位置的变化主要受风和水流的影响，本文水华位置的预测仅考虑这两种因素，基于相同时间间隔的多组风速、风向、流速、水华面积与水华空间坐标间的关系建立位移函数来实现水华位置的预测。论文水流方向视为由上游到下游，并以水流方向为x轴、其垂直方向为y轴建立相对坐标系，其公式如下：
x 2 = x 1 + ∫ t 1 t 2 V 风 s i n θ α 1 d t + ∫ t 1 t 2 V 水 b 1 d t
y 2 = y 1 + ∫ t 1 t 2 V 风 c o s θ α 2 d t         ( ∫ t 1 t 2 V 风 c o s θ α 2 d t < l ) y 1 + l                                                           ( ∫ t 1 t 2 V 风 c o s θ α 2 d t ≥ l )
式中： x 1 , y 1为 T 1时刻水华的坐标； x 2 , y 2为 T 2时刻水华的坐标； V 水为水流速度； V 风为风速； θ为风向； l为水华位置到河岸的距离； α 1 、 α 2 、 b 1为速度系数，是基于基本单元通过不同时段的多组数据训练得到。本研究最终训练得到 α 1 = 11.47 S T 1、 α 2 = 3.57 S T 1、 b 1 = 67.19 S T 1。

2.3 预测结果可视化

预测结果的可视化（图5）是通过 T 2时刻提取的水华位置，建立半径为r的缓冲区实现，其公式如下：
r = S T 2 π
式中： r为预测半径，是以预测点（ T 2时刻基本单元的中心）建立的缓冲区半径； S T 2为 T 2时刻水华的面积。
图5 预测结果可视化

Fig. 5 Visualization of prediction results

2.4 精度评价

引入混淆矩阵把水华预测结果分为水华和非水华2个类别，通过生产者精度、总体分类精度、Kappa系数3个指标来对预测结果进行评价。其中，生产者精度表示被预测为是水华的样本中，真正的水华所占的比例；总体分类精度表示预测结果与实际提取结果具有一致性的概率；Kappa系数是一种衡量一致性的多变量分析方法，常用于评价遥感影像分类的正确程度[25]。公式如下：
P p = m 11 m 11 + m 12
P 0 = m 11 + m 22 n
P e = m 11 + m 12 m 11 + m 21 + m 21 + m 22 m 12 + m 22 n 2
K a p p a = P 0 - P e 1 - P e
式中： P p为生产者精度； P 0为总体分类精度； P e为错分误差； m 11 、 m 22表示2个类别正确预测的样本数目； m 12表示未把水华预测为非水华的样本数目； m 21表示未把非水华预测为水华的样本数目；n为样本总数。Kappa系数的值域为-1~1，值越大说明预测结果越准确。

3 实验区概况及数据来源

3.1 实验区概况

鹤地水库位于湛江市河唇镇，跨越玉林、茂名、湛江3市，是广东省主要饮用水源地之一，总库容11.51×108 m3，水面面积108.4 km2，正常蓄水位40.5 m[26]。水库年均供水量约1.55×1010 m3，是一座以灌溉和供水为主，兼顾发电、防洪、航运的大型综合性水库。研究选取其下游石角站附近作为研究区，空间分布如图6所示，石角站水域历史水质常处于Ⅱ—Ⅲ类，是鹤地水库水华最易爆发的区域之一。
图6 研究区空间分布

Fig. 6 Spatial distribution of the study area

3.2 数据来源

3.2.1 遥感数据

无人机多光谱影像通过小型四旋翼航测无人机（P4 Multispectral和 MAVIC 2 Enterprise）（表1）采集获取。P4 Multispectral搭载一台一体式多光谱相机，包含1个用于可见光成像的彩色传感器和5个用于多光谱成像的单色传感器（Blue、Green、Red、Red Edge、NIR），每个传感器有效像素为208万，MAVIC 2 Enterprise搭载红外传感器和可见光传感器。为了避免水华垂直结构的变化对预测结果产生影响，本文选择晴空、少云、风速较小的天气进行数据采集，2台无人机各飞行6个架次，飞行高度设为150 m，飞行速度10 m/s，拍照间隔2 s，空间分辨率0.09 m，航向重叠率80%，旁向重叠率60%，采集区域面积为4 499.81 m²。
表1 无人机传感器参数及影像采集时间

Tab. 1 UAV sensor parameters and image acquisition time

传感器 波段 波段名称 波长范围 采集时间（2022年夏季某月）
P4 Multispectral多光谱相机 b1 Blue 450 nm ± 16 nm 21日10:30
b2 Green 560 nm ± 16 nm 21日14:00
b3 Red 650 nm ± 16 nm 21日17:30
b4 Red Edge 730 nm ± 16 nm 22日10:30
b5 NIR 840 nm ± 26 nm 22日14:00
MAVIC 2 Enterprise红外相机 b1 IR 8~14 μm 22日17:30
利用Pix4d软件对获取的无人机遥感数据进行几何校正、地理定位、图像拼接等预处理，通过标准标定板进行多通道光谱影像辐射校正和辐射定标，生成具有地表真实反射率的正射影像。研究区多光谱正射影像如图7所示。
图7 研究区多光谱正射影像

Fig. 7 Multispectral orthophotography of the study area

3.2.2 实测数据

在每次采集无人机多光谱影像时同步进行水质采样（表2），在研究区内均匀选取10个采样点（包括6个实验点和4个验证点），于水面以下50 cm位置进行水质采样，同时用GPS进行定位，共采样60 次，采样点的空间分布如图8所示。水质采样仪器使用AMT-YB101型便携式多参数水质检测仪，采集数据主要包含pH、溶解氧（Dissolved Oxygen, DO）、总氮（Total Nitrogen, TN）、总磷（Total Phosphorus, TP）等。为探究水华的移动规律，对研究区内12处水华区域进行标记（图8），并记录不同时刻被标记水华的位置信息。风速、风向、流速等数据来自鹤地水库各监测站和中国气象数据网（http://data.cma.cn/）。
表2 实测数据（部分）

Tab. 2 Measured data (part)

时间 pH DO/（mg/L） TN/（mg/L） TP/（mg/L） 风速/(m/s) 风向/° 流速/（m/s）
21日10：30 8.12 7.54 0.87 0.041 1.2 264.6 0.21
21日10：30 8.16 7.68 0.84 0.050 1.2 267.5 0.23
21日10：30 8.23 8.56 0.91 0.047 1.2 266.1 0.21
21日10：30 8.17 5.96 0.85 0.042 1.2 263.8 0.21
21日10：30 8.29 9.28 0.87 0.048 1.2 267.6 0.21
21日14：00 8.18 7.96 0.82 0.051 1.5 280.3 0.12
21日14：00 8.19 8.13 0.88 0.046 1.5 278.7 0.11
21日14：00 8.26 9.21 0.86 0.042 1.5 278.6 0.12
21日14：00 8.20 10.02 0.81 0.051 1.5 279.5 0.12
21日14：00 8.27 8.64 0.92 0.047 1.5 280.1 0.11
21日17：30 8.19 7.95 0.84 0.043 1.3 245.6 0.16
21日17：30 8.18 8.97 0.82 0.049 1.3 247.8 0.16
21日17：30 8.21 7.69 0.84 0.049 1.3 246.6 0.17
21日17：30 8.25 9.98 0.90 0.047 1.3 243.9 0.16
21日17：30 8.26 9.63 0.87 0.045 1.3 246.3 0.16
图8 水质采样点及水华标记位置

Fig. 8 Water quality sampling point and bloom marking position

4 结果及分析

4.1 水质参数反演结果

水质参数与水华的生长代谢有着紧密的关系，是水华面积预测的基础。本文基于获取的无人机多光谱影像，通过每个采样点的光谱反射率数据得到水体光谱曲线（图9）。可以看出各采样点光谱曲线有一定差别，但整体都呈现出一般内陆水体的光谱特征，在460~530 nm波段区间变化基本趋于一致，表现为反射率与波长呈正相关变化；在560 nm附近出现反射峰，这是因为水体中叶绿素和其他可溶性有机物在该波段内的吸收率较弱；在560~710 nm 反射率与波长呈负相关变化，710 nm附近出现反射率谷，主要归因于水体中藻蓝蛋白色素及叶绿素的强吸收。
图9 研究区水体光谱曲线

Fig. 9 Spectral profiles of water bodies in the study area

在构建水质参数反演模型时，需要选取与实测水质参数拟合度最高的波段或波段组合作为模型输入，本文基于6个采样点实测水质参数，共建立了47种不同的反射率波段或波段组合来进行相关性分析，最终筛选出与水质参数值最显著相关的 (b2+b3)/b5（r=0.453）、b1×b4（r=0.496）、(b3-b2)/b4（r=0.517）分别作为反演TP、TN、DO的特征波段反射率。通过选出的最佳波段作为自变量，分别建立线性函数模型、二次函数模型、指数函数模型和幂函数模型，最后根据R2选择出最优模型作为水质参数反演模型。3种水质参数的最优模型表达式如表3所示。
表3 水质参数模型表达式及精度

Tab. 3 Water quality parameter model expression and precision

水质参数 最佳波段 r 模型表达式 R2
TP (b2+b3)/b5 0.45 y=0.07839e1.0483x 0.78
TN b1×b4 0.50 y=0.04872x2-0.02873x+0.08913 0.80
DO (b3-b2)/b4 0.52 y=0.07254x2-0.00871x+0.96418 0.72
为了后期更好地提取和预测水华，需要降低背景信息的干扰，排除非水体区域及水体区域内水生植被、船舶等因素的影响。利用无人机多光谱数据构建NDWI，通过NDWI > 0.2的阈值实现水体的精准提取，并在此基础上完成水质参数的反演。其中，TP、TN、DO通过水质参数反演模型获取，pH通过样点空间插值获取，NDVI通过波段计算获取，水温数据借助DJI Thermal Analysis Tool 3工具处理红外影像得到。以21日14：00的数据为例，各水质参数的空间分布如图10所示，可以看出研究区水质整体呈弱碱性，pH为8.16~8.42，表现出由上游到下游递减；DO为1.402~14.13 mg/L，在近岸及水生植被周边浓度整体偏低；TP为0.0065~0.0804 mg/L，浓度呈现出由两岸向中心递减趋势，TN为0.032 ~ 1.094 mg/L，高浓度主要分布在上游西部和下游东部；NDVI为-0.747~0.728，空间分布特征与TN相似；水温为27.34~34.68℃，除下游西部略低外，整体分布较均匀。
图10 研究区水质参数空间分布（21日14:00）

Fig. 10 Comparison of cyanobacteria bloom area extraction results with true color orthophoto (10:30, 22nd)

通过6个时刻验证采样点的实测水质参数与相应位置处的反演结果进行线性拟合分析，以此验证TN、TP、DO的反演精度。由图11可以看出，3种水质参数反演结果与实测数据都有着较高的相关性，其R2分别为0.89、0.85、0.89，RMSE分别为0.018、0.003、0.649，说明水质反演结果具有较高的可信度，可以为水华预测提供数据支持。
图11 实测TN、TP、DO与反演结果散点图

Fig. 11 Scatter plot of measured TN, TP, DO and inversion results

4.2 水华面积提取结果

通过构建NDVI，采用阈值法结合人工目视判读对水华面积进行提取，以无人机多光谱数据合成的真彩色正射影像为底图，对提取结果进行验证。由图12（以22日10：30的数据为例）可以看出，研究区内水华提取结果的空间分布主要以近岸水面为主，在a、b、c、d区域内分布较为集中。通过将提取结果与正射影像进行叠加，可以看出提取结果无论是从水华爆发程度还是空间分布上都与正射影像目视结果具有较高的一致性，通过这种方法可以有效提取水华的空间分布与面积，并有较高的可信度。
图12 水华面积提取结果与真彩色正射影像对比（22日10:30）

Fig. 12 The results of bloom area extraction are compared with true color orthophoto (10:30, 22nd)

图13是通过阈值法结合人工目视判读提取的21日10:30—22日17:30共6个时刻的水华，表4为不同时段影像提取的水华信息。整体来看，水华的爆发具有一定规律性，表现为在10:30—14:00爆发程度要高于14:00—17:30，21日10:30—14:00水华面积增长11.25 m2，14:00—17:30水华面积增长7.67m2，22日10:30—14:00水华面积增长12.64 m2，14:00—17:30水华面积增长9.04 m2。水华的空间分布具有趋流性，因为其主要在水流和风力的作用下进行移动，21日10:30水华聚集区主要分布在研究区上游西侧，14:00水华逐渐向东南方向扩散， 18:00水华聚集区由上游转移到下游。22日10:30水华分布以上游近岸及下游西侧为主，14:00向东扩散，17:30主要集中分布在下游。
图13 水华提取结果（21日10:30—22日17:30）

Fig. 13 Bloom extraction results (10:30, 21st—17:30, 22nd)

表4 水华面积提取信息

Tab. 4 Bloom area extraction information

时间 像元数目/个 面积m2
21日10:30 1 354 627 109.72
21日14:00 1 493 579 120.97
21日17:30 1 588 194 128.64
22日10:30 1 030 591 83.47
22日14:00 1 186 588 96.11
22日17:30 1 298 194 105.15

4.3 水华预测结果

通过提取21日10:30—22日14:00共5个时刻的水华信息，结合水质参数及气象数据组成3组训练样本来构建水华预测训练模型，最终实现22日17:30的水华预测。在本文的预测模型中，预测窗口尺寸的选取直接影响预测结果的精度，分别采用8×8、10×10、12×12、14×14的预测窗口进行对比实验，遴选出最优尺寸窗口。图14为4种预测窗口的预测结果，可以看出当选用8×8的预测窗口时，预测模型生成的预测点最多，预测结果空间分布过于分散。当选用12×12、14×14的预测窗口时，预测点偏少，容易产生水华面积的过量预测，特别是在水华分布稀疏的区域预测效果不佳。相比其他3种尺寸的预测窗口，10×10预测窗口的预测点数目和预测半径大小都更加适宜，预测结果与提取结果在空间分布上也更为接近，尤其是在水体边界及下游水华聚集区的重合度很高，虽然水华分布分散的区域预测结果并不理想，但在一定程度上反映了水华的爆发趋势和空间分布，可为湖库水华爆发风险的短期预测提供参考。
图14 4种尺寸预测窗口的预测结果

Fig. 14 Prediction results of four sizes of prediction Windows

为了准确评估4种预测窗口预测结果的精度，以22日17:30提取的水华作为真实值，通过引入总体分类精度、生产者精度和Kappa系数对预测结果进行精度评价。由表5可看出，22日17：30提取的水体信息共30 834 683个像元，其中提取水华的像元共1 298 194个。使用4种预测窗口预测结果提取的水华像元个数分别为1 175 486、1 248 745、 1 021 482、847 268，其中，10×10预测窗口预测结果提取的水华像元数最接近真实值，水华面积预测误差仅为3.81%，表明其在水华总面积预测上有更高的准确性。图15为4种预测窗口水华预测正确像元和错误像元分布情况，可以看出10×10预测窗口正确预测像元的个数更多，特别是在下游水华聚集区更明显。4种预测窗口预测结果的3种精度指标如表6所示，综合来看，不同预测窗口预测结果整体精度的由高到低为10×10 > 8×8 > 12×12 > 14×14，4种预测窗口预测结果都具有较高的总体分类精度（0.97~0.98），表明预测模型对水华和非水华有着不错的辨别能力，生产者精度和Kappa系数均高于0.5，表明预测模型在水华的面积及空间分布上都有良好的预测能力。其中，10×10预测窗口预测结果与NDVI阈值提取结果最吻合，其总体分类精度、生产者精度、Kappa系数分别为0.98、0.77、0.77。
表5 水华提取像元个数（22日17:30）

Tab. 5 Number of pixels extracted from bloom (17:30 on 22nd)

数据 水华像元数/个 非水华像元数/个 水华面积预测误差/%
NDVI阈值提取结果 1 298 194 29 536 489 -
8×8预测窗口预测结果 1 175 486 29 659 197 9.45
10×10预测窗口预测结果 1 248 745 29 585 938 3.81
12×12预测窗口预测结果 1 021 482 29 813 201 21.32
14×14预测窗口预测结果 847 268 29 987 415 34.73
图15 水华预测正确像元和错误像元分布情况

Fig. 15 Distribution of correct and wrong pixels in the prediction of water bloom

表6 4种预测窗口预测结果的精度指标

Tab. 6 Precision index of prediction results of four prediction windows

预测窗口 预测类型 提取类型 生产者精度 总体分类精度 Kappa系数
水华 非水华
8×8 水华 860 661 314 825 0.66 0.98 0.68
非水华 437 533 29 221 664
10×10 水华 997 949 250 796 0.77 0.98 0.77
非水华 300 245 29 285 693
12×12 水华 708 751 312 731 0.55 0.97 0.59
非水华 589 443 29 223 758
14×14 水华 588 348 258 920 0.45 0.97 0.53
非水华 709 846 29 277 569
根据每个基本单元内初始时刻和预测时刻的水华面积，分别计算每个基本单元内水华的扩张或缩小面积，然后建立对应半径的缓冲区，以此表示每个基本单元水华的面积变化及位置变化。从 图16中可以看出，相比初始时刻，水华面积以扩张为主，主要分布在下游西侧的水华聚集区；水华面积缩小区域主要分布在上游、中游近岸水域。水华的位置变化整体上表现为由上游向下游，由东向西移动，大多聚集在下游水域，少量聚集在上游西侧近岸水域。
图16 水华面积变化（22日17:30）

Fig. 16 Change of cyanobacteria bloom area (17:30, 22nd)

5 讨论

本文基于水华生长机理和空间移动规律，借助遥感手段建立水华短时预测模型，实现了2022年夏季某月22日17:30鹤地水库石角站的水华预测。相比前人基于卫星遥感数据在深圳市水库[27]、鄱阳湖[28]、武汉东湖[29]的水质反演结果（表7），无人机遥感水质反演的监测范围具有局限性，但精度更高，时效性更强。且本文提出的模型不仅运用了主流的光学遥感方法[30-31]，还考虑了水质、水温、流速、风速、风向等对水华分布的影响，在鹤地水库取得了较好的水华预测效果，并且可以获取水华面积及位置变化的趋势，有助于确定水华暴发程度和时间，掌握水华暴发的规律。但预测模型还存在的不足之处，主要表现在以下几方面：
表7 本文与其他文献结果对比

Tab. 7 Existing water quality inversion research information

研究区 数据 反演精度（R²） 反演周期 文献来源
铁岗水库、西丽水库、石岩水库 GF-1 WFV影像 COD=0.83 1 年/次 [27]
鄱阳湖 GF-1 WFV影像 TP=0.74，TN=0.85 1 年/次 [28]
武汉东湖 Landsat7 ETM+SLC-OFF影像 TP=0.73，TN=0.71 1 月/2次 [29]
鹤地水库石角站 无人机遥感影像 TP=0.85，TN=0.89，COD=0.89 3.5 小时/次 本文
考虑的环境因素具有局限性。本预测模型仅考虑了pH、总氮、总磷、溶解氧、水温对水华爆发程度的影响及风速、流速、风向对其空间分布的影响，而水华的爆发是一个受水环境、气象、气候和生物行为等多重因素共同影响的复杂过程[32-33]，在以后的研究中应考虑更多的影响因素，来实现对水华预测模型的进一步优化，有助于加强对未来水华情势的精准估测。
缺乏对水华物种信息的考虑。本预测模型未考虑水华物种的差别对模型精度的影响，采用在一定程度上反映水华生物量的NDVI表征水华强度，而鹤地水库藻类群落结构的季节变化较大，空间差异性也很大，在今后模型的优化中，可将水华物种信息纳入预测模型中。

6 结论

本文使用无人机遥感影像和同步实测数据，基于MPP算法和GWR构建了水华短时预测模型，预测了水华的面积和空间分布，分析了预测窗口对预测结果的影响，结果表明：
（1）基于无人机多光谱遥感影像结合地面水质同步采样数据定量反演了总氮、总磷和溶解氧的浓度，3种水质参数反演结果与实测数据的决定系数分别达到0.89、0.85、0.89。NDVI阈值法提取的水华信息在水华爆发程度和空间分布上都与真彩色正射影像目视结果有着较高的一致性，可以反映出较真实的水华信息。
（2）本文提出的水华预测模型对水华面积的预测精度达到96.19%，能准确预测水华的空间分布，反映水华面积及位置变化特性，对水华和非水华的总体分类精度均大于0.96，生产者精度、Kappa系数均大于0.5。
（3）不同尺寸的预测窗口直接影响预测结果精度，相比其他3种尺寸，10×10预测窗口与无人机数据的空间分辨率更匹配，得到的水华面积预测误差最低（3.81%），产生的预测点数与实际更加吻合，在水华的空间分布上其总体分类精度、生产者精度、Kappa系数最高，分别为0.98、0.77、0.77，故最适合用于本文水华的预测。
本文提出的水华短时预测模型是耦合MPP算法和GWR，应用于水华的精细化监测、预测，实现了未来短时的水华面积及空间分布的模拟。通过对预测窗口的调整可适用于不同水域的水华预测，可为水华预警监测提供一定的技术支持，为相关部门提供借鉴与指导，提前制定合理的水华应急和防治措施。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
段恬媛. 湖库型饮用水水源地水安全评价及风险评估[D]. 四川: 西南交通大学, 2021.

[Duan T Y. Water safety assessment and risk assessment of lake-reservoir drinking water sources[D]. Chengdu: Southwest Jiaotong University, 2021.]

[2]
刘绍俊, 施艳峰, 翟竟余, 等. 星云湖微囊藻密度周年变化及其与水质指标的关系[J]. 环境化学, 2021, 40(7):2064-2072.

[Liu S J, Shi Y F, Zhai J Y, et al. Annual variations of Microcystis density and their relationships with water quality indices in Xingyun Lake[J]. Environmental Chemistry, 2021, 40(7):2064-2072.] 10.7524/j.issn.0254-6108.2020030401

DOI

[3]
史小丽, 杨瑾晟, 陈开宁, 等. 湖泊蓝藻水华防控方法综述[J]. 湖泊科学, 2022, 34(2):349-375.

[Shi X L, Yang J S, Chen K N, et al. Review on the control and mitigation strategies of lake cyanobacterial blooms[J]. Journal of Lake Sciences, 2022, 34(2):349-375.] DOI:10.18307/2022.0201

DOI

[4]
姚金忠, 范向军, 杨霞, 等. 三峡库区重点支流水华现状、成因及防控对策[J]. 环境工程学报, 2022, 16(6):2041-2048.

[Yao J Z, Fan X J, Yang X, et al. Current situation, causes and control measures of water bloom in the key tributaries of the Three Gorges Reservoir[J]. Chinese Journal of Environmental Engineering, 2022, 16(6):2041-2048.] DOI:10.12030/j.cjee.202112199

DOI

[5]
Yu M L, Hong G X, Xu H, et al. Effects of cyanobacterial blooms in eutrophic lakes on water quality of connected rivers[J]. Huan Jing Ke Xue, 2019, 40(2):603-613.] DOI:10.13227/j.hjkx.201804047

DOI

[6]
Rubén D M, Vojtech B, Antonio O, et al. Evaluation of daily and seasonal variations in a semi-closed photobioreactor for microalgae-based bioremediation of agricultural runoff at full-scale[J]. Algal Research, 2020, 47(1):2211-9264. DOI:10.1016/j.algal.2020.101859

DOI

[7]
秦伯强. 浅水湖泊湖沼学与太湖富营养化控制研究[J]. 湖泊科学, 2020, 32(5):1229-1243.

[Qin B Q. Shallow Lake limnology and control of eutrophication in Lake Taihu[J]. Journal of Lake Sciences, 2020, 32(5):1229-1243.] DOI:10.18307/2020.0501

DOI

[8]
Chaffin J D, Kane D D, Johnson A. Effectiveness of a fixed-depth sensor deployed from a buoy to estimate water-column cyanobacterial biomass depends on wind speed[J]. Journal of Environmental Sciences, 2020, 93(7):23-29. DOI:10.1016/j.jes.2020.03.003

DOI

[9]
Shi K, Zhang Y L, Qin B Q, et al. Remote sensing of cyanobacterial blooms in inland waters: Present knowledge and future challenges[J]. Science Bulletin, 2019, 64(20):1540-1556. DOI:10.1016/j.scib.2019.07.002

DOI PMID

[10]
嵇晓燕, 张迪, 李文攀, 等. 湖泊蓝藻天地一体化监测业务化运行平台的构建——以滇池流域为例[J]. 环境保护, 2016, 44(18):24-27.

[Ji X Y, Zhang D, Li W P, et al. Operational platform construction of heaven and earth integration algae monitoring in Lake Dianchi[J]. Environmental Protection, 2016, 44(18):24-27.] DOI:10.14026/j.cnki.0253-9705.2016.18.005

DOI

[11]
邱银国, 段洪涛, 万能胜, 等. 巢湖蓝藻水华监测预警与模拟分析平台设计与实践[J]. 湖泊科学, 2022, 34(1):38-48.

[Qiu Y G, Duan H T, Wan N S, et al. Design and practice of a platform for monitoring, early-warning and simulation of algal blooms in Lake Chaohu[J]. Journal of Lake Sciences, 2022, 34(1):38-48.] DOI:10.18307/2022.0102

DOI

[12]
段洪涛, 曹志刚, 沈明, 等. 湖泊遥感研究进展与展望[J]. 遥感学报, 2022, 26(1):3-18.

[Duan H T, Cao Z G, Shen M, et al. Review of lake remote sensing research[J]. National Remote Sensing Bulletin, 2022, 26(1):3-18.] DOI:10.11834/jrs.20221301

DOI

[13]
张运林, 张毅博, 李娜, 等. 利用陆基高光谱遥感捕捉太湖蓝藻水华日内快速变化过程[J]. 湖泊科学, 2021, 33(6):1951-1960.

[Zhang Y L, Zhang Y B, Li N, et al. Capturing the rapid intra-day change of cyanobacteria bloom by land-based hyperspectral remote sensing in Lake Taihu[J]. Journal of Lake Sciences, 2021, 38(6):1951-1960.] DOI:10.18307/2021.0628

DOI

[14]
张虎军, 宋挺, 朱冰川, 等. 太湖蓝藻水华暴发程度年度预测[J]. 中国环境监测, 2022, 38(1):157-164.

[Zhang H J, Song T, Zhu B C, et al. Annual forecast of the extent of cyanobacteria bloom in Taihu Lake[J]. Environmental Monitoring in China, 2022, 38(1):157-164.] DOI:10.19316/j.issn.1002-6002.2022.01.15

DOI

[15]
吴娟, 朱跃龙, 金松, 等. 三种机器学习模型在太湖藻华面积预测中的应用[J]. 河海大学学报(自然科学版), 2020, 48(6):542-551.

[Wu J, Zhu Y L, Jin S, et al. Area prediction of cyanobacterial blooms based on three machine learning methods in Taihu Lake[J]. Journal of Hohai University (Natural Sciences), 2020, 48(6):542-551.] DOI:10.3876/j.issn.10001980.2020.06.009

DOI

[16]
祁国华, 马晓双, 何诗瑜, 等. 基于多源遥感数据的巢湖水华长时序时空变化(2009—2018年)分析与发生概率预测[J]. 湖泊科学, 2021, 33(2):414-427.

[Qi G H, Ma X S, He S Y, et al. Long-term spatiotemporal variation analysis and probability prediction of algal blooms in Lake Chaohu (2009-2018) based on multi- source remote sensing data[J]. Journal of Lake Sciences, 2021, 33(2):414-427.] DOI:10.18307/2021.0204

DOI

[17]
刘东, 张民, 沈明. 基于遥感藻总量和气象因子的巢湖不同湖区藻华预测[J]. 湖泊科学, 2022, 34(4):1090-1101.

[Liu D, Zhang M, Shen M, et al. Predicting algal blooms in different zones of Lake Chaohu based on satellite-derived algal biomass and meteorological factors[J]. Journal of Lake Sciences, 2022, 34(4):1090-1101.] DOI:10.1830 7/2022.0404

DOI

[18]
Zhao H, Li J S, Yan X, et al. Monitoring cyanobacteria bloom in Dianchi Lake based on ground-based multispectral remote-sensing imaging: Preliminary results[J]. Remote Sensing, 2021, 13(19),3970. DOI:10.3390/rs13193970

DOI

[19]
丁铭, 李旭文, 姜晟. 基于无人机高光谱遥感在太湖蓝藻水华监测中的一次应用[J]. 环境监测管理与技术, 2022, 34(1):49-51,71.

[Ding M, Li X W, Jiang S, et al. An application of UAV-based hyperspectral remote sensing in monitoring cyanobacterial bloom in Taihu Lake[J]. The Administration and Technique of Environmental Monitoring, 2022, 34(1):49-51,71.] DOI:10.19501/j.cnki.1006-2009.2022.01.007

DOI

[20]
Lu B B, Hu Y G, Murakami D, et al. High-performance solutions of geographically weighted regression in R[J]. Geo-spatial Information Science, 2022, 25(4):536-549. DOI:10.1080/10095020.2022.2064244

DOI

[21]
Ying H T, Xia K, Huang X X, et al. Evaluation of water quality based on UAV images and the IMP-MPP algorithm[J]. Ecological Informatics, 2021, 61:101239. DOI:10.1016//j.ecoinf.2021.101239

DOI

[22]
潘鑫, 杨子, 杨英宝, 等. 基于高分六号数据提取太湖蓝藻方法的对比及适用性分析[J]. 湖泊科学, 2022, 34(6):1866-1876.

[Pan X, Yang Z, Yang Y B, et al. Comparison and applicability analysis of methods for extracting cyanobacteria from Lake Taihu based on GF-6 data[J]. Journal of Lake Sciences, 2022, 34(6):1866-1876.] DOI:10.18307/2022.0607

DOI

[23]
窦世卿, 张寒博, 徐勇, 等. TRMM降水数据在长江流域的降尺度分析与校正[J]. 中国农业气象, 2021, 42(5):377-389.

[Dou S Q, Zhang H B, Xu Y, et al. Research on downscaling and correction of TRMM cata in the Yangtze River Basin[J]. Chinese Journal of Agrometeorology, 2021, 42(5):377-389.] DOI:10.3969/j.issn.1000-6362.2021.05.003

DOI

[24]
覃文忠. 地理加权回归基本理论与应用研究[D]. 上海: 同济大学, 2007.

[Qin W Z. The basic theoretics and application research on geographically weighted regression[D]. Shanghai: Tongji University, 2007.]

[25]
叶章熙, 郭倩, 张健, 等. 基于无人机可见光影像与OBIA-RF算法的城市不透水面提取[J]. 农业工程学报, 2022, 38(4):225-234.

[Ye Z X, Guo Q, Zhang J, et al. Extraction of urban impervious surface based on the visible images of UAV and OBIA-RF algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering, 2022, 38(4):225-234.] DOI:10.11975/j.issn.1002-6819.2022.04.026

DOI

[26]
赵晶. 水库悬浮物遥感模型及其时空动态Markov评价[D]. 昆明: 昆明理工大学, 2020.

[Zhao J. Remote sensing model of reservoir total suspended solids and its Markov evaluation in spatiotemporal variation[D]. Kunming: Kunming University of Science and Technology, 2020.]

[27]
李俊, 张文志, 邓孺孺, 等. 高分一号WFV影像的深圳市水库CODMn浓度时空分布特征[J]. 遥感学报, 2022, 26(8):1562-1574.

[Li J, Zhang W Z, Deng R R, et al. Study of spatial—Temporal characteristics for CODMn in Shenzhen Reservoir based on GF-1 WFV[J]. National Remote Sensing Bulletin, 2022, 26(8):1562-1574.] DOI:10.11834/jrs.20219380

DOI

[28]
李怡静, 孙晓敏, 郭玉银, 等. 基于梯度提升决策树算法的鄱阳湖水环境参数遥感反演[J]. 航天返回与遥感, 2020, 41(6):90-102.

[Li Y J, Sun X M, Guo Y Y, et al. Remote sensing retrieval of water quality parameters in Poyang Lake based on the Gradient Boosting Decision Tree algorithm[J]. Spacecraft Recovery & Remote Sensing, 2020, 41(6):90-102.] DOI:10.3969/j.issn.1009-8518.2020.06.009

DOI

[29]
何报寅, 丁超, 杨小琴, 等. Landsat7 ETM+ SLC-OFF数据的修复及其在武汉东湖水质反演中的应用[J]. 长江流域资源与环境, 2011, 20(1):90-95.

[He B Y, Ding C, Yang X Q, et al. Recovering of Landsat7 ETM+ SLC-OFF data and its application on water quality retrieval in East Lake in Wuhan[J]. Resources and Environment in the Yangtze Basin, 2011, 20(1):90-95.] DOI:11870/cjlyzyyhj2011009006

DOI

[30]
晁明灿, 赵强, 杨铁利, 等. 基于Landsat8影像的蓝藻水华提取方法对比研究[J]. 大气与环境光学学报, 2021, 16(6):520-528.

[Chao M C, Zhao Q, Yang T L, et al. Comparative research of cyanobacteria blooms extraction methods based on Landsat8 images[J]. Journal of Atmospheric and Environmental Optics, 2021, 16(6):520-528.] DOI:10.3969/j.issn.1673-6141.2021.06.006

DOI

[31]
李一民, 谭振宇, 杨辰, 等. 基于多源卫星的滇池藻华提取机器学习算法研究[J]. 地球科学进展, 2022, 37(11):1141-1156.

DOI

[Li Y M, Tan Z Y, Yang C, et al. Extraction of algal blooms in Dianchi Lake based on multi-source satellites using machine learning algorithms[J]. Advances in Earth Science, 2022, 37(11):1141-1156.] DOI:10.11867/j.issn.1001-8166.2022.064.

DOI

[32]
杨素, 万荣荣, 李冰. 太湖流域水文连通性:现状、研究进展与未来挑战[J]. 湖泊科学, 2022, 34(4):1055-1074.

[Yang S, Wan R R, Li B. Hydrological connectivity research in Lake Taihu Basin: Status, progress and future chal-lenges[J]. Journal of Lake Sciences, 2022, 34(4):1055-1074.] DOI:10.18307/2022.0402

DOI

[33]
Xiang C H, Li Y, Ke Z X, et al. Effects of daily nitrogen and phosphorus input on planktonic community metabolism in a semi-enclosed bay by mesocosm experiment[J]. Acta Oceanologica Sinica, 2022, 41(8):99-110. DOI:10.1007/s13131-022-1986-4

DOI

Options
文章导航

/

〈 〉