近年来, 城市黑臭水体治理一直是国内外专家学者研究的重点难题, 2015年8月《城市黑臭水体整治工作指南》提出“控源截污、内源治理;活水循环, 清水补给;水质净化, 生态修复”的黑臭水体整治基本技术路线。在控源截污的基础上, “活水循环, 清水补给”是提高城市内河水环境容量, 提升水体自净能力, 长效改善并维持内河水质的关键措施^[1,2]。但众所周知, 我国在面临严峻的水环境污染问题的同时也面临着水资源和可再生能源日益短缺的问题, 在“节水节能”的治理理念下, 如何充分利用地区水力优势, 实现水资源高效利用, 有效进行黑臭水体整治是亟需解决的关键问题。

　　 潮汐作用是潮汐河网地区得天独厚的自然资源^[3], 将潮汐资源化利用进行纳潮引水, 形成潮汐河网地区特有的“纳潮引水”活水循环方式, 是潮汐河网地区黑臭水体治理中最经济、最环保、最具特色的关键措施。

　　 高潮位变化规律分析与典型潮位过程选取是影响纳潮引水方案决策的关键, 而高潮位变化规律分析存在以下两个问题: (1) 每日潮位过程虽呈现出一日潮或半日潮特征, 但相邻高潮位值差异较大, 无显著周期规律可循, 仅以日典型潮位过程作为纳潮引水调度依据, 缺乏科学性; (2) 潮汐作用受月球运动影响, 而不同月份同一农历日高潮位之间亦存在显著差异, 无法单单以不同农历日来区别分析。因此, 在分析纳潮引水方案时, 基于潮位周期变化规律的典型潮位过程选取存在一定困难, 对方案的科学性与经济性产生严重影响。

　　 研究高潮位周期及其趋势变化特征, 分析潮汐河网地区纳潮引水调度周期, 进而选取科学的、切实反映潮汐动力条件的典型潮位过程, 不仅对潮汐河网地区“纳潮引水”活水循环技术的改进与完善具有深刻意义, 还可以为潮汐河网地区潮汐能源开发、利用与管理以及防洪排涝提供科学依据。

　　 目前, 潮位周期识别方法研究发展相对滞后, 而降水量、年径流量等水文序列特征分析的研究发展较为成熟。其中, Huang等^[4]于1998年提出的经验模态分解 (EMD) 方法, 是一种对非线性、非平稳信号进行平稳化处理的时频分析方法, 不需要建立在任何先验假设基础上, 可完整保留数据本征, 被广泛应用于降水量与年径流量特征的分析中。但EMD存在模态混叠的现象, 因此, Wu等^[5]在EMD分解中随机添加白噪声, 提出集合经验模态分解 (EEMD) 方法, EEMD分解可更为准确地反映时间序列的周期规律与趋势特性, 在水文序列分析中效果更为显著^[6,7]。因此, 本文将EEMD引入潮位周期分析方法研究中, 提出基于EEMD方法的潮汐河网地区纳潮引水调度周期分析方法。

　　 EMD的基本原理为基于序列的特征时间尺度, 经验地确定固有振动模态, 将时间序列分解为一系列具有不同时间尺度的固有模态函数 (IMF) , 每个IMF都是平稳信号, 必须满足: (1) 极值点数-零点数≤1; (2) 极大值与极小值包络线均值小于允许误差上限ε。

　　 EMD分解步骤如下^[8]:

　　 (1) 计算原始序列Y (t) 中的极大值点Y_max (t) 与极小值点Y_min (t) , 通过3次样条插值法对极大值与极小值序列进行插值, 形成上下包络线, 并计算其均值m₁ (t) =[Y_max (t) +Y_min (t) ]/2。

　　 (2) 计算类距平序列h₁ (t) =Y (t) -m₁ (t) 。

　　 (3) 将h₁ (t) 作为新的时间序列, 重复 (1) ~ (2) 步, 直到h₁ (t) 同时满足作为IMF的2个条件, 将满足条件的h₁ (t) 作为IMF, 表示为c₁ (t) 。

