应对水文序列非一致性变化影响的溯源重构法研究

秦 毅，李 时

（省部共建西北旱区生态水利国家重点实验室西安理工大学，陕西西安 710048）

1 研究背景

近年来，大量研究成果表明实测水文序列存在明显的非一致性特征［1-4］。它带来两方面问题：（1）在非一致性条件下，继续使用以一致性样本为基础的预测分析方法，包括传统水文频率分析方法、多元统计分析的Coupla方法、机器学习等方法是否能获得合理的成果，以及对未来的预测应该如何评价。（2）水文模型参数的率定工作变得困难，因为参数率定采用资料的非一致性导致了不确定性的增加。相比较而言，第一个问题在设计分析方面显得更为突出，例如水利工程、水土保持工程及生态工程的设计，洪旱灾害的防治等都离不开对实测样本序列的统计分析与预测，除非研究对象的物理机制被真实地掌握，否则所做出的水文分析与预测的结果只有在一致性（平稳性）序列条件下才有意义。因此学者们投入了大量研究，包括：（1）对水文序列非一致性变异诊断的深入研究，给出了相对成熟的诊断方法：如相关系数检验法、Spearman秩次相关法、Kendall 秩次相关法、Mann-Kendall检验法［5-6］、R/S分析方法［7］、非参数Pettitt检验法［8］等。在此基础上，谢平等［9］构建了水文变异诊断系统，更全面地反映了水文序列的变异特征。（2）解决非一致性条件下的频率分析问题，这是工程设计迫切需要解决的问题［10］。为应对非一致性序列的影响，除了采用水文分析计算相关规范提出的还原或还现方法外，近年来有学者们将研究焦点放在了时变矩法［11-12］上，并据此建立时变概率分布，再通过重现期或可靠度与时变概率分布间的关系确定出工程设计值。当前比较有代表性的方法包括基于修正重现期的EWT［13-14］和ENE方法［15-16］，和基于可靠度的ER［17］、DLL［18］和ADLL［19］等方法。

然而，这些方法仍然存在一定的争议。首先是对趋势的看法，有学者认为观测到的水文序列的趋势并非是非一致性特征的表现，而可能是平稳过程中的周期性频率波动所致［20］，甚至认为趋势这个词都缺少明确的定义［21］；再者，实践中，在工程设计寿命内，一定可靠度下的设计值将随着计算时段的起始时间以及统计参数与协变量间关系拟合曲线的类型而变化，这就意味着未来设计值的可靠性严重依赖统计参数的时变特征，而这种时变特征却因为时间序列遍历性的缺失，造成对未来水文序列统计参数预测的不确定性［22］。正如Serinaldi等［23］指出的，如果模型结构无法通过演绎推断的方式获得，特别是描述分布参数非平稳变化模型是仅通过归纳推理拟合得来，则一个额外不确定性来源被无形中引入了模型结构中，使得所得到的非平稳模型无法在实际应用中提高设计值的可信度与准确性，而模型的这种不确定性将很容易导致不具物理意义的结果产生；此外，利用时变参数识别非平稳序列的分布函数过程中，所采用的经验频率计算方法也是引起不确定性的重要因素，因为计算经验频率的方法也同样要求样本为简单随机样本，基于非平稳序列直接进行经验频率计算显然不恰当。可见，对于统计特征随时间变化的单一随机变量，其分布研究尚存在如此多的问题，更何况解决多元时变随机变量的联合分布问题，难度可想而知。

很显然，造成“难度”的核心是样本序列的非简单性。如果能够将非一致性样本序列转换为一致性序列，则因为对此简单序列已经有了成熟的分析理论、技术和用好技术的方法，前述困难将不复存在。因此，本文提出溯源重构法，它是基于反映水文现象物理机理的源函数将非一致性序列重构为新的一致性序列的方法。依靠重构的一致性序列，解决需要基于一致性序列才能进行分析工作的问题。事实上，当前已存在将非一致性序列转变为一致性序列的方法，即还原或还现法、分解合成法，本文将在后面就这两种方法与溯源重构法进行比较与讨论。

