基于系统微分响应的暴雨强度公式参数优化

2022-09-29李琼芳许树洪周正模和鹏飞虞美秀陈启慧

水资源保护 2022年5期

李琼芳，许树洪，周正模，和鹏飞，杜尧，虞美秀，陈启慧

(1.河海大学水文水资源学院，江苏南京 210098；2.长江保护与绿色发展研究院，江苏南京 210098)

高精度的暴雨强度公式是城市防洪排涝基础设施建设规划设计的重要依据，事关城市地区人民生命财产安全和社会经济的可持续发展[1]。在过去几十年的探索与实践中，已经形成了暴雨选样、频率分析、公式拟合的一套完善的暴雨公式编制流程。随着新的优化算法的应用，公式参数的求解发展出多种方法[2-4]，叶姗姗等[5]分别用粒子群算法、高斯-牛顿法、遗传算法、麦夸尔特法推导得到淮安市暴雨强度公式，并对比分析了不同算法得到的结果；白庆芹等[6]选用麦夸尔特全局优化法、准牛顿全局优化法、模拟退火法、粒子群算法及遗传算法求解得到暴雨强度公式参数；张小潭等[7]等基于改进的实时编码加速遗传算法(real coding-based accelerated genetic algorithm, RAGA)求解暴雨强度公式参数，并与传统求解方法(牛顿迭代法、高斯-牛顿法)和优化算法(RAGA、蚁群算法)得到的结果进行了对比。尽管不同学者在求解暴雨强度公式参数时使用的优化方法有所不同，但在求解过程中均需要设定参数收敛的目标函数，如给残差分配相等权重的绝对均方误差、不等权重的相对均方误差等，这类目标函数都与误差平方和目标函数的本质相同。使用此类目标函数，一方面会由于非线性模型的参数平方而出现额外解[8]；另一方面此类目标函数的响应面存在大量局部最优值，使智能算法很容易陷入局部最优，无法得到真正的全局最优解[9]。因此，智能算法求解非线性问题时会普遍出现参数等价性、模糊性、非唯一性、病态性和不可识别性现象[10]，如暴雨强度公式参数求解中的“异参同效”问题[11]。针对以上问题，亟须论证智能算法求解暴雨强度公式参数能否收敛到全局最优[12]，并探索在不考虑其他可行参数集的情况下能有效获取唯一最优参数集[13]的新的暴雨强度公式参数求解方法。

本文以国家海绵城市建设试点城市——镇江市为研究对象，解析暴雨强度公式参数求解常用目标函数响应面上参数解的空间分布特征，论证分析传统智能算法在暴雨强度公式参数求解中的效果，提出基于系统微分响应的暴雨强度公式参数优化方法并分析其应用效果。研究成果可为优化求解暴雨强度公式参数和提升暴雨强度公式精度提供新的思路。

1 常用暴雨强度公式参数求解方法存在的问题

1.1 参数求解常用目标函数合理性分析

目前国内常用的4参数城市暴雨强度公式为

(1)

式中：I为暴雨强度；A1为雨力参数；C为校正参数；T为重现期；t为降雨历时；b为修正参数；n为暴雨衰减指数。

《城市暴雨强度公式编制和设计暴雨雨型确定技术导则》(下文简称《导则》)中规定需要按编制公式计算暴雨强度并进行精度检验，不同于绝大多数非线性函数参数优化是以计算与实测样本的误差平方和为目标函数；GB 50014—2006《室外排水设计规范》(2011年版)(下文简称《规范》)则要求把计算值与实测值的平均绝对均方误差(AMSE)和平均相对均方误差(RMSE)作为评价暴雨强度公式精度的指标。因此，为使所得暴雨公式参数反映样本统计特征，将AMSE或RMSE作为参数优选目标函数。

以一般p次非线性函数为例：

y=x+ax+(ax)2+…+(ax)p

(2)

式中：y为模型输出；x为模型输入；a为参数。

