湿地水文过程作为影响湿地生态系统的重要因素之一，不仅左右着湿地的物理、化学和生态过程，也对湿地发育演化和维持景观效益起到重要作用^1]，鄱阳湖水文过程同时受到赣、抚、信、饶、修五水及长江的影响^[2].按照鄱阳湖水利枢纽建设办公室提供的观测数据，自1950年以来，鄱阳湖的水位(星子水位)在7~23m之间波动，最高水位为22.52m(1998年8月2日)，最低水位为7.11m(2004年2月4日).涉及鄱阳湖湖泊水文节律的研究主要使用定位观测数据^{[3 — 12]}和遥感、GIS^{[12 — 14]}等技术.闵骞等自1990s开始，持续关注鄱阳湖观测水位数据，从多角度描述鄱阳湖水位变化特征和规律，并在此基础上讨论了围垦等人类活动对鄱阳湖水位的影响^{[9 — 11]}.李鹏等利用1989-2010年多期卫星影像信息提取湖泊面积，结合定位观测数据得到了鄱阳湖水面面积随水位变化的关系^[12].针对鄱阳湖地区的湿地植被，研究方法集中在采用传统的采样调查^{[15 — 21]}、使用遥感数据反演生物量^[22]、解译植被种类^{[19 — 22]}等.张全军等通过实地调查研究了南矶自然保护区的植被类型及分布^[18].胡振鹏等利用3S技术对鄱阳湖湿地进行分类，研究了主要植物群落结构和分布规律^[19].葛刚等通过长期实地观测和调查描述了鄱阳湖地区优势植被种群的分布格局^[20]，以及外来入侵种的科目和对于鄱阳湖湿地的影响^[21].

本文使用的谐波分析方法不仅在信号处理领域应用广泛，随着遥感数据的积累，也常见于对时间序列遥感数据的处理.谐波分析法能准确描述时间序列数据的变化规律，在地表覆盖分析方面的应用已经比较成熟^{[25 — 34]}.于信芳等利用谐波分析的方法对我国东北森林物候期进行监测，并获取了其空间分布格局^[33].本文使用的遥感数据具有高时间分辨率的优点，可以从多角度综合反映地表变化.谐波分析的方法则可以将反映地表变化的时间序列信号简化为不同的特征参量.鄱阳湖因其湿地结构波动明显，相对于瞬时数据，谐波分析结果能够更客观地反映鄱阳湖地表综合特征.

1 研究区概况

鄱阳湖位于江西省北部，属于吞吐型湖泊.湖泊水面面积在夏季丰水期可达3700km²以上，在秋、冬枯水期可缩小至不足1000km²，形成独特的“丰水一片，枯水一线”的景观.研究表明，在吴淞高程12~16m的水位升高过程中，每升高1m，主湖水面以超过300km²的面积扩展.其中12~13m为面积增加最大区间，在12m升高到13m阶段，水面扩展662.3km².随着水位升高，草洲、滩地的分布减少明显.在12~14m区间，草洲随着水位每升高1m，平均以超过300km²的幅度减少^①.这种随水位而导致湿地结构发生波动的变化规律，是自古鄱阳湖形成以来的典型特征^[2].研究样区选取鄱阳湖国家自然保护区(以下简称鄱阳湖自然保护区)和南矶湿地国家级自然保护区(以下简称南矶自然保护区)(图 1).两个保护区为典型的水体-滩地-草洲交错分布的湖泊湿地结构.

① 鄱阳湖水利枢纽工程对湿地与候鸟的影响及对策研究课题组，李文华、刘兴士等.鄱阳湖水利枢纽工程对湿地与候鸟的影响及对策研究报告，2010年10月.

2 数据处理 2.1 数据来源及预处理

本文采用的数据主要有：MODIS-NDVI 16日集成的标准产品数据(空间分辨率为250m)、水文站定位观测数据、DEM数据.标准产品数据的时间段为从2000年第49d(2月18日)至2013年第353d(12月19日)，以16d为间隔，共计319期.MODIS-NDVI数据产品由MODIS逐日数据计算，经过16日最大值合成法(MVC)生产.MVC(Maximum Value Composite)指的是用合成期内最大NDVI值代表这16日的地表NDVI.从NASA提供的数据中提取出NDVI波段，通过遥感图像处理的标准化流程，并实现WGS84/UTM50坐标系的投影转换.

DEM数据(图 2)是将1 :25000比例尺等值线DEM中的等高线离散化成为点数据，并与1 :10000比例尺等高点数据空间叠加后生成不规则三角网(TIN)，在TIN的基础上内插值成规则格网DEM.^①

2.2 计算谐波分量特征值

