江心洲是分汊河段中出露于水面的成型淤积体^[1]。作为长江中下游河流生态系统的重要组成部分，其具备净化水质、提供栖息地、维持生物多样性等重要功能^[2-4]。但在三峡工程运用的影响下，长江中游水沙情势剧变，江心洲生境结构发生变化，其上植被密度、活力及分布模式也随之改变^[5]。因此，有必要开展针对江心洲植物的长期观测，以制定合理的洲滩植被保护和修复方案^[6]。

相较于传统实地观测，卫星遥感技术可以大范围、周期性、低成本地对植被进行动态监测。在众多对地观测卫星中，Landsat系列卫星凭借着其长达40年的连续观测数据和30 m的空间分辨率，成为目前获取长时间、长河段植被动态的主要手段^[7]。大量研究利用Landsat数据解译了河口三角洲或湖泊型湿地中外来物种入侵或植物演替的过程^[8-10]。但对于物种丰富度相对较低且处于地理隔离的江心洲而言，厘清三峡工程运用对其上植物密度和活力的影响则更为关键。Zhou等^[11]应用Landsat数据评估了三峡运用前后长江中游江心洲植被面积的变化趋势；Zhu等^[12]则利用Landsat数据反演了荆江段洲滩的植被覆盖度，并进一步分析了植被覆盖度与淹没历时之间的关系。上述研究在方法层面存在2个问题：(1)其皆是通过比较前后两年汛期结束后的遥感影像以获取植被年际间变化，然而植物历年生长周期间的差异以及遥感影像中的非系统性误差将严重影响结果的可靠性；(2)其所提取的植被指标，包括植被指数和植被覆盖度等，无法反映植物一个生长期内活力、密度和长势等植被动态特征。而从植被指数序列中提取的植物物候指标则常用于反映年植被动态，也适用于分析植物年际间变化趋势^[13]。但就常被洪水淹没的江心洲而言，淹没前后的植被指数序列中含有不连续的植物物候信息。因此，若要构建适用于江心洲植物的物候指标，需基于精确的淹没数据，并结合植物特性提出适当的物候提取方法。通常情况下，用于获取湿地淹没范围的手段包括水动力学模型、原型观测和卫星遥测^[14]。鉴于前两者皆需开展大量实地勘察，卫星遥测依旧是目前最为高效的手段。利用卫星遥测数据反演洪水淹没范围的方法可分为两类：基于大量实测数据的监督分类方法和基于特定规则的非监督分类方法^[15]。目前诸多研究结合植被指数、水体指数^[16-17]或合成孔径雷达数据^[18]，构建了可用于识别湿地淹没情况的监督分类模型。但鉴于江心洲实测淹没数据较为匮乏，基于长江流域分界明显的汛期与非汛期，构建一种非监督变化检测法以替代普适性较差的静态阈值法成为一种新的途径。

本文选取长江中游4个典型江心洲作为研究对象，利用Landsat数据构建一系列江心洲淹没范围和植被动态反演方法，并定量分析：(1)近40年江心洲植被动态变化趋势；(2)三峡工程运用前后漫滩洪水对植被长势的影响；(3)三峡工程运用前后不同高程梯度下植被分布模式的调整规律。本文所提出的方法及得到的结论可为监测长江中游江心洲植被动态和制定相关修复策略提供支撑。

1 研究区域与数据 1.1 研究区域概况

长江中游是指宜昌至湖口之间的河段，该河段内有近42个江心洲，枯水期出露面积可达400 km^{2 [19]}。考虑到大部分江心洲处于非自然状态，本文选取城汉河段(城陵矶-汉口河段)4个洲体形态稳定且受人类活动影响较小的江心洲作为研究对象，包括南阳洲、南门洲、中洲和白沙洲(图 1a~d)。这4个江心洲位于三峡大坝下游447~542 km处，枯水期出露面积约为6~9 km²。实地调查和以往文献资料表明：4个江心洲可根据植物种类和多年淹没频率在垂向上划分为边滩、低滩和高滩3个部分(图 1e~f)^[20-21]。其中，边滩中水位时淹没，枯水位时露出水面，少量耐水植物会在退水后重新萌发；低滩在汛期会随着水位波动而呈现不同的淹没情况，其上分布有耐水性较强的虉草、狗牙根等禾本科植物；而高滩只在大洪水时才会被淹没，其上分布着耐水性较差的芦苇、南荻、白茅和少量灌木、乔木。