L个样本的函数计算值与实测值的AMSE为

(3)

式中xi、yi分别为每组样本的模型输入和实测值。

将AMSE作为目标函数，对式(3)求一阶偏导，得到

(4)

令上式等于0时会额外增加了p-1个参数解，因此将AMSE作为目标函数，求解此目标函数最小值时会显著增加参数求解的不确定性。

假设暴雨强度公式参数值分别为A1=5、C=1、b=1、n=0.7，该暴雨强度公式可以表示为

(5)

以式(5)得到的暴雨强度为实测暴雨强度。假设存在4个样本，分别为(T=10 a，t=5 min)、(T=10 a，t=6 min)、(T=100 a，t=5 min)、(T=100 a，t=6 min)，并假定在参数C和b已知的条件下以AMSE为目标函数确定参数A1和n。假定A1在[0，10]内每间隔1取1个值，n在[0，1]内每间隔0.1取1个值，生成11×11的网格，计算网格点目标函数值，绘制得到网格化目标函数曲面,如图1(a)所示，可以看到局部最优解沿深色区域呈曲线延展至全局最优，局部最优解共7个。

(a)2个未知参数

类似地，假定参数C已知，通过寻找AMSE最小的参数组合来确定未知参数A1、b、n。假定A1和b都在[0，10]内每间隔1取1个值，n在[0，1]内每间隔0.1取1个值，生成11×11×11的网格，计算网格点AMSE，得到的网格化目标函数曲面如图1(b)所示，相应的局部最优解共21个。

随着参数个数的增加和参数组合样本的增多，目标函数曲面上的局部最优值数量越来越多，可以发现大量局部最优的目标函数值近似相同。因此可以推断，当参数范围越大，间隔步长越小(目标函数曲面越光滑)，参数组合越多，局部最优参数解的个数会趋于无穷多个，并且在参数真值附近的局部最优解的目标函数值会与参数真值的目标函数值近似相等，“异参同效”现象非常显著。

这些问题的存在表明暴雨强度公式参数求解常用的目标函数有缺陷，无法保证推求得到的参数解为真值。下面将结合具体案例进一步分析暴雨强度公式目标函数曲面上参数解的分布特征，并结合一些常用智能算法评价其求解能力。

1.2 基于具体案例的目标函数响应面参数解分布特征

1.2.1研究区选取及设计暴雨强度推求

选取国家海绵城市试点城市——镇江市作为研究区。相关研究表明镇江市短历时暴雨发生次数呈显著增长趋势[14]，因此编制高精度的城市暴雨强度公式对城市防洪至关重要。基于镇江市国家基本气象站丹徒站1981—2016年逐分钟实测雨量资料，采用年最大值选样方法构建5 min、10 min、15 min、20 min、30 min、45 min、60 min、90 min、120 min、150 min、180 min共11个历时的暴雨资料统计样本。考虑到P-Ⅲ型分布相较于耿贝尔分布和指数型分布更适用于研究区[15]，选用P-Ⅲ型分布对各历时年最大降水量样本进行频率分析，得到100 a、50 a、30 a、20 a、10 a、5 a、3 a、2 a共8个不同重现期不同历时的设计暴雨强度(表1)。

表1 不同历时不同重现期的暴雨强度

1.2.2目标函数响应面参数解分布特征

由于样本数量较大，目标函数较为复杂，很难通过对公式及目标函数的理论分析来获取参数解数量和位置信息，因此采用随机抽样的方法分析参数解的信息。根据镇江市暴雨强度公式现有的参数取值：A1=38.362 3、C=1.017 3、b=19.137 7和n=0.975，将各参数初始范围设置为A1∈(0,50)、C∈(0,1.1)、b∈(0,30)、n∈(0,1)，然后在参数范围内通过均匀分布生成10万组随机样本。根据《导则》要求编制的暴雨强度公式，其精度应满足AMSE小于0.05 mm/min，RMSE小于5%，因此分别统计满足精度要求的样本组数，并确定相应的参数取值范围，在此范围内继续进行抽样，重复上述步骤直至满足精度要求的参数范围不再变化为止，最后点绘出该参数范围内样本参数值-目标函数值关系图，并分析参数解的分布特征。