2 溯源重构法的理论与方法

2.1 源函数的概念一般来说，一个水文变量Y的变化是多个影响因子X1，…，XN共同作用的结果，见图1。在排除其他影响因子影响的条件下，单独考察影响因子Xi对水文变量Y的作用，将这种作用规律记为fi（Xi），则fi（Xi）是Xi对Y的物理作用机制，因此这里称fi（Xi）是Xi作用于Y的源函数。作为物理机制，源函数是不会改变的，就像河道流量Q=AV，其中A是过流断面的面积，V是流速，则V和A各自对Q的作用规律或说机制，也即源函数就是fV（V）=V，fA（A）=A。当V或A或两者共同随时间变化时，流量也会随着V（t）或A（t）或两者在不同时刻的取值而表现出随时间变化的Q（t）。但是，不论V和A怎样取值，源函数fV（V）=V，fA（A）=A保持不变，而流量等于这两个源函数相乘的机制也不会改变。即物理机制不随时间变化，变化的永远是Xi的状态。

图1 各影响因子对研究变量Y的源函数作用关系

严格的说，源函数的确定有两种方法：理论推导（数学解析）和实验分析。实际上，由于水文现象变化的复杂性，即便我们认可源函数的存在，也会因为人类对自然的有限认知而很难得到绝对准确的物理机制数学表达，至少当前很少看到对单一因子影响的研究。由于事物之间的作用规律可以隐含于统计规律中，所以当下可以借用Y和Xi间的统计关系来近似估计源函数。

既然源函数或由源函数构成的函数不随时间变化，我们就可以利用这个特征将非一致性序列转变为一致性序列。

2.2 溯源重构方法的基础模型欲利用源函数时不变的特点，就要表达出多个影响因子各自对Y作用后的一般综合关系。鉴于数学解析方法暂时行不通，故采用建立统计模型的方法来表征研究对象Y与解释变量源函数之间的关系。众所周知，数学的基本运算包括加（减）与乘（除），结合运算法则，Y与解释变量源函数之间的关系有两种一般性广义表达方式，一是各解释变量间相互独立时采用的叠加模式：

另一个是可以考虑解释变量间存在相关性的乘积模式：

其中：N指作用于研究变量Y的所有影响因子的个数，这些影响因子包括人类认识到的和尚未认识的；fi（Xi（t））为影响因子Xi的源函数。假设在N个影响因子中，有n个因子是人类目前已认识到的，它们的“源”函数分别为f1（X1（t）），…，fn（Xn（t））；未知影响因子及其“源”函数分别为fn+1（Xn+1（t）），…，fN（XN（t））。则式（1）和式（2）改写为：

分别记δ=fn+1(Xn+1(t))+…+fN(XN(t))和δ*=fn+1(Xn+1(t))…fN(XN(t))，他们代表人类未知的影响因素对Y的作用，通常可以认为是自然随机变量，具有概率分布，并具有对应的均值μ、μ*，和方差σ2、σ*2。令Ks(t)=f1(X1(t))+…+fn(Xn(t))和Kp(t)=f1(X1(t))…fn(Xn(t))。在生产实践当中，由于我们更关心影响因子取给定值时，Y的变化情况，与回归分析的考虑相同，在计算Y之前，解释变量X1(t)，…，Xn(t)的取值已经确定，可以作为可控变量处理［24］。于是Ks(t)和Kp(t)在t时刻为常数，故Ks(t)和Kp(t)是已知的时间函数。重新写上述式（3）、式（4）为：

对比两种表达式下Y（t）的条件数学期望和条件方差如下。

叠加模型：

乘积模型：