1.2 Landsat卫星遥感数据及预处理

Landsat系列卫星是美国航空航天局自1972年发射的一系列观测地球表面资源与环境的卫星。其具有16 d的单星重访周期(8 d的双星重访周期)和30 m的空间分辨率(可见光、近红外和中红外波段)。本文从美国地质调查局(USGS)地球勘测数据库(https://earthexplorer.usgs.gov/)中检索并下载了1986—2020年覆盖4个江心洲且具有最高质量(Tier 1)的所有Landsat二级(Level-2)产品，共计780幅(表 1)。对于扫描矫正器(Scan Lines Corrector, SLC)失效的Landsat 7 ETM+数据，使用原始SLC-off图像以避免数据填充对后续物候指标提取的影响。Landsat二级产品中包括可见光和红外波段的地表反射率以及质量评估(Quality Assessment, QA)波段。QA波段主要为CFMask算法生成的云层数据^[22]，其用二进制的方式存储有云、云阴影的分布情况。依据Landsat 4-7以及Landsat 8官方产品指南，利用开源程序(https://github.com/shixi2ng/Satellite_RS_process_module)对Landsat影像进行去云处理。

1.3 谷歌地球高分辨率遥感影像

谷歌地球是一款融合多源遥感数据的地图软件。本文收集了谷歌地球中汛期(5—10月)高分辨率遥感影像，用于评估淹没提取方法的精度。但高分辨率卫星影像重访周期长且受云层遮挡影响，很难找到与Landsat同一天探测的数据。因此，依据以下2点要求对数据进行初步筛选：(1)所选高分辨遥感影像前后各16 d内存在不受云层影响的Landsat数据；(2)Landsat与高分辨率影像观测日期之间的绝对水位差不超过1 m。经筛选，共计11幅高分辨率影像满足要求(表 2)，但其中仅有1幅(2020-08-16/南门洲)捕捉到了高低滩植物群落被洪水淹没的情况，而其余影像只拍摄到边滩淹没情况。随后，利用eCognition^®商业软件中的多尺度分割工具对谷歌高分辨率遥感影像进行处理，并手动赋予分割结果不同的淹没情况，以此得到淹没样本数据。

1.4 水沙数据

本文选用长江中游莲花塘、石矶头水位站和螺山、汉口水文站的日均水位数据，4个站点的地理位置详见图 1。随后沿干流利用线性内插的方法计算出每个江心洲中点处的水位。同时，收集了螺山站的含沙量数据用以分析该河段含沙量变化。

2 研究方法 2.1 植被和水体指数序列的构建

由于江心洲土壤裸露度和含水量随季节波动较大，本文选用优化型土壤调节植被指数(OSAVI)^[23]以减轻土壤条件波动对植被指数序列的影响。另一方面，构建了改进的归一化水体指数(MNDWI)^[24]序列用于识别淹没区域：

$ O S A V I=1.16 \times \frac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{Red}}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}}+0.16} $

(1)

$ M N D W I=\frac{\rho_{\text {Green }}-\rho_{\text {MIR }}}{\rho_{\text {Gren }}+\rho_{\text {MIR }}} $

(2)

式中，ρ_NIR、ρ_Red、ρ_Green和ρ_MIR分别代表近红外、红光、绿光和中红外波段的地表反射率；常数0.16为优化的土壤调节参数。

2.2 淹没区域的识别