　　 (4) 从原始序列Y (t) 内去掉高频分量c₁ (t) 特征, 得到新序列r₁ (t) , 因此, r₁ (t) =Y (t) -c₁ (t) , 将r₁ (t) 重复 (1) ~ (3) 步, 得到第二个IMFc₂ (t) 。

　　 (5) 重复 (1) ~ (4) 步, 得到新IMF分量, 记为c₃ (t) 、c₄ (t) , 其中:

　　 设置门限值ε, 判断当SD≤ε时, 其中, 迭代停止。原始序列, 其中r_n (t) 为趋势项, 反映原始序列的变化趋势, IMF分量表示原始序列中不同时间尺度上的波动, 周期由小变大。

　　 模态混叠^[9]指的是通过EMD方法分解得到的IMF分量中, 同一个IMF分量包含不同时间尺度特征成分, 或者一个时间尺度特征成分存在于不同的IMF分量中的现象。模态混叠对EMD分解结果产生的影响包括: (1) IMF分量呈现出锯齿状曲线分布, 使IMF中包含的物理意义失效; (2) IMF的分解结果失去物理单一性, 相同的两个原始序列EMD分解结果不一致。

　　 基于EMD分解方法存在的问题, EEMD通过随机加入频率分布均匀的高斯白噪声^[10], 以均衡不同信号特征, 避免序列时间尺度的跳跃性变化。其零均值的特性在有效抑制模态混叠现象的同时, 完整保留了原始序列的时间尺度特征。

　　 EEMD分解步骤如下^[11]:

　　 (1) 在原始序列Y (t) 的基础上加入高斯白噪声序列n_i (t) , 转换为Y_i (t) , 其中Y_i (t) =Y (t) +n_i (t) 。

　　 (2) 对序列Y_i (t) 进行EMD分解, 得到相应的IMF分量。

　　 (3) 重复 (1) ~ (2) 步, 每次添加的高斯白噪声序列需保持不一样, 重复N次, 为消除人工添加白噪声带来的影响, 取得到的IMF分量平均值作为最终IMF分量c_i (t) , 其中

　　 快速傅里叶变换^[12]的基本原理为基于旋转因子W_N=e^-2j/N的周期规律与旋转对称特征, 对离散傅里叶变换 (DFT) 的计算过程逐步分解, 得到新DFT点, 对于一组原始序列x (n) , 其DFT变换为:

　　 基于旋转对称的特征以及周期规律, 将表达式整理为x (k) =A (k) +W_N^kB (k) 其中, , 再对A (k) 、B (k) 进行以上分解, 直到将N点最终分解为两点, 通过信号频谱分析识别序列的周期性规律。

　　 利用FFT变换提取IMF分量的平均周期信息, 步骤如下^[12]:

　　 潮位大小是潮汐作用与流域水文循环过程的共同作用结果, 河川径流量的大小是决定潮位值的高低的关键因素, “纳潮引水, 活水循环”的水源补给是影响纳潮引水方案的重要因素。因此, 本文针对流域径流量变化规律的分析, 分别对丰、平、枯水年纳潮引水的调度周期进行基于EEMD的全面分析, 以保障调度周期计算结果的科学性, 提高计算结果的应用价值。

　　 对时间序列进行EEMD分解时, 序列端点信息会失真, 随着分解次数的增加, 端点效应会扩散至序列内部。因此, 本文采用镜像延拓法^[13]对原始序列两端进行延拓, 以剔除因EEMD分解而产生的端点效应。

　　 EEMD提取的周期信息包含多个不同的时间尺度, 为使纳潮引水调度周期分析结果具有更广泛的有效性, 本文利用T检验^[14]与F检验^[15], 对EEMD方法提取的IMF信号中的周期信息进行有效检验, 同时通过T检验与F检验的平均周期为合理周期值, 以此确定纳潮引水调度周期, 提高周期分析结果的科学性。

　　 福州市南台岛位于闽江北港、乌龙江交汇处的冲积平原, 总面积120km², 内河水系属感潮河段, 四面环江, 受到潮水双向顶托作用。南台岛内河水体污染严重, 区域内黑臭水体占水系总数的58.5%, 缺乏水体自净能力。而与南台岛相连的闽江与乌龙江水质基本维持在地表水Ⅱ类水及以上, 南港水质优于北港水质。结合区域天然地理优势, 在控源截污的基础上, “纳潮引水, 活水循环”为南台岛黑臭水体治理关键措施。