不难看到，两种模型下Y(t)的条件均值均是时变函数，一旦Xi发生趋势性变化，Y(t)的条件数学期望也就出现趋势性变化；与均值不同，在叠加模型下的方差为常数，而乘积模型下方差是时变函数，且会随Xi的趋势性变化而产生更强的趋势性。

现在来考察水文现象变化的普遍特点。水文现象是在水文系统内各要素间普遍存在相互影响的背景下发生的。各解释变量间基本不存在绝对的独立性，他们的相互作用使水文系统呈现出高度复杂的自适应非线性状态。随着近年来的环境变化，大量的研究成果表明不仅水文变量的量级表现出了非平稳变化，而且水文极端事件也急剧增多［25-26］。从统计的角度来看，这种现象表明了水文序列的一、二阶矩均发生了变化［27-28］，特别是极端事件的极端特性越强，二阶矩的变化就越不能被忽视［29］。因此，可以得到结论，乘积模式相较于叠加模式更适用于描述具有复杂相关性的水文系统。故这里采用乘积模型（11）作为溯源重构法的基础模型。

2.3 溯源重构函数与重构序列考虑到Y（t）的非一致性变化是因为解释变量的非一致性变化引起，那么如果将所有带有非一致性变化的影响因子的作用从Y（t）中剔除，见式（12），则剩余量将表现出一致性。

在式（12）中，m是非一致性变化的影响因子个数。其右端为Y（t）剔除全部非一致性因子影响后剩余一致性影响因子以及尚未被人类认知的影响因子源函数的乘积，它是平稳的随机过程。如果只剔除了部分非一致性变化因子X1，X2，…，Xl（l

理论上，当引起Y非一致性变化的影响因子全部被识别出来，即l=m，且这些因子的源函数fi（Xi（t））已有准确的表达，那么由溯源重构函数计算出来的溯源重构序列RSt={rs1，rs2，…，rsn}为一致性序列，可用于需要一致性序列才能进行工作的情况。生产实践中，因人类认知和现有方法的局限，我们不可能得到绝对的平稳过程，这也不是我们的目标，但随着非一致性变化的解释变量对应的正确源函数不断加入溯源重构函数中，则溯源重构序列将逐渐趋向并最终趋向于平稳，接近于一个随机噪声。

需要提醒的是，由于估计源函数存在不确定性，所以使用中必须配合溯源重构序列RSt的一致性检验。若检验结果为一致，则同时接受所确定的影响因子及其源函数的估计形式，否则，仍需进一步研究影响因子及其源函数。

由于研究变量Y和影响因子间的因果关系，研究变量Y的各种非一致性变化（线性、非线性、突变）与影响因子的相应非一致性变化一致，这使得溯源重构法去除非一致性时并不受非一致性变化类型的影响，即具有普遍性。相较于对非一致性变化的Y进行纯数学模拟的时变矩法，溯源重构方法更注重找到引起Y变化的非一致性变化影响因子Xi和它对Y的作用规律，并以此规律将非一致性变化从Y序列中排除，这正是溯源重构法能够将非一致性变化的Y（t）变为一致性序列RS（t）的关键所在。

3 实例研究

3.1 研究区域与数据利用陕西省榆林市佳芦河流域申家弯水文站年最大洪峰系列资料作为实例研究，考察溯源重构法构建一致性序列的效果。由于溯源重构法的重点是依据非一致性变化影响因子对水文变量的作用（即源函数）来确定溯源重构函数的，因此在研究中，首先应研究确定非一致性变化的影响因子。实践中，可以依据水文现象的相关理论知识并结合实地踏勘调研等手段，初步确定研究变量的主要影响因子，之后对这些主要影响因子进行趋势分析，找出造成水文变量非一致性变化的影响因子。