由于缺乏足够的实测资料，本文基于长江中游流域界线明显的汛期和非汛期及江心洲较稳定的地貌类型构建了一种非监督变化检测法(图 2)。该方法的基本原理是在像素尺度构建MNDWI动态阈值进而判别淹没状态。首先，利用Zhu等^[25]提出的连续变化检测和分类(continuous change detection and classification，CCDC)算法处理所有Landsat数据，生成1986—2020年江心洲地貌类型的变化情况。经分析后发现，江心洲大部分区域未发生明显地貌变化，仅有少量处于低滩和边滩交界的区域在1995—2006年期间发生过一次转变。随后，针对每个像素，构建其处于同种地貌类型下的MNDWI序列直方图，对于低滩和少部分高滩区域，其MNDWI序列的直方图呈现明显的双峰分布(图 2)，基于这一特征，先利用Zhang等^[26]提出的改进的双峰判别法确定双峰直方图的两峰分界线T_inti。随后将MNDWI序列划为汛期(5—10月)和非汛期(上年11月-当年5月)两组，非汛期序列中小于T_inti的数据服从高斯或偏态分布，因此利用下述公式计算出该像素当前地貌类型下的动态阈值MNDWI_i:

$ \mu_{\mathrm{nf}_i}=\operatorname{mean}\left(M N D W I_{\mathrm{nf}_i}\right) \quad M N D W I_{\mathrm{nf}_i}>T_{\mathrm{int}_i} $

(3)

$ \delta_{\mathrm{nf}_i}=\operatorname{std}\left(M N D W I_{\mathrm{nf}_i}\right) \quad M N D W I_{\mathrm{nf}_i}>T_{\text {int }_i} $

(4)

$ M N D W I_i=\mu_{\mathrm{nf}_i}+\Delta \cdot \delta_{\mathrm{nf}_i} $

(5)

式中，下标i代表该像素的第i种地貌类型；μ_nfi、δ_nfi分别代表该像素非汛期第i种地貌类型下MNDWI序列的平均值和标准差；Δ为标准差的系数，根据3δ法则选定为3^[27]；MNDWI_i为该像素第i种地貌类型下的MNDWI动态阈值，在MNDWI汛期序列中超过该动态阈值的数据被认定为淹没。而对于永久性水体和淹没频率较低的高滩，其MNDWI序列呈偏态分布，因此直接利用所有非汛期数据即可构建动态阈值MNDWI_i。应用上述方法判别所有影像的淹没范围后，将同一年内的淹没分布图合成为年淹没分布图，并利用多年的年淹没分布图计算出年淹没频率。

2.3 江心洲高低滩的划分

江心洲包括边滩、低滩与高滩3个部分，由于江心洲实测地形数据较为稀缺且空间分辨率较低，本文利用江心洲年淹没频率以及日水位数据划分江心洲高滩、低滩与边滩。首先，利用平滩水位划分边滩与低滩，主要步骤包括：(1)选取研究时段退水期内(长江中游为10—12月)日水位处于螺山站平滩水位±0.5 m以内(22.93±0.5)m的所有遥感影像；(2)利用前文中提出的非监督变化检测法识别所有影像中未淹没的区域；(3)计算出整个滩体与临近河道的未淹没频率，将未淹没频率大于95%的像元视为高低滩区域，其余像元则为边滩或河道。另一方面，高低滩的划分则主要依据年淹没频率方法为：(1)利用螺山站水位数据分析所选年份内极端洪水的出现频率；(2)将年淹没频率小于年极端洪水频率的区域视为高滩，其余则为低滩。

2.4 植物物候指标的提取

在利用非监督变化检测法剔除植被指数序列中受淹没影响的数据后，进一步应用Elmore等^[28]提出的七参数双逻辑曲线对植被指数序列进行平滑处理：

$ v(m, t)=m_1+\left(m_2-m_7 t\right)\left(\frac{1}{1+\mathrm{e}^{\frac{m_3-t}{m_4}}}-\frac{1}{1+\mathrm{e}^{\frac{m_5-t}{m_6}}}\right) $

(6)