表2给出了均匀随机抽样结果。在没有有效确定参数初始值范围的情况下，10万组样本中仅有2组样本满足AMSE小于0.05 mm/min，即随机抽样法求解参数能满足精度要求的概率仅为0.000 02。经过第一次抽样，缩小了参数初始值范围，在第二次采用随机抽样法时，满足精度要求的概率则达到了0.108 69，有了明显提升，但满足精度要求的参数范围与前一次抽样确定的参数范围基本相同，说明能满足AMSE小于0.05 mm/min的稳定的参数范围已基本确定。想要通过进一步约束参数范围来获得最优解，需要加强约束条件。序号3给出了在第2次抽样确定的参数范围内随机抽样，统计满足AMSE小于0.045 mm/min和RMSE小于4%的样本个数及参数范围。可以看出满足AMSE小于0.045 mm/min的样本数为3 118，即满足精度要求的概率为0.031 18，但也发现满足精度要求的参数范围较初始范围没有较明显的变化。此外，案例还表明满足AMSE精度要求的样本均能满足RMSE的精度要求，因此下文均以AMSE作为目标函数进行分析。

表2 均匀随机抽样结果

图2为第3次抽样获得的10万组样本的AMSE与各参数的关系。从图2可以看出，约束参数范围后抽样所得到的样本中仍有一部分无法满足AMSE的精度要求；但满足精度要求的参数解数量众多，且不同参数样本的目标函数值近乎相同；同时满足AMSE最小的样本有无穷多个，“异参同效”特征非常明显，进一步证实了第1.1节的结论。此外，未发现能够决定目标函数最优解的参数，但相对而言，AMSE对参数A1、n较为敏感，其值的合理界定能够有效过滤导致目标函数值相对较大的无用样本。当A1∈(37,38)、n∈(0.94,0.95)时，AMSE变化区间明显缩窄，由样本信息可知，当A1=37.877 9、n=0.943 8时，AMSE最小，因此A1、n范围的有效界定对模型参数的优化求解尤为重要。通过上述抽样分析方法可以不断缩小参数范围，最终选取对应于样本中目标函数最小的参数组合作为最优解，但此方法相对繁琐，建议慎用。

(a)A1-AMSE散点图

从以上分析发现，在参数范围未有效确定的情况下，满足精度要求的有效样本数目较少，而在缩小参数范围后，满足精度要求的参数解数量又较多，信息混杂。因此利用随机搜索类优化算法推求暴雨强度公式参数无法保证所获得的参数值是其真值。

在最终确定的参数范围内生成多维度参数组合的离散化网格，并计算网格点目标函数值。将A1、C、b、n参数中任意两个作为x和y轴，目标函数作为z轴，绘制参数-目标函数值三维空间图(图3)；将A1、C、b、n参数中任意3个作为x、y和z轴，绘制参数-目标函数值切片图(图4)。基于三维空间图和切片图，分析参数变化对目标函数值的影响特征以及满足精度要求的参数解数量及分布情况。

(a)A1、C

(a)A1、C、b

从图3中可以看出，任意两参数组合对应的目标函数响应面呈现U形，满足精度要求的参数解分布在响应面底部深色区域，深色区域呈线性延展至全局最优点。响应面底部存在间隔分布的凹槽，凹槽底部为局部最优值。从图4中可以看出，目标函数值较小的参数组合样点分布在图中深色区域，存在无穷多个，在参数变化范围内均有分布。从其分布位置与参数的相互关系来看，A1、C、b参数中，目标函数值随A1变化而变化，而基本不随C、b变化而变化；A1、C、n参数中，参数A1、n对目标函数值影响程度大于参数C，尤其是当A1、n同时减小时，目标函数值也减小；A1、b、n参数中，A1、n对目标函数值影响程度大于参数b，目标函数值也是随着A1、n的减小而减小，但不随b变化而变化；C、b、n参数中，n对目标函数值的影响程度大于C、b，目标函数值随着n的减小而减小，但不随C、b变化而变化。总体来看，目标函数值对A1和n的敏感程度较高，此结论与采用Sobol法[16]的敏感性分析结果(表3)一致。