佳芦河位于黄河中游的河龙区间，毛乌素沙漠南缘，陕北黄土高原北段，属于黄土高原沟壑区。河流由西北流向东南，至佳县佳芦镇木厂湾村注入黄河。河流全长93 km，流域面积1134 km2，出口控制站为申家湾水文站。流域洪水由降雨引发。由于气候较为干旱，自然条件较差，该地区容易发生风沙灾害，水土流失现象普遍。为了控制水土流失，流域内修建了数十座淤地坝。淤地坝一般由两大件组成——坝和泄洪卧管。为了使入库泥沙得到沉积，实现多拦沙少拦水的功效，淤地坝运行方式设计成：发生洪水时，高含沙洪水首先被拦截在库内，后由卧管分数天将洪水排出库外，且要求卧管排水速度缓慢，排水流量一般为1～1.5 m3/s。结果导致出口断面洪峰较建坝前明显减小，库坝泄流对流域出口的洪峰流量的贡献可以忽略不计。这种效果等同于流域产生洪峰的产流面积减少，减少的数量为库坝控制面积。这里称流域面积与库坝控制面积的差为有效产流面积；另一方面退耕还林以及淤地坝的蓄水，使得坡面较大范围内的土壤含水率增加，促进了植被增长，导致流域植被截留、下渗、汇流等条件均发生改变，也造成了洪峰流量序列的非一致性变化。因此，选择有效产流面积与植被覆盖度作为佳芦河流域年最大洪峰流量序列非一致性变化的主要影响因子。

洪峰流量序列采用申家湾水文站1957—2010年实测数据进行分析。1957—2010年有效产流面积序列则根据流域面积与逐年修建淤地坝的控制面积计算确定。植被覆盖度（FVC）是利用landsat4-5 TM 30 m×30 m遥感资料，通过归一化植被指数（NDVI）的分析计算得到［30］。由于植被覆盖度序列长度受限于遥感资料的年限，仅有1988—2010年23年资料。但在1980年代前，流域人类活动影响较弱，且退耕还林等生态工程尚未起步，因此1980年代以前的植被状态较为稳定，本文将1988年前的植被覆盖度均取值为1988年数据。

图2和图3分别展示了佳芦河流域出口断面水文站申家湾站的年最大洪峰流量序列的变化和流域有效产流面积的变化。可以看到，流域洪峰流量和有效产流面积在1975年左右均发生了突变。图4显示了流域植被覆盖度的变化情况。对各序列采用Pearson 相关系数检验判断序列的线性趋势，Spearman 秩相关系数检验判断包括非线性趋势在内的一阶矩趋势，采用Breusch-Pagan 检验［27］判断二阶矩的趋势。因1988年前的植被覆盖度数据为作者外延处理确定的，对植被覆盖度趋势性检验时采用1988—2010年的数据进行分析。由于淤地坝因逐年修建而增多，且建坝本身不属于随机事件，所以有效产流面积序列呈现显著单调下降的趋势，也没有一般随机序列的波动变化。

图3 佳芦河流域有效产流面积序列

图4 佳芦河流域植被覆盖度序列

检验结果（表1）表明年最大洪峰流量序列的均值和方差都出现了显著减少趋势，但植被覆盖度的均值出现显著上升趋势。由于年最大洪峰流量一般在汛期发生，本文对佳芦河流域申家湾站的各时段致洪降雨序列进行了趋势性检验。考虑到黄土丘陵沟壑区降雨过程一般历时较短，不超过12 h，因此最大12 h降雨能够代表相应洪水的致洪暴雨，故在此处仅展示最大12 h降雨变化情况（图2）。同时将与佳芦河流域较近的榆林气象站最大日降雨序列作为参考，一并进行检验。检验结果（表2）显示，佳芦河流域致洪降雨序列未发生显著变化（包括线性和非线性一阶矩、二阶矩）。因此排除了降雨变化造成洪峰流量非一致性变化的可能。

图2 申家湾站年最大洪峰流量序列与致洪暴雨序列

表1 佳芦河流域洪峰流量序列及其影响因子的趋势性分析