式中，t为年积日；m₁表征冬季植被指数均值；m₂表征夏季植被指数峰值与 m₁间的差值；m₃与m₅分别表示植物开始生长与完全枯死的年积日；m₄和m₆分别表征夏季植物生长速率与秋冬季衰退速率；m₇则表征夏季衰退速率。利用七参数双逻辑公式处理后的平滑曲线可划分为4个阶段：快速生长期(3—5月)、成熟期/缓慢衰退期(5—10月)、快速衰退期(10—11月)和休眠期(12月-次年3月)(图 3a)。根据Karunaratne等^[29]的研究，这种两段式的衰退主要与禾本科植物茎和叶不同的生长趋势有关，叶生物量在快速生长期结束后缓慢下降，而茎生物量在秋冬季才开始减少。基于这种植物特性，本文选用年植被指数峰值与成熟期平均植被指数2个物候指标来综合反映江心洲植被活力与茂密程度的变化(后文中用植被动态代替)。年植被指数峰值定义为快速生长期结束、成熟期开始这段时间内植被指数的最大值，常用于表征植被动态。然而，受制于Landsat 16 d的单星重访周期以及长江中游流域的梅雨季节，大部分情况下Landsat并不能捕捉到真实情况下植被指数的最大值^[13]，因此该指标只能用于分析植被动态长期变化趋势或空间分布规律。而另一指标，成熟期平均植被指数定义为年植被指数峰值至快速衰退期开始这段时间内植被指数的平均值，如若当年该区域被洪水淹没则将范围缩短到峰值至洪水事件开始这一段时间。此外，在应用成熟期平均植被指数分析年际间植被动态时，如果两年峰值至洪水事件开始这段时间范围长度不等，则只比较两年中较短范围内的植被指数均值。

3 结果分析与讨论 3.1 淹没识别方法精度评价

引入动态地表水范围(dynamic surface water extent, DSWE)^[30]和全球地表水(global surface water, GSW)^[31]两种水体分布产品，对本文提出的非监督变化检测进行精度评价，样本集选用高分辨率遥感数据生成的淹没范围。首先，统计了各种淹没识别方法的分类精度(表 3)，结果表明非监督变化检测法拥有最高的总体分类精度且表现较为稳定，4个江心洲的总体分类精度在94.4%~98.9%之间。而其他两个基于多光谱遥感数据的识别方法，DSWE和GSW在南门洲表现出了明显误差，总体分类精度相较非监督变化检测法分别降低了5.8%和3.9%，其中差异主要来源于捕捉到漫滩洪水淹没植被的一景影像。

由于评价淹没识别方法精度的一个重要依据为其能否正确分辨植被的淹没情况。因此进一步生成了3种识别方法绘制的淹没范围(图 4)。图 4a中蓝色的区域代表利用高分辨率遥感数据生成的淹没范围(已被重采样到Landsat 30 m网格中)。从图中挑选了2个具有代表性的区域：区域α和β。区域α主要由少量植物冠层和大面积淹没区域构成；区域β则包含4种地物类型：裸土、未淹没植被、浅灰色淹没区和水绿色淹没植被(图 4a副图)。在30 m分辨率的Landsat影像中，由于存在光谱混合，区域α内像元的相关指标可能更偏向于未淹没的情况从而导致误判，例如GSW(图 4d)。然而，对于非监督变化检测法而言，区域α中像元的动态MNDWI阈值主要集中在-0.11~-0.08这一范围内，当该区域被淹没时，MNDWI会立即提高至-0.04~0.13，从而精准分类。对于区域β，DSWE将未淹没或是土壤含水量较高的植物像元错误地识别为淹没，而GSW将区域β右下角的大部分洼地错误地识别为未淹没区域。相比之下，非监督变化检测法精准地划定了水陆边缘。另一方面，分析了式(5)中Δ对识别精度的影响，当Δ∈[2, 3]时，总体识别精度变化较小，而错分和漏分误差相互转换，因此Δ的取值对整体识别精度影响不大。综上所述，本文所提出的非监督变化检测法相比目前主流水体分布产品表现得更为精确且稳定，适合于分析长江中下游江心洲淹没动态。