谐波分析指的是由MODIS-NDVI数据产品得到每一个像元对应地表NDVI随时间变化的离散波信号，通过傅里叶变换将时间域的波形变化到频率域，用一系列正弦波来表示NDVI的波动^{[23, 28]}.不同的地表覆盖类型具有不同的NDVI变化模式，差别在正弦波的中心频率及振幅上.假设f(x)是[1，L]上的连续函数，则f(x)在闭区间[1，L]上的傅里叶展开为：$f\left( x \right) = \sum\nolimits_{n{\rm{ = 1}}}^\infty {\left( {{a_n}{\rm{cos}}\frac{{{\rm{2}}{\rm{\pi }}nx}}{L}{\rm{ + }}{{\rm{b}}_n}{\rm{sin}}\frac{{{\rm{2}}{\rm{\pi }}nx}}{L}{\rm{ + }}\frac{{\rm{1}}}{{\rm{2}}}{a_{\rm{0}}}} \right)} $, 式中，L为一个时间序列信号的时间长度.对于j≥1，f(x)的第j次傅里叶展开为${a_j}{\rm{cos}}\frac{{{\rm{2}}{\rm{\pi }}nx}}{L}{\rm{ + }}{b_j}{\rm{sin}}\frac{{{\rm{2}}{\rm{\pi }}nx}}{L}$，将第j次展开变换为单独余弦波来表示：

$ \begin{array}{l} {a_j}{\rm{cos}}\;\frac{{{\rm{2}}{\rm{\pi }}nx}}{L}{\rm{ + }}{b_j}{\rm{sin}}\frac{{{\rm{2}}{\rm{\pi }}nx}}{L} &{\rm{ = }}\sqrt {{a_j}^{\rm{2}}{\rm{ + }}{b_j}^{\rm{2}}} \left( {\frac{{{a_j}}}{{\sqrt {{a_j}^{\rm{2}}{\rm{ + }}{b_j}^{\rm{2}}} }}\cos \frac{{{\rm{2}}{\rm{\pi }}jx}}{L}{\rm{ + }}\frac{{{b_j}}}{{\sqrt {{a_j}^{\rm{2}}{\rm{ + }}{b_j}^{\rm{2}}} }}\sin \;\frac{{{\rm{2}}{\rm{\pi }}jx}}{L}} \right)\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ & = \sqrt {{a_j}^2 + {b_j}^2} \left( {{\rm{cos}}\;{\phi _j}\cos \;\frac{{{\rm{2}}{\rm{\pi }}jx}}{L}{\rm{ + }}\sin {\phi _{\rm{j}}}\cos \;\frac{{{\rm{2}}{\rm{\pi }}jx}}{L}} \right)\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ & = {c_j}\cos \left( {\frac{{2{\rm{\pi }}jx}}{L} - {\phi _j}} \right) \end{array} $

式中，${{c}_{j}}\sqrt{{{a}_{j}}^{2}+{{b}_{j}}^{2}},{{\phi }_{j}}=\arctan \frac{{{b}_{j}}}{{{a}_{j}}}+\rm{ }\!\!\pi\!\!\rm{ }$，令${{c}_{0}}=\frac{1}{2}{{a}_{0}}$，则有$f\left( x \right)=\sum\nolimits_{n=1}^{\infty }{{{c}_{n}}\left( \rm{cos}\frac{2\rm{ }\!\!\pi\!\!\rm{ }nx}{L}-{{\phi }_{j}} \right)+{{c}_{0}}}$，c_n表示第n次展开的振幅，φ_n表示第n次展开的相位.对于有限时间序列数据集，应用离散傅里叶变换(DFT)得到离散信号的分解展开式.展开后，关注每一个谐波分量的振幅和周期.振幅表示NDVI时间序列相对于平均值的波动，谐波分量的振幅越大，对应能量越高，说明原始信号中这一谐波分量代表的变化模式越明显.取振幅最大的谐波分量，计算其对应的变化周期(T).T反映对应地表NDVI在[1，L]内的变化趋势和极值的个数.

3 结果和分析 3.1 时间序列遥感数据谐波分析

分别对全部14 a与每一水文年(9月至次年9月)两个不同尺度的遥感数据进行谐波分析，得到地表NDVI主要变化谐波分量对应周期(T)，如图 3、4、5.从图 3可知14 a间鄱阳湖的湿地结构.根据谐波分析原理和不同地表类型在红光波段及近红外波段反射特点，将无明显周期性变化特征的像元判定为水体.T＝1a的区域为水体、含植被水体、滩地、水-滩过渡带等；T＝0.5a的区域是鄱阳湖的草洲带.鄱阳湖湿地主要的植被群落种类有芦苇群落、南荻群落、苔草群落等，这些群落存在不同程度的镶嵌情况，一般呈现出春季及秋季两次生长周期^[18]，即T＝0.5a.因此选择NDVI变化周期为半年的区域特征值表征湿地植被群落分布.湿地植被群落多分布于子湖周边以及入湖河流的各级支流周边.不同的地表覆盖具有不同的高程分布特征，水-滩过渡带主要在11~12m，草洲带主要分布在12m以上的高程范围.