表2 佳芦河流域致洪降雨序列趋势性分析

3.2 溯源重构序列的生成通过对流域实测资料的分析，我们可以得到这样的认识：佳芦河流域的降雨未发生显著的非一致性变化，申家湾站年最大洪峰流量序列的非一致性变化主要产生于流域有效产流面积和植被覆盖度的显著变化。因此，根据溯源重构法的思想，重构时不考虑降雨因子，仅考虑有效产流面积与植被覆盖度两个因子。在对洪峰流量序列进行重构之前，首先要分析出有效产流面积和植被覆盖度对洪峰流量作用的源函数。

3.2.1 源函数的建立方法 源函数是溯源重构法的关键，也是提升对水文过程认知的重要基础，是不应该被忽略的学术问题。它的确定应该来自于理论推导或实验，但现阶段，介于水文过程的复杂性，基于理论和实验手段成功确定出的源函数几乎看不到。因此，本着影响因子的作用规律必然体现在统计规律中的认识，我们可以通过回归的方法来逐个估计源函数的近似值。为叙述方便，这里令Y表示洪峰流量，X1和X2分别表示流域有效产流面积和植被覆盖度。源函数建立的具体做法：首先建立第一个因子X1与Y的关系f1（X1（t）），并作为X1的源函数作用于Y，之后将该因子的影响从Y序列中排除，即，再建立剩余序列与第二个因子X2的源函数f2（X2（t））。推广而言，若选择了m个有趋势变化的影响因子，则第3个影响因子的源函数由与第3个因子的关系估计，依次进行，直到m个因子的源函数被估计出。

显然，这种做法会带来因影响因子的引入顺序改变各因子源函数形式的问题。理论上，真实的源函数将不会受重构顺序的干扰。而实际中，由于源函数是估计函数，重构顺序的干扰是可能存在的，为避免这个问题，可根据各影响因子与研究变量之间相关性的强弱顺序确定重构顺序。

3.2.2 有效产流面积与植被覆盖度的单因子源函数 所谓单因子指只有一个非一致性性影响因子。众多的研究成果证明洪峰流量和流域面积之间呈现幂函数关系［31-32］。权琦泽［33］的研究结果表明榆林地区洪峰流量与流域面积的0.73次方成正比，故假定有效产流面积的源函数的估计为

植被对洪峰流量的作用关系研究不多见。根据周宏飞等［34］在径流试验场的实验研究成果，径流深与植被覆盖度FVC间呈现指数函数关系：

式中：R为径流深，mm；k为经验参数。借助洪峰与洪量间的幂函数关系、洪量与径流深间的关系，便可由式（15）得到洪峰流量与植被覆盖度间的作用规律：Qm∝eαFVC，因此植被覆盖度FVC对洪峰流量作用的源函数估计式为式（16）：

通过实测洪峰流量和eαFVC之间的回归，得佳芦河流域的参数α为-0.071，即：

3.2.3 年最大洪峰流量的溯源重构序列 当确立了有效产流面积和植被覆盖度这两个变化因子的源函数后，根据式（13）定义的溯源重构函数，建立佳芦河流域的溯源重构函数。为了比较溯源重构方法的有效性，分别建立单个因子和两因子的溯源重构函数如下：

利用有效产流面积的单因子重构：

利用植被覆盖度的单因子重构：

利用有效产流面积和植被覆盖度两因子重构（方法见3.2.1节）：

其中fFVC(FVC)=eαFVC中的α=-0.059则由RSA，t与eαFVC的回归得到。

理论上讲，若源函数是通过理论分析与实验分析的方式确定，则无论是单因子重构还是两因子重构，源函数不应发生变化。这里α 取值的差异正是源函数的估计受到重构顺序影响的体现。因此，实际应用中建议根据各影响因子与研究变量之间相关性的强弱顺序确定重构顺序。将洪峰流量序列Qm，t及其影响因子序列Ae，t和FVCt的值分别代入式（18）—（20）中，得到相应的溯源重构序列RSA，t、RSFVC，t和RSdouble，t。

