基于约束机制粒子群算法的水库库容曲线修正方法与流程

本发明涉及水利工程中的已建水库库容曲线修正方法，尤其是涉及基于约束机制粒子群算法的水库库容曲线修正方法。
背景技术：
：水库库容曲线是水库工程规划设计和运行管理的必备参数和基础成果之一，其精度对水库调控方式及效益产生重大影响。目前，我国大约有9.8万余座已建成水库工程，仅有少数大型的、特别重要的、泥沙淤积问题严重的水库会经常性的对库区地形进行测量，但由于测量代价较大，大多数水库库容曲线极少进行复核测量，若在运行管理期间库区地形形态发生显著性变化，对水库运行调度决策将产生重大影响。因此，如何采用成本较低的方式复测水库库容曲线具有广阔的现实需求。通常获取水库库容曲线有两种途径：一是直接测量法：对库区陆地和水下地形进行全面、高精度的测量；二是水量平衡法：基于水库运行资料，采用水量平衡方程对库容要素进行推算。直接测量法又分为地形图法和横断面法。地形图法通过全库区周边及水下范围的测量，“自下而上”进行累计计算，该方法所得成果可靠、精度高，但外业工作量巨大，对测量设备要求高，需投入大量人力、物力和财力，作业周期较长，因此该方法适用于库区范围较小、库区冲淤变化不显著的水库或库容影响重大的水库。水量平衡法根据水库不同水位高程下的入库、出库流量、降雨、蒸发、渗漏等观测资料，依据水量平衡方程计算对应库容差，综合形成库容曲线。该方法几乎不需开展外业工作，对库区范围、形态也没有限制，仅对运行资料的精度、时段、完整性等有一定要求，因而适用范围相对较广，成本相对较低。但是，水量平衡法也存在一定的局限性和问题：首先，由于水文测验精度有限，库周汇流难以精确测算，入库流量计算可能存在一定的误差；其次，水库蓄泄运用时，库岸坍塌、库区泥沙冲淤将引起不同高程下的水面面积发生变化，影响了库面降雨直接入库水量及库面蒸发水量的测算，进而增大入库、出库水量的计算误差；再次，由于库容值是由库容差值由低水位向高水位累加求和，导致库容误差被不断累计和传播；最后，水量平衡法仅得到库容值，尽管通过试算法可反推对应的水面面积，但会出现高水位水面面积小于低水位水面面积这种违背常理的不协调现象。因此，如何克服上述水量平衡法存在的问题，提高水库库容曲线精度是一个亟需解决的问题。技术实现要素：本发明目的在于提供一种基于约束机制粒子群算法的水库库容曲线修正方法，以不同水位对应的水面面积作为决策变量，通过构建符合变量内部相邻成员间变化规律的约束机制，运用改进的粒子群算法，求解与运行资料偏差最小(即离差平方和最小)的库容曲线。为实现上述目的，本发明采取下述技术方案：本发明所述基于约束机制粒子群算法的水库库容曲线修正方法，包括下述步骤：步骤1，构建决策变量内部约束机制：通过分析水库库区地形的一般变化规律，根据水库水位～面积函数的微分特性及实际应用时水库水位～面积～库容关系为离散形式，以不同水位对应的水面面积作为决策变量，构建符合变量内部相邻成员间变化规律的约束机制；步骤2，水库运行资料预处理：为便于分析和比较求解成果与运行资料的匹配程度，根据目标水库的运行资料、气象资料、库区地质资料以及近期库容曲线，由水量平衡原理，逐时段计算初、末水位及其对应的时段库容差k＝1，2，……，p，并以此系列作为基准库容差；步骤3，目标函数设计：针对某一备选库容曲线，在曲线上逐时段查得运行资料中的初、末水位及其对应的时段库容差值dvk，以此系列作为备选方案库容差；以时段备选方案库容差与基准库容差的离差平方和最小为所述目标函数，即：步骤4，基于约束机制的改进粒子群算法求解目标函数，步骤如下：子步骤4.1，设定算法基本参数，初始粒子数量m、粒子维度n、更新迭代周期上限s，学习因子c1、c2均取2，速度限制常数vmax取值为1/(2n)，惯性权重w采用线性递减自适应调整策略；子步骤4.2，建立粒子与面积函数决策变量的映射关系；子步骤4.3，基于步骤1中相关的约束机制，随机生成满足要求的初始粒子群体；子步骤4.4，对第s迭代周期，逐个粒子i计算相应的适应度值得到本时刻各粒子的自身最好位置及群体最好位置子步骤4.5，根据改进粒子群算法，粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置；子步骤4.6，重复子步骤4.4、4.5，直至粒子群更新迭代周期上限s时，输出群体最好位置，经换算后即为求解结果；步骤5，输出条件判别：设定允许距离阈值为δ*，计算步骤4中求解的库容曲线结果与步骤2中水库运行资料中库容曲线的欧氏距离δ，若δ≤δ*，则输出步骤4计算的库容曲线结果；反之，用步骤4的库容曲线替代步骤2中近期库容曲线，并重复步骤2、3、4，直至满足输出条件为止。所述步骤1中，构建符合变量内部相邻成员间变化规律的约束机制步骤为：子步骤1.1，通过分析水库库区地形的变化规律，提出水库水位～面积函数的微分特性为：设水库水位变量z，对应水面面积函数a＝f(z)，根据库区地形的一般变化规律可知，随着水库水位的抬高，相应水面面积越大，即面积函数的一阶导数大于零，有：a′＝f′(z)＝da/dz＞0式(1)另外，根据水位越高水面面积增长率也越来越大，即面积函数的二阶导数也大于零，因此有：a″＝f″(z)＝d2a/dz2＞0式(2)式2可以看出，面积函数为严格凹函数，由其性质可知，对于x1＜x2，存在0＜α＜1，有：f(αx1+(1-α)x2)＜αf(x1)+(1-α)f(x2)式(3)子步骤1.2，基于子步骤1.1的面积函数微分特性，等差水位相邻离散点距的相互制约关系为：设库容曲线中离散水位值为zj，(j＝0，1，2，……，n)，从小到大排列的序列，对应水面面积系列为aj＝f(zj)，由式(1)有：(aj-aj-1)/(zj-zj-1)＞0，因zj＞zj-1，有aj-aj-1＞0，即：0≤aj-1＜aj(j＝1，2，……，n)式(4)由式(2)有：对等差水位序列，简化为：子步骤1.3，针对某些库区的水位～面积曲线存在面积函数局部区间为非凹函数时，约束条件变更判别式及非凹函数情况下的制约关系为：对于遵循凹函数规律的水位～面积曲线而言，式(5)是严格成立的，有(2×aj-1-aj-2)＜(aj+1+aj-1)/2，整理得：aj+1＞(3×aj-1-2×aj-2)(j＝2，3，......，n-1)式(6)因此，式(6)是判别aj有无满足凹函数性质解的条件；对于部分特殊位置的地形，存在不满足凹函数的情况，即式(5)无解或式(6)不成立，所以式(6)也是约束条件变更的判别式，若式(6)不成立，说明在该库区高程范围对应的面积增长率非正值，此时的面积函数应变更为非凹函数约束，即：(aj+1+aj-1)/2≤aj＜aj+1(j＝1，2，......，n-1)式(7)子步骤1.4，初始化水位～面积曲线时，由于aj+1未知，无法利用式(5)～式(7)来约束aj，根据式(3)，其离散表达形式为：aj＜a0+(an-a0)×(zj-z0)/(zn-z0)式(8)对等差水位序列，进一步简化为：aj＜a0+(an-a0)×j/n(j＝1，2，……，n-1)式(9)由此，在初始化水位～面积曲线时，利用式(8)或式(9)来代替式5中的约束条件上边界。所述步骤2中，所述目标水库的运行资料包括库水位、入库流量、出库流量；所述气象资料包括库区降雨量、蒸发量。本发明优点体现在以下方面：1、本发明提出以不同水位对应的水面面积作为决策变量，并作用于符合实际库区地形变化规律的约束机制，计算结果更为合理，解决了常规水量平衡法求解结果中可能出现的高水位水面面积小于低水位水面面积这种违背常理的问题。2、本发明提出以不同水位对应的水面面积作为决策变量，解决了常规水量平衡法计算过程中无法精确计算时段蒸发量和渗漏量的问题，在一定程度上有效保证了库容曲线求解精度。3、本发明通过运行资料的预处理，分离了部分水位和部分库容，阻断了误差累计和传播途径，将误差限定到一定影响范围内；通过筛选水库蓄水和消落这种对称过程的资料作为输入，有效降低了流量资料中系统误差对库容曲线计算结果的影响，进一步提高了库容曲线求解精度。附图说明图1是本发明方法的流程框图。具体实施方式下面结合附图对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述实施例。如图1所示，本发明所述基于约束机制粒子群算法的水库库容曲线修正方法，包括以下步骤：步骤1，构建决策变量约束机制本发明以水库不同水位对应的水面面积作为决策变量，因此满足决策变量内部相邻成员间的约束条件是确保求解结果符合库容曲线自然变化特性的基础和前提。本发明通过分析水库库区地形的一般变化规律，总结提出水库水位～面积函数的微分特性，考虑实际应用时水库水位～面积～库容关系通常为离散值，提出等差水位条件下相邻离散点距的相互制约关系；具体包括以下子步骤：子步骤1.1，通过分析水库库区地形的一般变化规律，总结提出水库水位～面积函数的微分特性；设水库水位变量z，对应水面面积函数a＝f(z)，根据库区地形的一般变化规律可知，随着水库水位的抬高，相应水面面积越大，即面积函数的一阶导数大于零，有：a′＝f′(z)＝da/dz＞0式(1)另外，地形越向高处通常也越开阔，也就是说水位越高，水面面积增长率也越来越大，即面积函数的二阶导数也大于零，因此有：a″＝f″(z)＝d2a/dz2＞0式(2)式2可以看出，面积函数为严格凹函数，由其性质可知，对于x1＜x2，存在0＜α＜1，有：f(αx1+(1-α)x2)＜αf(x1)+(1-α)f(x2)式(3)子步骤1.2，基于子步骤1.1的面积函数微分特性，提出一般情况下等差水位相邻离散点距的相互制约关系；设库容曲线中离散水位值为zj(j＝0，1，2，……，n)，为从小到大排列的序列，对应水面面积系列为aj＝f(zj)，由式(1)有：(aj-aj-1)/(zj-zj-1)＞0因zj＞zj-1，有aj-aj-1＞0，即：0≤aj-1＜aj(j＝1，2，……，n)式(4)由式(2)有：对等差水位序列，可简化为：子步骤1.3，针对某些库区的水位～面积曲线可能存在局部“反常”情况，即面积函数局部区间为非凹函数时，提出约束条件变更判别式及非凹函数情况下的制约关系；对于遵循凹函数规律的水位～面积曲线而言，式(5)是严格成立的，有(2×aj-1-aj-2)＜(aj+1+aj-1)/2，整理得：aj+1＞(3×aj-1-2×aj-2)(j＝2，3，......，n-1)式(6)因此，式(6)是判别aj有无满足凹函数性质解的条件；对于部分特殊位置的地形，存在不满足凹函数的情况，即式(5)无解或式(6)不成立，所以式(6)也是约束条件变更的判别式，若式(6)不成立，说明在该库区高程范围对应的面积增长率非正值，此时的面积函数应变更为非凹函数约束，即：(aj+1+aj-1)/2≤aj＜aj+1(j＝1，2，......，n-1)式(7)子步骤1.4，初始化水位～面积曲线时，由于aj+1未知，无法利用式(5)～式(7)来约束aj；根据式(3)，其离散表达形式为：aj＜a0+(an-a0)×(zj-z0)/(zn-z0)式(8)对等差水位序列，可进一步简化为：aj＜a0+(an-a0)×j/n(j＝1，2，......，n-1)式(9)由此，在初始化水位～面积曲线时，可利用式(8)或式(9)来代替式5中的约束条件上边界；步骤2，水库运行资料预处理为便于分析和比较求解成果与运行资料的匹配程度，根据目标水库的运行资料(库水位、入库流量、出库流量)、气象资料(库区降雨量、蒸发量)、库区地质资料及近期的库容曲线，由水量平衡原理，逐时段计算初、末水位及其对应的时段库容差，如式(10)所示，并以此系列作为基准库容差；的计算公式如下：其中：k为当前计算时段标识，k＝1，2，......，p，p为总时段数；为由水量平衡方程计算的第k时段库蓄水量的变化，即时段基准库容差；分别为第k时段初、末的水库蓄水量；qs，k、qe，k分别为第k时段初、末的入库流量；qs，k、qe，k分别为第k时段初、末的出库流量；δtk为第k时段的时段长度，时段划分时应保证时段内入库、出库流量呈线性变化；wpre，k、wevap，k、wseep，k分别为第k时段库区的降水量、蒸发量及渗漏量；步骤3，目标函数设计首先，针对某一备选水位～面积曲线，根据椎体体积计算公式，计算与之对应的库容曲线，计算公式如下：其中：c、j均为库容曲线离散点序列标识，c，j＝1，2，......，d，d为离散点总数；zc、ac、zj分别为第c个序列的水位值、面积值和第j个序列的库容值；再从该库容曲线上逐时段查得运行资料中初、末水位及其对应的时段库容差值dvk，以此系列作为备选方案库容差；dvk的计算公式如下：dvk＝ve，k-vs，k式(12)其中：dvk为从备选库容曲线上查得的第k时段库蓄水量变化，即时段备选方案库容差；vs，k、ve，k分别为从式(11)计算的备选库容曲线上查得的第k时段初、末的水库蓄水量；以时段备选方案库容差dvk与基准库容差的离差平方和最小为本发明的目标函数，即：f＝min(φ(dv))式(13)其中：f为目标函数；p为运行资料的时段数量；φ(dv)为备选方案库容差dvk与基准库容差的离差平方和函数；min()为最小值函数；步骤4，基于约束机制改进粒子群算法的目标函数求解具体包括以下子步骤：子步骤4.1，设定算法基本参数，初始粒子数量m、粒子维度n、更新迭代周期上限s的具体取值应视求解问题而定，学习因子c1、c2均取2，速度限制常数vmax取值不宜过大，建议取值vmax＝1/(2n)；惯性权重w采用典型的线性递减自适应调整策略，计算公式如下：w＝wmax-(wmax-wmin)×s/s式(15)其中：wmax、wmin分别为设定的最大、最小惯性权重；s为当前迭代周期；子步骤4.2，建立粒子与决策变量间的映射关系；设第i个粒子xi＝(xi，0，xi，1，…，xi，j，…，xi，n)与决策变量ai＝(ai，0，ai，1，…，ai，j，…，ai，n)的映射关系如下：ai，j＝a0+(an-a0)×xi，j式(16)其中：i＝1，2，......，m，xi，j∈[0，1]；式(16)反函数为：xi，j＝(ai，j-a0)/(an-a0)式(17)子步骤4.3，基于步骤1的相关约束机制(式(5)及式(8)或式(9))，随机生成满足约束要求的初始粒子群体；初始化的粒子内各维度变量间应满足如下关系：其中：为第i个粒子初始化时(迭代周期为0)的第j个维度变量；因此，建议通过如下方式生成第i个粒子第j个维度的初始位置和速度：其中：sgn为符号函数，取值为{-1，0，1}；rnd为随机函数，且rnd∈[0，1]；子步骤4.4，逐个粒子i计算相应的适应度得到迭代周期s时各粒子的自身最好位置及群体最好位置对每个粒子，首先根据式(16)和式(11)转换成对应的备选库容曲线，再由式(12)、(14)及其他数据计算对应的以备选方案库容差，以备选方案库容差与基准库容差的离差平方和作为该粒子的适应度值。因此，迭代周期s时各粒子自身最好位置对应的适应度为：其中：min()为最小值函数；。迭代周期s时群体最好位置对应的适应度为：子步骤4.5，根据改进粒子群算法粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置；改进粒子群算法粒子的速度和位置更新公式如下：其中：r1、r2为介于[0，1]区间的随机数；根据式4、式5和式16，更新后的粒子各相邻维度位置之间应满足如下两个公式条件：根据式6和式16，更新后的粒子当前维度位置存在满足凹函数解区间的判别公式如下：式(23)用于确保水位～面积曲线的合理性，须强制性满足；式(24)用于确保水位～面积曲线符合凹函数要求，存在解空间时应强制性满足；式(25)是粒子当前维度是否存在凹函数解空间的判定条件；因此，粒子当前维度位置约束过程如下：式(23)成立时，判断式(24)是否成立：若式(24)成立，说明粒子当前维度位置满足要求，可转入下一步骤；若式(24)不成立，需进一步判断式(25)是否成立：若式(25)成立，说明粒子当前维度存在凹函数解空间，当j＝1或j＝n时，采用下式重新生成粒子当前维度的位置和速度：当j＝2，3，......，n-1时，采用下式重新生成粒子当前维度的位置和速度：若式(25)不成立，说明粒子当前维度不存在凹函数解空间，此时应变更粒子维度位置的约束机制，当j＝1或j＝n时，采用式(26)重新生成粒子当前维度的位置和速度；当j＝2，3，......，n-1时，采用下式重新生成粒子当前维度的位置和速度：式(23)不成立时，判断式(25)是否成立：若式(25)成立，说明粒子当前维度存在凹函数解空间，采用式(26)或式(27)重新生成粒子当前维度的位置和速度；若式(25)不成立，说明粒子当前维度不存在凹函数解空间，此时应变更粒子维度位置的约束机制，采用式(26)或式(28)重新生成粒子当前维度的位置和速度；子步骤4.6，重复子步骤4.4、4.5，直至粒子群更新迭代周期至上限s时，群体最好位置g(s)经式(16)和式(11)转换后即为求解结果；步骤5，输出条件判别设定允许距离阈值δ*，计算步骤4中求解的库容曲线结果与步骤2水库运行资料中的库容曲线的欧氏距离δ，计算公式如下：若δ≤δ*，则输出步骤4库容曲线结果；反之，用步骤4的库容曲线替代步骤2中近期库容曲线，并重复步骤2、3、4，直至满足输出条件。本发明方法各步骤说明如下：步骤1说明：1)水库库容曲线的离散格式水库库容曲线在理论上应是连续函数，即a＝f(z)和v＝g(z)在实际应用当中，由于这种函数关系难以精确获取，因此通常情况下，水库库容曲线均采用离散格式来标识，如表1所示。表1某水库库容曲线变量名称zav序号水位(m)面积(km2)库容(百万m3)01100.0011203.31121307.362314014.3169415025.4365表中，z0＝110，z1＝120，z2＝130，……；a0＝0.0，a1＝3.3，a2＝7.3，……v0＝0，v1＝11，v2＝62，……这种离散形式也满足a1＝f(z1)，v1＝g(z1)，a2＝f(z2)，v2＝g(z2)，……2)水库库容曲线具备性质的依据水库库容曲线所具备的函数性质，主要依据天然河谷的地形特性：a、水都是从高处向低处流，河道都是正比降；b、由于风化作用和山体稳定的要求，随着高度的增加，多数山体的自然坡度会先增加再变小。在上述两方面地形特性的作用下，反映到库容曲线形式上是水位～面积函数f(z)的1阶导数大于零，2阶导数决大多数情况下大于零。由此根据微分函数的导数特性，推导出离散形式下的面积变量aj(j＝0，1，2，……，n)的相互制约关系，包括式(4)、(5)和式(9)。例如上述库容曲线中，a2＞a1，即7.3>3.3。以及(2×a1-a0)＜a2＜(a3+a1)/2，即(2×3.3-0)<7.3<(14.3+3.3)/2。以上均用于说明对面积离散点aj，与其邻近点之间存在本专利提出的这种约束关系，为建模寻优提供了理论和物理基础。步骤2说明：本步骤主要为资料预处理，核心为水量平衡原理，利用该原理求解分段库容，得到(时段初水位～时段末水位～库容差)这样一系列数组；不直接累加，可以避免误差的累积。需要说明的是这里的库容差是利用水量平衡原理计算出来的，并非从库容曲线查出来的。步骤3说明：本步骤主要是构建目标函数，若已知水库库容曲线，根据步骤2中的时段初、末水位，可以直接在曲线中查得库容差值，如果所有测量数据没有误差，那么步骤2与步骤3这2个库容差值应是等值的。但由于测量误差的存在，2个库容差值不相等，在此以离差平方和最小为目标函数，认为哪个库容曲线代入目标函数的离差平方和最小，即为目标库容曲线。步骤4说明：本步骤主要为目标函数求解；子步骤4.1：设定算法基本参数举例说明设置如下：粒子数量m＝500，更新迭代次数k＝40，粒子维度n＝7，wmax＝1.3，wmin＝0.35，vmax＝1/14，最低水位z0＝600m，对应面积a0＝0km2，最高水位z7＝628m，假定对应面积上限为a7＝50km2。z1～z7分别为604、608、612、616、620、624、628，等4m间距。求解目标为最优的一组a1～a7。子步骤4.2：建立粒子与决策变量间的映射关系，如说明公式所示。子步骤4.3：初始化粒子每1个粒子代表1条备选水位～面积曲线，首先初始化粒子位置和速度上标“0”表示第0次更新，下标i表示第i个粒子，下标j表示第j个维度。初始化时应满足式18的关系。例如为0～1/7之间的随机数，为之间的随机数，……，为之间的随机数。例如随机生成的初始化粒子位置：……初始化粒子速度为0。子步骤4.4：逐个粒子i计算相应的适应度根据式(16)和式(11)，转换为对应的备选水位～面积曲线：……由式(11)转换成对应库容曲线后，再根据(12)、(14)计算对应的目标函数值，作为粒子的适应度值。根据式(20)、式(21)得到各粒子自身最好位置(第0代为自身)，群体最好位置：g(0)＝(0，0.016，0.117，0.138，0.235，0.237，0.299，0.950)，对应适应度值为：子步骤4.5：根据标准pso算法粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置。根据式(22)，更新粒子的位置和速度，对于速度超限制，调整至边界上。对的每个维度的位置值，代入式(23)～(25)进行验证，均满足可直接采用该位置值，若不满足则采用子步骤5中的调整方式进行调整，得到符合约束要求的位置和初始化的速度值。如式(23)不成立，采用式(26)重新生成当前维度的位置和速度值；式(23)、(24)、(25)均成立，保留其位置和速度信息；式(23)、(25)成立，式(24)不成立，采用式(27)重新生成当前维度的位置和速度值；式(23)成立，式(24)、(25)均不成立，采用式(27)重新生成当前维度的位置和速度值；……经约束调整后，各维度位置和速度如下例所示：其余粒子更新和约束调整方式类似，不再赘述。子步骤4.6:重复步骤(4)、(5)，得到：第1代群体最好位置：g(1)＝(0，0.017，0.036，0.060，0.103，0.224，0.428，0.931)，对应适应度值为：第10代群体最好位置：g(10)＝(0，0.013，0.047，0.090，0.140，0.203，0.353，0.614)，对应适应度值为：第16代群体最好位置：g(16)＝(0，0.029，0.059，0.097，0.135，0.176，0.356，0.601)，对应适应度值为：继续更新迭代群体最好位置不变，说明算法的收敛速度很快。得到的库容曲线如下：水位(m)面积(km2)库容(百万m3)600006041.451.936082.9510.566124.8526.006166.7549.106208.880.1162417.8132.2662830.05226.90步骤5说明：分析步骤4得到的库容曲线与目标库容曲线(上一次计算结果)的欧式距离：若欧氏距离在设代的阈值范围内，说明此次计算的库容曲线满足要求，可直接输出作为结果；若在阈值范围之外，可用本次计算得到的库容曲线代替目标库容曲线，重复步骤2、3、4，直至满足输出条件为止。当前第1页12