图 4、图 5分别为两个保护区NDVI在每一水文年内的变化模式，例如T＝3~4月.在鄱阳湖自然保护区，湿地植被分布主要集中在蚌湖东侧带状区域、大湖池周边以及保护区东南众多小型子湖周边.当湿地植被分布减少时，大湖池周边的分布首先减少；当湿地植被分布增加时，其分布向主湖区蔓延.在南矶自然保护区，湿地植被在高程12m分布广泛，多集中于各子湖周边.

根据每一水文年和14a的谐波分析结果中T＝0.5a所占比例(表 1、图 6)，在两个保护区中，鄱阳湖自然保护区内湿地植被分布随着水位波动更大.以14a结果以±15 %为分析阈值，即鄱阳湖自然保护区结果阈值区间为[21.21 %，51.21 %]，南矶自然保护区结果阈值区间为[8.68 %，38.68 %].在阈值范围外的年份有2000-2001、2005-2006、2006-2007、2010-2011、2012-2013水文年，其余年份与14a结果相近.

根据鄱阳湖国家自然保护区管理局提供的2009年10月蚌湖实地调查结果(图 7、表 2)、张全军、胡振鹏、葛刚等学者关于鄱阳湖地区主要湿地植被种类物候期的描述^{[18 — 20]}(图 8)以及各点高程得出在2009-2010水文年各点NDVI的变化模式.A、B、C、D点主要植被覆盖类型为苔草科和芦苇，高程均大于13m.根据2009-2010水文年水位数据，A~D点所呈现的NDVI变化模式主要由湿地植被的物候期决定，因此谐波分析结果为T＝0.5a.而E、F、G点的高程小于13m，根据水位数据，在2009-2010水文年内，仅在其秋季生长期有长时间的露出，因此，对应的NDVI变化模式是以1a为周期的，与谐波分析得到的结果一致.

3.2 水位变化规律

鄱阳湖主要水文站有都昌站、棠阴站、星子站和吴城站等，其中吴城站水位取吴城赣江站和吴城修水站水位数据平均值.将都昌水位与其他水文站水位进行拟合，均呈现较好的拟合效果(图 9)，故选择都昌站水位数据进行分析.按3.1节中分析结果将这14a数据分成两组分别观察水位变化规律(图 10).

在谐波分析结果异常的水文年中9、10月和次年4、5月都昌水位相对于14a平均水位出现明显偏移.例如在分布比例超过60 %以上的2006-2007水文年，9-10月呈现低水位、小幅度减少趋势；在次年4-5月呈现低水位、小幅度增加趋势.而其他年份的水位偏移主要集中于11月至次年2月.

在2000-2001水文年的秋季生长期，即10月及11月，都昌水位比平均水位高出3m左右.而2005-2006水文年，水位波动与14a平均水位波动模式相似.但根据都昌每日水位数据，在该水文年丰转枯和枯转丰这两个变化过程中的高程10~13m区间，出现多个极值水位，造成水与滩之间湿地结构的反复变化，使得这些区域湿地植被的NDVI变化模式变成了以3~4个月为周期.在2010-2011水文年，10月至次年5月的平均水位均在9~10m.比较草洲带这一水文年与其他年份同一地表区域对应象元的Z谱曲线，在这一水文年的12月至次年2月，即主要湿地植被的枯萎期，NDVI没有变小，因此导致谐波分析结果出现绝大区域均为T＝1a的情况.类似的规律也曾出现在2012-2013水文年.

与之形成对比的是2006-2007水文年的分析结果.水位波动虽与2010-2011水文年相似，枯水位较低，但谐波分析结果与2010-2011水文年恰好相反，湿地植被分布高达60 %以上.2006-2007水文年秋季生长期水位较低，不仅为高程在13m以上的草洲带提供了湿地植被生长的必需条件，同时为更低高程区域的植被提供了足够的生长时间.该水文年的春季2-3月水位开始出现明显上涨，为湿地植被的春季生长提供了充足的水环境.因此在两个保护区内，一年两次生长期的湿地植被群落分布广泛.出现类似结果的还有2011-2012水文年.该水文年的9月及10月水位与14a平均水位相差3m以上，其湿地植被分布比例明显大于总14a的结果，尤其在鄱阳湖自然保护区.

3.3 湿地植被分布面积与水位的相关性分析

根据各水文站的地理位置，分别选择吴城站和棠阴站观测水位分析鄱阳湖自然保护区湿地植被分布面积与水位变化之间的关系(表 3).