4 讨论

按照本文提及的溯源重构理论建立的溯源重构方法，其实质是一种转换方法，所获得的重构序列RSt理应是一个具有一致性的随机序列。作为验证，本文以第3节中的溯源重构序列RSA，t、RSFVC，t和RSdouble，t为对象，通过考察其趋势性变化，探究溯源重构理论的正确性和溯源重构序列的合理性。同时，鉴于以下两个原因，对重构序列的特征进行探讨同样显得必要。①虽然RSt（如RSA，t、RSFVC，t和RSdouble，t）不是直接观测变量，但因它是由观测变量计算出来，因此有物理意义：相当于模数，应该具有它的特征；②实践中，在溯源重构时，随趋势性影响因子识别的完整性不同，所获得的重构序列会具有不同的特征。因此重构序列的特征分析对于溯源重构序列RSt的应用起着理论支撑作用。

4.1 溯源重构序列的趋势特征采用同样的假设检验方法，对上述溯源重构序列的趋势性加以检验，检验结果见表3。可以看出，具有显著趋势性变化的原洪峰流量序列经单因子溯源重构后，重构序列的各种相关系数均比原序列相应相关系数减小，说明不论是线性趋势还是非线性趋势均得到一定的消除；特别是双因子重构后，Pearson相关系数和Spearman秩相关系数均小于相应的临界值，说明溯源重构序列RSdouble，t已无显著趋势存在，且Breusch-Pagan 检验统计量c也由25.406 减小到5.157，尽管仍大于临界值，但洪峰流量序列二阶矩的非平稳性得到很大改善，说明引入非平稳变化的影响因子越合理，重构序列的平稳性越好。不仅如此，原序列存在的突变也得到消除，见图5。

图5 溯源重构序列与年最大洪峰流量序列的对比

表3 年最大洪峰流量溯源重构序列一致性检验

检验结果表明，溯源重构方法在消除非一致性序列的非一致性（包括线性、非线性趋势和突变）方面具有统计意义上的有效性。

4.2 溯源重构序列的分布特征分布特征是随机变量取值特性的完整描述，是研究随机现象的核心。这里讨论年最大洪峰流量溯源重构序列的分布函数，并与基于时变矩获得的洪峰流量序列（原序列）的分布函数进行比较。同时，针对当前工程师们无奈之选，即无视非平稳的存在，直接进行传统频率分析得出的分布函数也一并列入比较。

针对3.2 节得到的佳芦河年最大洪峰流量的溯源重构序列RSdouble，t，本文选取水文领域常用的Weibull（WEI），Log-normal（LN），Gamma（GA），Gumbel（GU）和GEV 五种分布类型作为备选分布，利用线性矩法［35］估计分布参数，依据纳什效率系数［36-37］、均方根误差RMSE等拟合优度准则确定最优分布。申家湾站溯源重构序列的最优分布为GEV，对应的分布参数分别为μ=9.657，σ=83.525，κ=0.878。对于非一致性序列的时变分布则基于GAMLSS模型［38］采用时变矩法确定，备选分布仍选取上述五种类型，以有效产流面积与植被覆盖度作为协变量，评价准则为AIC准则［39］，具体见表4。而无视非一致性存在，直接以非一致性原洪峰流量序列估计的总体分布也列于表4。

表4 不同方法确定的年最大洪峰流量序列的分布类型

