湖泊水位是衡量湖泊水文状况的一项最基本要素，其时空变化不仅受到气候、水文、地貌等自然因素的影响，与近年来日益增强的人类活动也存在紧密联系^[1-2].湖泊水量收支平衡直接决定着湖泊水位的丰枯，以降水和融雪作为主要补给源的湖泊，其水位波动具有明显的季节性和年际变化特征.自然状态下的水位波动是湖泊的固有属性，对维系湖泊生态系统健康及完整性至关重要^[3-4].此外，由于湖盆形态和区域水文条件的差异，湖泊水位的空间分布格局往往也呈现出异质性.

洞庭湖处于复杂的“江-湖-库”系统中，是与长江互联互通的重要调蓄湖泊，江湖水情变化对流域水资源利用、极端水情应对、水环境与水生态等产生了一系列影响^[5-6].学者们针对洞庭湖水情变化特征及成因分析^[7]、三峡工程运行对水情变化的影响^[8-9]、湿地生态系统对水位波动的响应^[10-11]等方面开展了广泛而深入的研究.詹泸成等^[12]利用稳定同位素技术研究了洞庭湖湖区河水、湖水和地下水稳定同位素的时空变化规律.赖锡军等^[13]运用江湖耦合水动力模型揭示了三峡汛末蓄水对洞庭湖各水域的影响具有明显的空间异质性，呈现出“北高南低，东强西弱”的格局.此外，洞庭湖水情演变还遵循着季节性规律，不同水文季节的时空演变特征也不完全相同.

目前，针对洞庭湖水位在不同水文季节以及不同水域的时空变化特征的研究相对较少，从水情关系上揭示影响湖泊水位时空演变特征的驱动因素仍有待进一步深入研究，研究手段也主要依靠遥感技术、同位素技术、水文模型等.因此，本文以洞庭湖水位为研究对象，主要通过统计学知识来研究洞庭湖水位在不同水文季节的时空演变特征及其驱动因素，从而为更加深刻地认识洞庭湖水位的时空分异规律提供科学依据.

1 材料与方法 1.1 数据源

水文数据由长江水利委员会和湖南省水文水资源勘测局提供，包括洞庭湖城陵矶、鹿角、营田、杨柳潭、南咀、小河咀6个主要水位站水位数据（吴淞冻结基面），长江干流宜昌站、螺山站流量以及“三口”（松滋、太平、藕池，对应水文站分别为新江口和沙道观、弥陀寺、康家岗和管家铺）总入湖流量、“四水”（湘江、资水、沅江、澧水，对应水文站分别为湘潭站、桃江站、桃源站、石门站）总入湖流量和城陵矶出湖流量，水位和流量均为1985-2014年逐日数据.宜昌站位于三峡水库下游，可以直接反映三峡出库流量；螺山站位于洞庭湖出湖口下游，可以直接反映长江荆江段流量和洞庭湖出湖流量，2个水文站的流量变化均可间接反映洞庭湖水位的变化情况.此外，为了克服以往研究中仅用某一水位站水位或多个水位站水位的算术平均值代表洞庭湖水位的局限性，本研究通过构建泰森多边形来计算洞庭湖水位.洞庭湖周边主要水系及各水文站分布如图 1所示.

1.2 研究方法 1.2.1 小波分析

小波分析是基于信号时频局部分析的一种方法，具有多时频分辨率功能，广泛应用在气象和水文要素长时间序列的周期性特征研究方面^[14-15].湖泊水位随时间变化而不断波动是多种因素共同作用的结果，属于连续非平稳序列，并包含了多时间尺度特征^[16]，通过小波分析能较好地揭示洞庭湖水位序列的周期性特征.小波函数是具有震荡性且能够迅速衰减到零的一类函数，常用的基小波函数有Haar、Daubechies、Mexican hat、Morlet、Meyer等.这里选取水文研究中最常用的Morlet连续复小波函数：

$\begin{split} φ(t)=\text{e}^{-\frac{t^2}{2}}·\text{e}^{i·w·t} \end{split}$

(1)

式中，ω为常数, i为虚数.