相关分析结果表明，在鄱阳湖自然保护区和南矶自然保护区，湿地植被的分布与对应水文年9月的水位有紧密的负相关性，与10月水位相关性稍差，均呈显著相关，与其他月份的水位相关性较小.在2006-2007以及2011-2012水文年，9月及10月的水位比14a平均水位低2m以上，因此对应水文年的湿地植被分布比例高.在2000-2001、2005-2006水文年，9、10月吴城和棠阴水位均明显大于14a平均值，因此对应水文年的湿地植被分布比例低.鄱阳湖自然保护区湿地植被分布与水位的相关系数大于南矶自然保护区，说明鄱阳湖自然保护区的湿地植被分布面积对于水位变化更加敏感.结合鄱阳湖地区的水文特征，9月及10月的水位反映了对应水文年的退水时间，次年3、4月水位则反映了鄱阳湖地区次年涨水时间^[2].因此退水时间对湿地植被影响较大.

3.4 不同高程湿地植被分布面积对水位变化的响应

结合DEM数据(图 2)，在大于12m(吴淞高程)的区域，以1m为间隔统计湿地植被占保护区面积的比例结果(图 11)表明，在鄱阳湖自然保护区，湿地植被在12m以上的高程区间均有分布，且各高程分布相对南矶自然保护区比较均匀.而在南矶自然保护区，湿地植被主要集中于12~14m的高程区间，在大于14m高程上分布较少.

分别计算不同高程鄱阳湖自然保护区湿地植被面积与吴城水位的相关关系和显著性系数(表 4)，南矶自然保护区湿地植被面积与棠阴水位的相关关系和显著性(表 5).在不同高程，湿地植被的分布与水位的相关关系不同.

1) 12~13m：湿地植被分布面积主要受9、10月水位影响，且均显著相关；

2) 13~14m：鄱阳湖自然保护区的湿地植被分布面积同时与9、10月水位呈负相关，相关程度相近；南矶自然保护区的湿地植被主要受10月棠阴水位影响，与次年4月水位呈弱相关；

3)＞14m：湿地植被分布同时受到9、10月和次年4月水位影响.湿地植被的分布与9、10月都昌水位的相关关系在不同水平上(0.05~0.10)显著，但与次年4月水位关系不显著.相对于鄱阳湖自然保护区，南矶自然保护区的湿地植被分布与水位相关关系更显著.

4 结语

本文采用时间序列遥感数据的谐波分析，根据鄱阳湖地区主要湿地植被种类的物候特征，以NDVI变化周期为半年的区域特征值表征湿地植被群落分布.结合鄱阳湖的DEM数据与吴城站、棠阴水位数据，分析了湿地植被对水位变化的特征与规律，得出以下结论：

1) 两个保护区中，鄱阳湖国家级自然保护区湿地植被群落分布对水位的变化更为敏感.鄱阳湖自然保护区的湿地植被分布面积与9、10月吴城水位相关系数分别为－0.71和－0.66，均大于南矶自然保护区与棠阴水位的－0.66和－0.62.

2) 两个保护区范围内的湿地植被分布与对应水文年9、10月的水位呈现比较强的负相关性，且在0.05水平上显著.

3) 不同高程的湿地植被群落分布与水位关系不同.在吴淞高程12~13m的湿地植被群落分布与9月水位相关性最强，且显著相关；在13~14m高程的湿地植被分布与10月水位相关关系的显著性强于9月水位.南矶自然保护区在各高程的湿地植被分布与水位相关关系的显著性强于鄱阳湖自然保护区.

4) 鄱阳湖地区呈现一年两季生长的湿地植被受退水时间影响最大，其次是次年的涨水时间，与枯水期水位无明显的相关关系.

本文所用方法可以得到对于每一年地表变化的综合评价，避免了瞬时观测数据的偶然性误差，为不同年份之间进行动态对比提供了标准化的指标.采用谐波分析得到湿地植被的分布与水位关系结果，与叶春等使用生物量作为指标获取的结果^{[35 — 36]}变化趋势一致.

将T＝0.5a作为鄱阳湖主要湿地植被的表征，并不代表其他区域没有湿地植被分布，只是在250m空间分辨率的尺度上没有显现.

由于本文采用的MODIS-NDVI数据产品起始于2000年，采用谐波分析研究湿地植被与水位变化，需要在进一步研究中补充更长时间序列的高时间分辨率遥感数据.

鄱阳湖水利枢纽建设办公室提供的资料显示，鄱阳湖受人为扰动因素影响，湖区的61个碟形子湖中，有44个设有水闸进行水位控制，仅有17个子湖处于与鄱阳湖自然连通的状态.在人工调控背景下，对于水文节律改变的影响以及湿地植被对水文节律的响应关系等变化，需要在日后长期的研究中加以深化.

致谢: 鄱阳湖水利枢纽工程建设办公室提供水位观测数据，鄱阳湖国家级自然保护区管理局提供实地采样数据和吴城赣江、修水水位观测数据，在此一并表示感谢.