理论上，年最大洪峰流量序列属于极值序列，它的分布应该是极值型。从表4可以看到不同属性序列和不同分析方法所挑选出的分布函数不同，其中只有溯源重构序列的分布符合极值型，本文作者对陕西省榆林地区大理河等6个流域的年最大洪峰流量溯源重构序列进行分布函数分析，结果均为极值类分布，这说明溯源重构法可以体现原序列的极值属性。这种既平稳又符合抽样特点的序列表明溯源重构法的理论正确，所得序列可用。对于非一致性原序列，相较于传统方法，时变矩法因考虑了原序列非一致性的特征，所以挑选出了与极值类型相接近的GA分布。有趣的是若直接应用非一致性原序列，采用传统方法分析出的分布函数却是不同于极值型分布的LN分布。这一结果从侧面反映了水文序列分布的确定会受到参数估计方法的影响，即估计方法越适应序列属性，所选择出的分布越可靠。

利用频率分析推求水文设计值是水资源利用、管理与保护领域的最基本问题。合理的设计值理应由合理的分布确定。Li等［22］在研究淤地坝影响下的设计洪水计算问题时，所采用的去非平稳性方法实质上正是溯源重构思想的体现。其获得的合理结果同样证明溯源重构方法获得的概率分布合理，同时说明该方法在非一致性序列频率分析方面具有良好的适应性。

4.3 与其他构建一致性序列方法比较溯源重构法的思想是通过溯源重构将非一致性序列转变为一致性序列，为以一致性序列作为基础的分析计算工作提供样本资料。这种想法同样体现在还原或还现法和分解合成法中。但由于方法所依据的原理不同，方法的应用效果随之不同。

还原或还现法依据的原理是通过对水文变量系统性变化的增量回加到原序列或对水文变量进行修正，达到使水文变量序列一致的目的。在增量的量化或修正值的量化工作中，由于要求资料的种类多、数量大，且对资料的准确性依赖程度高，使用不易。同时牵扯的一些概念难以精确界定［40］，存在方法失真或失效［41］等问题。即便用分布式水文模型直接计算增量，也会面临资料的非一致性导致模型参数率定（优化法）的不确定性加大问题［41］。正因如此，本文不做具体的分析对比。

分解合成法［42］认为，“无论水文现象的变化多么复杂，水文序列总可以分解成确定性成分和随机性成分这两部分。一般来说，水文序列的确定性成分主要受人类活动的影响，但并不排除气候因素（如气候转型期）和下垫面因素（如火山爆发、地震等）的影响，其变化规律可以在较短的工程年代里发生缓慢的渐变或剧烈的突变，因此水文序列中确定性成分的变化规律往往是非一致的。而水文序列的随机性成分则主要受气候、地质等因素的影响，其变化规律需要一个漫长的地质年代才能改变，因此水文序列中随机性成分的统计规律是相对一致的。”于是谢平等［42］建议采用成因分析法与统计分析法分别对确定性成分和随机性成分进行识别与检验，并对确定性成分进行拟合计算，对随机性成分采用传统方法进行频率分析计算。后者意味着分解合成法获得的随机变量样本是一个简单随机样本。下面将继续以佳芦河流域年最大洪峰流量序列资料，研究分析分解合成法生成的随机样本特征。

根据分解合成法原理，非一致性水文序列Xt可按下式表示：

式中：Yt为确定的非周期成分（包括趋势、跳跃等）；Pt为确定性的周期成分（包括简单的或复合周期的成分等）；St为随机成分。

为避免周期对趋势的影响，这里首先进行周期分析，利用小波分析［43］得到周期序列存在21年、9年两个显著周期，由此获得周期函数如式（22）：

剔除周期后发现序列仍存在明显的趋势，通过回归分析得到式（23）所示的趋势成分：

根据式（21）剔除周期成分与趋势成分后得到随机序列St，见图6，对该随机序列进行趋势检验，结果见表5。

图6 分解合成法构建的佳芦河年最大洪峰流量一致性序列

表5 分解合成法构建的佳芦河年最大洪峰流量序列的一致性检验

分解合成法获得的佳芦河年最大洪峰流量的重构序列十分有效地消除了一阶矩的显著线性趋势，但Spearman秩相关系数依然大于临界值，表明其一阶矩的非线性型非一致性仍然显著存在；方差的非一致性检验统计量χ为23.064，几乎没有任何改善，说明分解合成法重构的随机成分并非简单样本。