对于时间序列f(x)∈L²(R)，其连续小波变换为：

$\begin{split} W_f(a,b)=|a|^{-\frac{1}{2}}∫^{+∞}_{-∞}f(t)\left(\frac{t-b}{a}\right)\text{d}t \end{split}$

(2)

式中，W_f(a, b)为小波变换系数；a为尺度因子，反映频域特征；b为时间因子，反映时域特征；φ(t)为φ(t)的共轭复数.

小波方差反映了水文或气象要素波动能量的强弱，其随时间尺度的变化可以用来确定主周期.小波方差计算式为：

$\begin{split} \text{var}(a)=∫^{+∞}_{-∞}W_f(a,b)^2\text{d}b \end{split}$

(3)

1.2.2 聚类分析

聚类分析是将所有研究对象划分为有限的不同类别，从而使类内的相似度最大，而类间的相似度最小，常用方法有层次聚类、K-means聚类、模糊C均值（FCM）聚类等^[17].由于层次聚类适用于样本量少的分析且能充分识别研究对象的属性特征，故本文选择该方法对洞庭湖6个主要水位站进行聚类研究.

1.2.3 地统计

地统计是研究空间结构模型及其相关性的重要理论方法，目前已广泛应用于水文学、地质学、环境科学等领域^[18-19].区域化变量的空间自相关性可用变差函数γ(h)来描述，其表达式为：

$\begin{split} γ(h)=\frac{1}{2N(h)}∑\limits^{N(h)}_{i=1}[Z(x_i)-Z(x_i+h)]^2 \end{split}$

(4)

式中，Z(x_i)和Z(x_i+h)分别表示在空间位置x_i和x_i+h的变量值, N(h)为空间场中距离为h的样本点对数.

洞庭湖水位是随空间位置变化的一种区域化变量，具有连续性、随机性和结构性，利用变差函数能较好地揭示其空间特征.常用的变差函数模型有球状模型、高斯模型和指数模型，经过比较验证，洞庭湖主要水位站所形成的空间结构利用高斯模型进行拟合的效果最优.高斯模型通式为：