　　 本文选取闽江解放大桥站1956~2015年逐日高潮位资料, 根据闽江年径流特征^[16], 以2009年、1975年、1966年为典型丰、平、枯水年, 对南台岛纳潮引水调度周期进行分析。

　　 闽江潮水为典型半日潮过程, 即一日出现2次高潮位, 本文将一整年中出现的高潮位作为一个时间序列, 利用EEMD方法进行周期识别。

　　 2009年共出现708次高潮位, 经镜像延拓法将原始序列长度延拓至1 024次, 利用EEMD方法对潮位序列逐步分解, 得到6个IMF分量与1个趋势项, 分别去除延拓部分, 得到原始序列内在IMF信息, 如图1所示。

　　 由图1可知:

　　 (1) 丰水年典型高潮位序列分解得到6个IMF分量, 表示丰水年高潮位序列包含多个时间尺度特征。

　　 (2) IMF2与IMF3保持了较高的相似性, 2月与7月高潮位波动较大, 3~5月波动相对较小, 8月之后波动较为稳定。

　　 (3) 随着分解阶数的增加, IMF分量频率逐渐减小, 周期增大, 波动趋于稳定。趋势项R呈先上升后下降趋势, 表示高潮位随时间的推移, 呈现出先上升后下降的趋势, 在7、8月达到最高, 这与福州呈现7、8月主汛期的特征相符。

　　 对IMF1~IMF6进行FFT变换, IMF1的FFT频谱如图2所示, 由图2分析可知, IMF1的平均周期为2.08 (周期电位为12h, 下同) 。

　　 同理, 分别对IMF2~IMF6进行周期提取, 并对平水年与枯水年高潮位序列通过EEMD方法进行周期识别, 结果如表1所示。

　　 采用T检验与F检验分别对丰、平、枯水年的EEMD方法周期识别结果进行检验, 显著性水平选为0.05, 考虑T检验对小样本检验更有效, 本文根据平均周期大小与样本长度对T检验的自由度进行适应性调整, T检验与F检验结果分别如表2与表3所示。

　　 由表2与表3分析可知, 丰、平、枯水年高潮位平均周期计算结果中同时通过T检验与F检验的周期结果分别为:28.16、28.24、28.24。

　　 根据以上分析可知, 解放大桥站高潮位平均周期在闽江丰、平、枯水年之间差异不大, 均呈现出28次高潮位一个周期的特征。因此, 福州市南台岛活水循环技术中纳潮引水调度周期应取为14d, 即根据潮水涨落规律, 以14d为一个周期, 选取相应的典型潮位过程, 作为纳潮引水调度模型潮位边界条件, 制定纳潮引水调度方案, 更为科学。

　　 (1) 本文创新性地将EEMD信号处理分析方法引入潮位序列分析方法中, 并考虑了流域径流量特征对潮位序列的影响、时间序列周期分析中的端点效应以及周期的合理性检验, 提出了基于EEMD分解的潮汐河网地区纳潮引水调度周期识别方法。EEMD分解可得到时间序列多个时间尺度特征, 综合体现了原始序列的局部特征与整体特征, 周期分析结果更为全面;

　　 (2) 以福州市南台岛为例, 利用本文提出的周期识别方法对解放大桥站历史实测潮位资料进行分析, 丰、平、枯水年高潮位序列均分解得到6个IMF分量与1个R趋势项, 利用FFT变换提取平均周期, 并进行T检验与F检验, 得到丰、平、枯水年平均周期为:28.16、28.24、28.24, 进而得到调度周期为14d, 分析结果合理可行, 进一步验证了周期识别方法的科学性与可行性。

　　 (3) 潮位序列既受水文循环规律影响, 又受月球运动规律影响, 变化规律错综复杂, 本文提出的纳潮引水调度周期识别方法仅为初步探究, 通过不同地区的周期分析检验、不同方法的对比检验以及纳潮引水实际调度工作检验将是下一步研究工作的重点。