3.2 近40年江心洲植被动态与淹没频率变化趋势

图 5给出了近40年4个江心洲上年植被指数峰值的分布情况及变化规律。图中每个箱形图的最大、最小值分别代表年植被指数峰值的5%和95%分位值，并用其进一步构建了上、下包络线。此外，利用历年江心洲年植被指数峰值的中位数，分别在三峡运行前、运行后以及全部时段3个范围内进行线性回归分析，以回归分析所得斜率k表征多年植被动态变化趋势。结果表明：(1)近40年内，4个江心洲的植被动态皆呈现上升趋势(k>0)，中洲和南阳洲在三峡运行后的上升趋势大于三峡运行前，南门洲的植被动态则主要在三峡运行前有明显的上涨(k_三峡前=0.0091)，而三峡运行后保持稳定(k_三峡后=-0.0020)；(2)南门洲、南阳洲以及中洲上植被指数最小值分别在1996、2006、2018年之后发生明显的抬升，从接近于0的水平抬升至0.3~0.4；(3)白沙洲长时间植被指数在三峡修建前后都稳定在0.4~0.6，近40年仅有小幅度的增长(k =0.0018)。其中，3个江心洲植被指数最小值发生抬升表明江心洲植被覆盖范围有所增长，这种现象往往与江心洲淹没频率和范围的变化相关。

因此，进一步生成了4个江心洲三峡运行前后年淹没频率的分布情况(图 6)。对比三峡运行前后淹没频率分布图可以发现，三峡运行后南阳洲与中洲低滩出现明显的扩张。综合其植被动态变化进行分析：在河道局部性淤积以及三峡大坝削弱洪峰的共同作用下，部分江心洲原有边滩淹没频率下降并逐渐满足耐水禾本科植物的生长条件，促使江心洲植被动态最低值发生抬升。而另一方面，白沙洲和南门洲2个江心洲低滩和高滩的范围没有出现明显的调整，但三峡运行后高、低滩淹没频率相较运行前有明显的下降，这会导致原有高滩上不耐水但长势茂密的禾本科植物逐渐向这些区域蔓延，促使江心洲植被动态整体上的小幅抬升。

3.3 江心洲植被对三峡运行前后洪水的响应过程

以2002年与2016年两场漫滩洪水为例，分析了江心洲植被对三峡运行前后漫滩洪水的响应过程。这两场漫滩洪水具有相似的流量过程(双样本K-S检验下，D_mn=0.034 m，P＜0.05)，且前序1~2年皆没有其它洪水漫滩，保证植物根系有充足的时间从土壤中获取必要的营养物质。随后利用年际成熟期平均植被指数的变化量表征植被动态变化，如图 7所示。图 7a~b表明白沙洲2002年被淹没区域在2003年呈现出长势更好的绿色，而未被淹没区域则呈现长势变差的红色；图 7c~d、g~h表明2016年白沙洲和南门洲被淹没区域的植被动态在次年呈现出减弱的趋势，而未被淹没的区域则呈现上涨趋势。初步结果表明：相较2002年洪水，2016年的洪水事件对江心洲植物的生长更具不利影响，这种现象可能与气候条件、漫滩洪水间及植物群落的差异有关。

为剔除气候条件对植被的影响，进一步统计了2场洪水前后江心洲淹没和未淹没区域年际间成熟期平均植被指数变化量的分布情况，如图 8所示。结果表明，2002年白沙洲和南门洲淹没区域的变化量相较未淹没区域存在整体的抬升，平均增长量分别为0.023和0.039；而2016年淹没区域相较未淹没区域呈现下降的态势，平均减少量分别为-0.018和-0.009。由于两场洪水具有相似的规模，同时气候条件对植物的影响也已经通过比较同年淹没和未淹没区域间的差异而抵消，这一现象可能与河流和江心洲间能量、物质交换相关。理论上，漫滩洪水很难正向促进高滩不耐水禾本科植物的生长，而洪水中所携带的泥沙则有可能通过改变植物生境影响植物长势^[32]。过往大量研究表明生源物质含量越高，禾本科植物的密度、高度、生物量则越大^[33-34]。分析螺山站的含沙量数据表明，在水库拦沙的作用下，2016年洪水的日均含沙量为0.08 kg/m³，相较2002年的0.39 kg/m³下降了80%，而吸附在泥沙上的磷、氮也随之减少^[35]，导致输送至江心洲的生源物质大幅降低。因此，三峡工程运行前，含沙量较大的洪水通过向江心洲输移更多生源物质促进植物生长；而运行后，含沙量较小的洪水则有抑制植物生长的作用。