$\begin{split} γ(h)=\left\{\begin{array}{l} 0,&h=0\\ C_0+C(1-\text{e}^{-\frac{h^2}{a^2}}),&\ \ h>0 \end{array}\right. \end{split}$

(5)

式中，C₀为块金常数，C为拱高，C₀+C为基台值，a为相关距离，高斯模型的变程A=$\sqrt 3$a.其中，块金常数、基台值、变程是解释区域化变量空间相关性的3个关键参数.

2 结果与讨论 2.1 洞庭湖水位的时间变化特征 2.1.1 年内变化特征

洞庭湖水位的季节性变化特征十分明显，年内呈单峰变化且波动幅度巨大，表现出了典型的洪-枯交替变化规律（图 2）.在一个水文年，洞庭湖水位变化过程大致可分为4个阶段：12月至翌年3月为枯水期；4-5月为涨水期；6-9月为丰水期；10-11月为退水期.湖泊水位在枯水期降到年内最低，最低水位出现在1月(24.35 m)，随后湖泊水位开始缓慢上升，枯水期平均水位为24.76 m.随着汛期的到来，湖泊水位在涨水期持续上涨，上升速率为1.34 m/月.进入丰水期后，长江干流和“四水”入湖径流量明显增大，湖泊水位始终保持在较高水平，在7月达到峰值(31.16 m)，7月末开始小幅回落，丰水期平均水位为29.97 m，比枯水期高出5.21 m，可见洞庭湖水位年内落差悬殊.在退水期，由于长江干流水位下降对洞庭湖的拉空作用比较突出，导致洞庭湖水位持续下降，下降速率为1.44 m/月.显而易见，洞庭湖水位在退水期的下降速率快于在涨水期的上升速率，说明湖泊水位在不同水文季节的变化机理是存在差异性的.此外，不同水位站水位涨落差异也十分显著，落差较大的有城陵矶、鹿角和营田，其月平均水位变化超过8 m，尤其是城陵矶水位落差甚至高达9 m；相比较而言，杨柳潭、小河咀和南咀的年内水位变化幅度较小，其月平均水位落差约为4 m.这些水位涨落的差异除与洞庭湖西高东低的地势有关外，还受到三峡水库增减泄流量以及与长江干流水力联系紧密程度的影响^[20].

2.1.2 年际变化特征

用Morlet连续复小波函数分析洞庭湖年水位序列周期性变化特征，小波变换系数实部时频分布如图 3所示.图中实线表示小波变换系数实部为正值，对应水位升高；虚线表示小波变换系数实部为负值，对应水位降低.洞庭湖年水位序列在21~22、27~28 a尺度上的波动比较明显，均出现了3个高低水位交替震荡中心.其中，27~28 a尺度上的水位波动具有全时域分布的特点，且21~22 a尺度的低值中心嵌套在27~28 a尺度的低值中心内.

图 4a为洞庭湖年水位序列的小波方差随时间尺度的变化过程，有2个波峰，分别对应着22和28 a尺度的主周期.最大峰值出现在28 a尺度处，表明该尺度下的震荡周期最强，即28 a为第一主周期，22 a为第二主周期.根据小波方差分析结果，分别绘制28和22 a尺度对应的小波系数实部变化过程（图 4b）.在28 a尺度下，小波系数实部在1985-1987、1996-2004和2014年为负相位，表示水位降低，在1988-1995、2005-2013年为正相位，表示水位升高，湖泊水位经历了“偏低-偏高-偏低-偏高-偏低”的变化过程，平均周期约为18 a.在22 a尺度下，小波系数实部的振幅变化较大，也即湖泊水位变化幅度由高到低再到高，波动的平均周期约为14 a.

2.2 洞庭湖水位的空间变化特征 2.2.1 水位的空间分布格局

考虑到洞庭湖复杂的水系条件和特殊的湖盆形态，不同水域水位的空间变化特征必然存在差异性.位于南洞庭湖南部的杨柳潭和位于西洞庭湖的小河咀、南咀的水位变化趋势基本一致，峰型较为矮胖；而位于东洞庭湖的城陵矶、鹿角和位于南洞庭湖东部的营田的水位变化趋势基本一致，峰型较为高瘦(图 2). 也就是说，洞庭湖水位在空间上可能表现为西洞庭湖和南洞庭湖中西部水域（杨柳潭、南咀、小河咀）以及东洞庭湖和南洞庭湖东部水域（城陵矶、鹿角、营田）两类不同的水位演变特征.

为了验证上面的推论，通过层次聚类分析对洞庭湖6个主要水位站的水位进行空间聚类分析，结果如图 5所示.层次聚类分析的结果显而易见，即洞庭湖主要水位站分为两大类别，与图 2中得到的推论是一致的.其中，杨柳潭、南咀、小河咀水位相似度较高，在空间上聚为一类（Group 1）；而城陵矶、鹿角、营田水位相似度较高，在空间上聚为一类（Group 2）.在Group 1中，小河咀与南咀的水位相似度极高，可能是因为这2个水位站都位于西洞庭湖，其间通过洪道直接相连，水力联系更加紧密.在Group 2中，营田和鹿角的水位相似度更高.由于城陵矶位于洞庭湖出口，其水位受长江干流的影响更大，而位于东洞庭湖南部的鹿角和南洞庭湖东北部的营田距离长江相对较远，他们之间通过宽阔的洪道联通且受长江的影响比城陵矶小，从而导致了鹿角和营田的水位相似度更高.

2.2.2 水位变化的空间相关性

根据地统计学理论对洞庭湖6个代表水位站的水位进行分析，从而揭示其在不同水文季节的空间变异规律.各水位站水位在不同水文季节差异显著，且同一水文季节的水位差异也明显不同，枯水期水位差异最大，最低水位（城陵矶）比最高水位（小河咀）低7.34 m，丰水期水位差异最小，最低水位（城陵矶）比最高水位（南咀）低2.61 m，这与标准差检验的结果是一致的（表 1）.在洞庭湖水位由丰转枯的过程中，当城陵矶水位<26 m时，其与南咀水位变化的关联性变差；当城陵矶水位<27 m时，其与杨柳潭水位变化的关联性变差；当城陵矶水位<24 m时，其与鹿角水位变化才表现出差异性^[21]，尤其是在枯水期时各水位站的关联性为最弱.枯水期和涨水期的偏度系数为负值，而丰水期和退水期的偏度系数为正值，说明水位分布曲线随着时间的推移由向左倾斜转变为向右倾斜，且这种变化趋势由弱到强再到弱.

由于样本数据不满足正态分布，首先对其进行对数转换，然后通过高斯模型拟合变差函数，各项参数如表 2所示，其中残差平方和（RSS）的数量级均小于10^-6且R²>0.90，表明拟合效果较好.各水位站水位在不同水文季节均无明显的各向异性，恰恰表现出很强的各向同性（长轴变程与短轴变程之比近似为1），这主要是因为洞庭湖内部的水流方向无论何时都不变，即由地势较高的西洞庭湖流经南洞庭湖并最终汇入地势较低的东洞庭湖.块金常数C₀常用来衡量空间随机场的粗糙程度，而在实际情况中仅仅依靠块金常数无法比较不同空间随机场的粗糙程度，因此还需要结合块金常数与拱高之比（C₀/C）来说明空间随机场的相对粗糙程度，若C₀（C₀/C）越大，则绝对粗糙程度（相对粗糙程度）越大，反之则越小^[22].各水位站在不同水文季节的块金常数基本一致，但C₀/C却完全不同，具体表现为丰水期的相对粗糙程度最大，其次为涨水期和退水期，枯水期的相对粗糙程度最小.基台值（C₀+C）描述了区域化变量的变异强度，各水位站在不同水文季节的基台值大小依次为枯水期>退水期>涨水期>丰水期，说明各水位站水位在枯水期差异相对较大，而在丰水期则比较相近.变程A反映了区域化变量的空间延展性及空间相关性随距离变化的强弱，与基台值相反，各水位站的变程大小依次为丰水期>退水期>涨水期>枯水期，说明各水位站水位在丰水期的空间延展性和相关性最强，而在枯水期最弱.总而言之，洞庭湖6个水位站水位在丰水期所构成的空间随机场的可认知度最高，其次为退水期和涨水期，而在枯水期的可认知度最低.

洞庭湖内部距离长江越近且水力联系越紧密的水域受长江干流的影响就越大^[13].城陵矶作为洞庭湖的唯一出湖口，与长江干流直接发生水力联系，而东洞庭湖与南洞庭湖东部之间通过洪道连接，水力连通性相对较好，使得长江干流水位涨落的影响一直延续到南洞庭湖东部区域，所以城陵矶、鹿角和营田在空间上具有一定的相关性.随着距离的增加，长江的影响也会相应削弱，因而位于东洞庭湖南部的鹿角和南洞庭湖东部的营田具有更高的相似度.长江“三口”分流经由南咀进入西洞庭湖，加上长江洪水携带的大量泥沙导致洞庭湖淤积严重^[23]，洲滩快速发育，甚至在整个西洞庭湖和南洞庭湖中西部大片区域也被洲滩所演替，导致其水力联系较差，使得小河咀、南咀和杨柳潭在空间上聚为一类.由于洞庭湖水位的丰枯变化受到长江和“四水”季节性来水的影响，其在各水文季节的演变也处于动态变化之中.在枯水期，湖泊水面严重萎缩，洲滩完全出露，各湖区之间依靠狭小的河道相连通，水力联系发生阻断，因而导致各水位站水位在空间上不具有高度相关性.随着长江流域汛期的到来，湖泊水位逐渐上升，水面不断扩大.进入丰水期，湖泊水面持续扩大，7月达到最大值，各子湖泊也相互联系在一起组成一个大湖泊，湖泊水面比降减小，各湖区之间的水力联系加强，所以在丰水期各水位站水位所形成的空间相关性最大.在退水期，湖泊水位由丰水期末期的缓慢下降转而迅速下降，水面大幅缩减，各水位站之间的水力联系也随之减弱，直至枯水期达到最弱.所以，各水位站水位在涨水期和退水期的水力联系弱于丰水期，但强于枯水期，在涨水期和退水期所形成的空间随机场的可认知度介于丰水期和枯水期之间.

2.3 洞庭湖水位时空演变特征的驱动因素分析 2.3.1 逐步回归模型的构建及检验

湖泊水位演变与径流量的来源和变化息息相关，而洞庭湖水情变化与江湖水量交换存在密切联系，两者在不同水文季节的响应关系具有差异性^[24].因此，针对洞庭湖水位时空分异规律，通过逐步回归分析方法，分别建立两类水位（W）在不同水文季节与长江干流宜昌站流量（YC）、螺山站流量（LS）、“三口”总入湖流量（SK）、“四水”总入湖流量（SS）、城陵矶出湖流量（CLJ）之间的多元回归模型来探究洞庭湖水位时空演变特征的主要驱动因素.该模型的数学表达式为：

$\begin{split} W=\text{a}·YC+\text{b}·SK+\text{c}·SS+\text{d}·CLJ+\text{e}·LS \end{split}$

(6)

式中，a、b、c、d、e为回归模型的拟合系数.

由于自变量和因变量的量纲不同，数量级差别较大，所以在建模前对数据进行Z标准化处理，以消除数值大小和单位不同的影响.回归模型通过R语言构建，其优势是根据AIC（Akaike Information Criterion）准则^[25]选取最优模型，并同时给出了t检验和F检验的显著性结果. AIC准则是以保证数据拟合优度并尽可能避免过度拟合为前提，其值越小表明模型拟合优度越好.通过逐步回归方法，从模型包含所有变量开始，每次去除一个变量，直到比不去除该变量时的AIC值更高为止，从而在保证自变量数目的同时达到最优拟合度.基于以上原则对洞庭湖不同水文季节的两类水位分别进行逐步回归分析，拟合结果如表 3所示.

各回归方程系数的显著性由t检验给出，所有方程的系数都通过了0.05显著性水平，且部分系数通过了0.01或0.001显著性水平(表 3).回归方程的显著性由F检验给出，两类水位在各个水文季节拟合的回归方程均达到了0.001显著性水平.从决定系数R²来看，除Group 1水位在枯水期的R²为0.70外，其余各回归方程的R²>0.90.从RSE（Residual Standard Error）来看，除Group 1水位在枯水期的RSE为55%外，其余各回归方程的RSE均不超过31%.由此可见，自变量的选取具有代表性，回归方程的拟合精度较好.回归模型检验数据为2014年的Group 1水位、Group 2水位、宜昌站流量、螺山站流量、“三口”总入湖流量、“四水”总入湖流量和城陵矶出湖流量.根据表 3中的逐步回归模型，分别计算出两类水位在不同水文季节的回归值，并与实测值进行比较，结果如图 6所示. Group 1和Group 2水位的模拟值与实测值基本吻合，其相关性均超过98%，表明回归模型的构建是合理的，拟合精度较高.

2.3.2 结果分析与讨论

洞庭湖两类水位在相同水文季节以及每类水位在不同水文季节的驱动因素均存在显著差异(表 3).在枯水期，Group 1水位变化与宜昌流量、“四水”总入湖流量和螺山流量有关，Group 2水位变化与宜昌流量、“三口”总入湖流量、城陵矶出湖流量和螺山流量相关，其中螺山流量是两类水位的最主要影响因素.在涨水期，影响Group 1和Group 2水位变化的驱动因素完全相同，但同一因素的影响程度不同. Group1水位变化主要受城陵矶出湖流量、螺山流量和“三口”总入湖流量的影响，影响比例分别为39%、27%和22%，而Group 2水位变化主要受螺山流量的影响，其影响比例占到73%.在丰水期，Group 1水位的影响因素与涨水期相同，不同的是受螺山流量的影响大于城陵矶出湖流量的影响，Group 2水位与各因素均有相关性，并主要受到螺山流量的影响.在退水期，除螺山流量外，各影响因素对Group 1水位变化的贡献大小依次为“三口”总入湖流量>城陵矶出湖流量>“四水”总入湖流量>宜昌流量；Group 2水位变化仅与螺山流量和“三口”总入湖流量有关，其中螺山流量的影响比例为89%. Group 1水位在4个时期均受到“四水”总入湖流量的影响，在不同水文季节的最主要影响因子也有所不同. Group 2水位在任何水文季节都受到“三口”总入湖流量和螺山流量的影响，且“三口”总入湖流量表现出负影响，螺山流量则是最重要的影响因子.

一方面三峡工程的运行调度影响了长江中下游的水文情势，另一方面长江中下游水文情势的变化又影响着洞庭湖的水文状况.三峡工程的运行改变了长江中下游的输沙量和径流量，从而影响了长江“三口”和洞庭湖的冲淤特征，引起了湖泊形态和江湖关系发生改变^[26]，导致正常水文年的城陵矶水位在5-6月上升了1.36 m，在12月至次年4月下降了1.89 m^[27].螺山径流包含了长江上游来水和城陵矶出湖流量，其中长江径流在枯水期和退水期约占65%.研究表明，长江干流水文演变特征是城陵矶水位变化的主导因素，其贡献率达到65%^[28].由此可见，长江干流水文（螺山流量）变化特征对洞庭湖水位的影响十分显著，尤其是对Group 2水位的影响更加突出.对西洞庭湖和南洞庭湖西部而言，其水位的变化直接与“三口”和“四水”入湖径流的多寡相关^[29].当湖泊高水位退去后，该水域被分割成小湖泊群和泄流洪道，内部水力联系被大大削弱.这导致了位于长江“三口”入湖处的南咀水位变化主要受“三口”分流的影响，而位于资水入湖口附近的杨柳潭则受到洞庭湖流域来水的影响较大，从而使得Group 2水位演变特征在不同水文季节的主要影响因素也不完全相同.

综上所述，洞庭湖水位演变具有明显的时空分异规律，具体表现为：（1）西洞庭湖和南洞庭湖中西部水域（Group 1）的水位演变在不同水文季节的主要影响因素各不相同，在枯水期和丰水期主要受螺山流量的影响，在涨水期主要受城陵矶出湖流量的影响，在退水期主要受“四水”总入湖流量和城陵矶出湖流量的影响，这既与该水域同时接纳长江“三口”和湖南“四水”入湖径流有关，还与季节性的洲滩出露和淹没改变了该区域的水力联系有很大联系；（2）东洞庭湖和南洞庭湖东部水域（Group 2）的水位变化在各个水文季节主要受螺山流量的影响，其平均影响比例为75%，即由长江干流水情变化引起的对洞庭湖的拉空和顶托作用是其水位演变的主要影响因素.

3 结论

1）从时间尺度来看，洞庭湖水位年内呈单峰变化，具有典型的季节性变化特征，各水位站的水位涨落差异也十分显著；湖泊水位年际变化具有多时间尺度特征，其中28 a尺度为第一主周期，22 a尺度为第二主周期，分别对应着18和14 a的波动周期.

2）从空间尺度来看，洞庭湖6个主要水位站呈两大类分布格局，其中小河咀、南咀和杨柳潭水位变化特征一致，在空间上呈现聚类，而城陵矶、鹿角和营田的水位变化特征一致，在空间上也聚为一类；各水位站的空间自相关性随着水文季节的变化而变化，强弱依次为丰水期>退水期>涨水期>枯水期，并主要受江湖水力关系变化的影响.

3）近30年来，洞庭湖水位演变呈现出明显的时空分异规律，东洞庭湖和南洞庭湖东部水域（城陵矶、鹿角、营田）的水位演变在各个水文季节都主要受长江荆江段下游径流（螺山）的影响；西洞庭湖和南洞庭湖西部水域（杨柳潭、南咀、小河咀）的水位演变在枯水期和丰水期主要受长江干流水情变化的影响，在涨水期主要受出湖径流量的影响，在退水期主要受“四水”总入湖流量和出湖径流量的影响，这与该水域接纳长江“三口”和湖南“四水”来水以及广袤的洲滩发育有很大联系.