产生上述现象的原因可以概化为两个方面。一是由于定义两种成分的物理机制均不够清晰，例如分解合成法原理述及的“水文序列的确定性成分主要受人类活动的影响，但并不排除气候因素（如气候转型期）和下垫面因素（如火山爆发、地震等）的影响”与“随机性成分则主要受气候、地质等因素的影响，其变化规律需要一个漫长的地质年代才能改变”，其中气候转型期以及需要漫长的地质年代才能改变的气候因素无法清晰界定，更不能明确表达，因此确定性成分的拟合具有较大的不确定性，由实测资料推求出的随机性成分为实测资料与确定性成分之差，就会具有较大的不确定性；二是该方法假定水文序列Xt的各个成分满足线性叠加特性（即叠加模型）［44］。因此，根据上文叠加模型与乘积模型的对比分析，分解合成法构建的“相对一致”性序列仅能实现均值一致，而方差非一致的现象就不难理解了。这也成为了叠加模型无法描述序列的二阶矩变化的例证。

相比之下，溯源重构法认为水文变量Y的变化是一系列影响因子自身的非一致性变化引起的。这些因子的非一致性影响通过它们各自的源函数作用于Y。如果将全部非一致性变化因子的作用从Y序列中剔除，则剔除后的序列RSt就是一致性序列。与分解合成法不同，溯源重构法中的随机成分是指未被人类认知的次要因子的作用，无需专门率定。同时溯源重构法的基础模型采用乘法模式，更符合水文系统的非线性相依特征。溯源重构法的物理含义明确，分析的不确定性相对较小。

5 结论

众所周知，任何现象都有其因果关系，以数学术语来说，即因变量与自变量的关系。作者认为，水文变量Y的非一致性变化是由于自变量（影响因子）自身的非一致性变化所引起的。据此提出溯源重构法，即从源函数概念出发，依据因变量与自变量间的物理机理所建立起的统计模型，从因变量Y序列中剔除非一致性变化的自变量影响，实现非一致性序列Yt向新的一致性序列RSt的转化，从而可以服务于以一致性序列为基础的分析计算工作。溯源重构法在统计原理的基础上充分考虑了水文现象的物理成因，在原理上不同于分解合成法和还原或还现法。溯源重构法采用源函数描述影响因子Xi作用于Y的规律。源函数是一种机理函数，具有时不变的特性。同时溯源重构法的基础模型采用乘积模型，更符合水文系统的非线性相依特征。该方法物理含义明确，分析的不确定性相对较小。

依据佳芦河流域申家湾水文站的实测资料，应用溯源重构法实现了非一致性年最大洪峰流量序列向一致性的溯源重构序列的转化，结果表明溯源重构法能够有效消除年最大洪峰流量序列存在的突变及一阶矩的线性或非线性趋势，也能够消除二阶矩的非一致性。与分解合成法相比，该方法的目标也是试图获得一致性序列，但分解合成法仅能去除序列一阶矩的线性趋势，无法去除序列非线性趋势以及方差的非一致性。相比于时变矩法确定出的时变概率分布，以及直接应用非一致性原序列强行估计出的总体分布，溯源重构序列的分布函数更符合原序列极值抽样的统计属性。因此，溯源重构法所构成的序列类似于原序列的模数序列，除可以用于水文频率分析外，还可以通过对溯源重构序列构建预测模型，实现溯源重构法在水文预测预报领域的应用，并获得不确定性相对较小的预报成果。事实上，只要科学合理的表达影响因子X，例如本文第3节中，若佳芦河流域内淤地坝发生漫坝或溃坝情形时，只需要对X，即有效产流面积，进行合理修正，就可以继续沿用重构函数获得性质相同的序列，因此溯源重构法具有实际应用前景。