表3 Sobol法灵敏度分析结果

1.3 常用智能算法的暴雨强度公式参数求解能力

从以上分析结果来看，当采用随机搜索的智能算法求解暴雨强度公式参数时存在以下问题：①在参数范围没有得到有效界定的情况下，基于随机生成的参数个体/种群来筛选信息进而寻优，很难随机生成包含全局最优参数信息的种群；②即使在参数范围确定的情况下，由于局部最优值数量过多，当初始生成的参数种群不包含全局最优信息时，算法会陷入局部最优，避免陷入局部最优的办法是增加新的参数样本更新种群，但如果新的参数样本仍然不能包含全局最优信息，算法又会陷入下一个局部最优；③由于全局最优附近的局部最优和全局最优的目标函数差别太小，而数学软件能识别的小数位数有限，算法将局部最优等同于全局最优处理，最终得到的结果可能只是目标函数近似最小的局部最优参数解，即“异参同效”。

结合前面的实际案例，本文采用粒子群算法和SCE-UA(shuffled complex evolution algorithm)算法优化求解了镇江市暴雨强度公式参数。表4给出了在参数初始范围A1∈(0,50)、C∈(0,1)、b∈(0,30)和n∈(0,1)条件下两种方法寻优得到的最终参数解及目标函数值。可以看出粒子群算法寻优效果相对较好，5次参数寻优对应的目标函数均满足精度要求(<0.05 mm/min)，但每次参数解却不尽相同。需要指出的是，当目标函数最小值落在参数范围边界时，算法不能够自动折返，而会陷入局部寻优。SCE-UA算法在求解具有区间约束的非线性问题时较为有效[17]，但在参数范围没有得到有效约束的情况下，寻优效果相对较差。

表4 优化算法参数解及目标函数值

2 暴雨强度公式参数系统微分响应求解方法

上述分析表明，以AMSE/RMSE为目标函数，在目标函数响应面上进行参数寻优存在诸多问题。由于暴雨强度公式本身在参数邻域内连续且可求导，因此可通过样本(Xi,Yi)与参数函数y(Y)=f(θ,Xi)相交的位置获取参数真值信息，从而使得在参数函数曲面上寻优的系统微分响应方法[18]具有较大的应用价值。

假设非线性函数模型可以表示为

Y=f(θ,X)

(6)

其中

θ=(θ1,θ2,…,θq)T

X=(x1,x2,…,xk)T

式中：θ为模型参数；X为模型输入自变量。

在模型参数优化时，将模型参数看作自变量，对函数模型线性化，即一阶泰勒级数展开。其表达式为

(7)

其中

Y(j)=f(θ(j),X)

式中：j为迭代次数；θ(j)为参数θ第j次率定值；e(j)为观测值与计算值的偏差。

假定共有m组观测样本，即(X1,Y1)、(X2,Y2)、…、(Xm,Ym)，函数模型向量形式为

Y=Y(j)+S(θ-θ(j))+E(j)

(8)

其中

Y=(Y1,Y2,…,Ym)T

要使函数计算值与观测值误差最小，校正后的参数θ满足：

(9)

因此，式(8)变换为Y-Y(j)=S(θ-θ(j))，由此可得

θ(j+1)=θ(j)+(STS)-1ST(Y-Y(j))

(10)

通过上述过程可由初始参数值一步步迭代计算得到参数最优值。但是由于式(7)只是对实际函数的线性近似，忽略了高阶泰勒级数展开项。因此引入校正系数c(c>1)，将式(8)修正为