3.4 不同高程梯度下的植被动态变化趋势

为研究不同高程梯度下植被动态变化，分别计算了2000—2020年高滩和低滩植被指数峰值分布的平均值，以表征高、低滩植被动态(图 9)。此外，图 9给出了上一年各江心洲的年最高水位柱状图，其中橘色代表上一年江心洲遭遇大洪水。结果表明：对于形态较为稳定的白沙洲和南门洲，高、低滩植被动态具有相似的变化规律：在三峡工程运行之前，白沙洲和南门洲高滩较低滩植被动态分别高0.03和0.06；三峡工程运行后，高低滩差异逐渐缩小；三峡工程运行10年后，高低滩差异完全消失，江心洲植被分布模式产生改变。将江心洲视为一个完整的系统，三峡工程运行下水文节律的变化超过了原有系统的承受能力，迫使白沙洲和南门洲从三峡运行前的稳态逐渐转变为运行后的新稳态。而对于南阳洲，三峡运行后其洲体形态出现明显的调整。因此，其高、低滩植被动态在三峡运行后伴随地形出现明显的调整。

此外，进一步构建了白沙洲高、低滩植物动态分布柱状图并对其进行了归一化处理，柱状图的重叠率用于表征高、低滩植被分布模式的同质性，如图 10所示。这8年涵盖了白沙洲经历的5个阶段：三峡运行前旧稳态、转换阶段、三峡运行后新稳态、三峡运行后新稳态对极端洪水的响应、人类活动影响下的江心洲。结果表明：高、低滩植被分布模式的同质性从2002年的66.5%上升至2013年的93.6%(图 10a~e)；高、低滩植被分布模式的空间异质性在2013年时已基本消失(图 10e)；2016年漫滩洪水致使低滩植被动态分布发生了明显的偏移，高低滩植被动态分布均值间的差异从2013年的0.005升至0.028，大于三峡修建前2002年的高低滩差异(0.09)，这一现象表明三峡运行后新稳态对极端洪水更为敏感(图 10f)；2020年高、低滩植被动态分布模式发生了较为明显的改变(图 10h)。对比分析2013年6月13日与2020年11月9日的高分辨率遥感影像后发现，白沙洲在这段时间内遭受了人为开垦(图 10i)，进一步加剧了江心洲植被的退化。

4 结论

受三峡工程运用的影响，长江中下游洲滩植被活力、密度和分布模式发生变化，亟需开展持续观测以揭示其对变化水沙条件的响应机制。本文利用长时间Landsat遥感数据，构建了一系列江心洲淹没范围和植被动态解译方法，并定量分析江心洲植被动态变化趋势，主要结论如下：

1) 三峡工程运行后，江心洲整体淹没频率降低，原有边滩逐步被植物覆盖，促使江心洲植被覆盖面积增加。同时，高滩不耐水植物向低滩蔓延，导致江心洲植被茂密程度呈现增长的趋势。

2) 三峡工程运行前，含沙量较大的漫滩洪水有促进植物第二年生长的趋势；而三峡工程运行后具有同样规模但含沙量减少近80%的漫滩洪水则有抑制植物生长的作用。

3) 部分形态稳定的江心洲，其高、低滩植被分布模式之间的异质性在三峡工程运用后逐渐减弱，并在约10年后彻底消失。这部分江心洲遭受严重的功能性退化，并对外界干扰表现得更为敏感。