Y=Y(j)+cS(θ-θ(j))+E(j)

(11)

相应地，式(10)可变化为迭代式：

θ(j+1)=θ(j)+B(STS)-1ST(Y-Y(j))

(12)

式中0

(13)

计算迭代过程中止应同时满足两个原则：

a.函数收敛原则。当算法在一次或多次迭代中不能明显改善目标函数值时，迭代过程停止，即|Y(j+1)-Y(j)|<ε，其中ε为收敛性控制常数，ε≥0。

基于系统微分响应优化方法的暴雨强度公式参数优化求解步骤为：①参数A1、C、b、n的初始化，即随机生成一组参数初始值θ(0)；②将初始参数θ(0)与样本结合，计算灵敏度矩阵S以及每个样本的Y-Y(0)；③根据式(10)确定搜索方向，根据式(12)(13)确定参数B；④利用式(12)计算得到θ(1)；⑤判断是否满足收敛性原则，若满足则终止寻优，若不满足则令θ(1)=θ(0)转步骤②。

3 实例应用与讨论

3.1 实例应用

为进一步分析系统微分响应参数优化方法的应用效果，利用该方法推求镇江市暴雨强度公式。图5给出了按照第2节步骤进行参数优化求解的迭代收敛过程(不同颜色代表参数不同初始值的优化过程)。由图5可知，对于不同参数初始值，最后都能收敛到同一结果，即参数真值，且各参数迭代次数均在15次以内，表明算法优化效率较高。当参数初始值接近参数范围下边界时，在优化求解过程中，参数取值逐渐增加，逐渐逼近参数真值，且算法能够通过灵敏度矩阵自动缩小步长，逐渐收敛到真值；而当参数取值接近参数范围上边界时，经过寻优，参数取值会向下边界逼近，直至小于参数真值，再逐渐增大向真值逼近。由此可知，该方法更擅长于从参数取值范围下边界逼近参数真值，但当参数取值靠近上边界时，算法都能有效折返，可以避免陷入局部最优。

(a)A1优化过程

综上，暴雨强度公式参数线性化优化方法能够寻找到暴雨强度公式参数的真值，且迭代次数较少，效率较高，该算法不需要通过目标函数反推参数，且能够解决参数范围的确定对求解结果的影响问题，能够避免陷入局部最优。推求的5～180 min镇江市暴雨强度公式为

(14)

式(14)所得暴雨强度公式计算的AMSE为0.042 63 mm/min，RMSE为3.360 1%，精度均满足《规范》要求。

3.2 讨论

系统微分响应法的理论基础是系统对参数扰动的局部敏感性[19],该方法在优化求解效率和精度上明显优于粒子群、SCE-UA算法。该方法通过一阶泰勒级数展开将非线性模型转化为线性模型，同时利用最小二乘法求解线性逆问题，由暴雨强度的计算误差推求参数的估计误差，并以此进行参数修正，具有较为清晰的物理解释，结构较为简单，性能较为稳定。针对暴雨强度公式参数求解中遇到的“异参同效”难题，许多优化算法可能仅仅得到了满足精度要求的一种解，而非唯一解，而系统微分响应法可以高效获得暴雨强度公式参数真值，表现出了明显的优势。从推导过程来看，该方法应适用于其他地区暴雨强度公式参数的优选求解。此外，如何构造新的目标函数来提高参数优化方法的优化性能[20]，同时反映模型系统参数本身特性，这是参数优化领域的一大难题。在暴雨强度公式推求应用中，新目标函数的确定也具有研究价值，李增永等[21]在推求暴雨强度公式参数时，采用对数离差绝对值和平均对数离差绝对值作为目标函数，此目标函数反映出双对数坐标系中暴雨强度-历时关系曲线呈现阶段式的线性特征[22]，因此从理论上更能反映暴雨衰减指数随历时变化的特点，但此目标函数响应面的参数解数量、分布特征以及目标函数与参数的关系值得进一步研究论证。