3

基于韧性曲线的安全韧性城市模型研究

3.1基于统计数据回归分析的复杂系统韧性建模研究
3.1.1概念与模型
3.1.1.1韧性三要素

韧性研究多采用韧性曲线。传统的韧性曲线通过图形表达的方法，描述系统在外界冲击下的功能变化情况，如图31所示。


图31传统的韧性曲线表达形式

（改自： Francis R,Bekera B,2014）



在实际建模过程中，为了量化“韧性”这一概念，需要推导出韧性的表达公式，而体现“韧性”的作用过程为： 外界冲击→系统韧性→功能变化，因此，仅有功能曲线的韧性曲线对于韧性的量化还不充分。为了正确地量化表征“韧性”，在给出韧性的功能曲线时，需要给出相应的外界冲击表征曲线。即完整的韧性曲线形式为： 功能曲线与相应的冲击曲线，如图32所示。


图32完整的韧性曲线表征


为此，提出“韧性三要素”的概念。“韧性三要素”的内容见表31。


表31韧性三要素



要素定义单位

冲击表征外界冲击作用总和的物理量，为过程量冲击单位，为模型建立的基础量
韧性复杂系统在突发事件的冲击下，能够依靠稳定性或适应能力来维持自身的功能，并能有序、快速地从非正常状态下恢复甚至改进的能力导出量
功能衡量系统功能水平的物理量，即系统韧性指标的表征量，为状态量功能单位，模型建立的基础量


应当指出，在实际应用时，需要通过冲击作用于系统引起的功能改变进行标度。为了定量化研究“韧性”特征，需要将外界冲击固化为一个客观存在、可定义的物理量，并假设外界冲击拥有基本单位，进而推导出韧性的标度单位。


功能单位可通过可观测的物理量(人口、GDP等)进行衡量，并归一化。
3.1.1.2韧性物理量
在定义了韧性三要素之后，对于韧性系统的定量化，可以推导出一系列的韧性相关物理量。经整理后，列于表32。


表32韧性系统的相关物理量



名称定义单位性质

冲击量——P表征外界冲击作用总和的物理量冲击单位基础量，过程量
冲击强度——I表征冲击量瞬时程度的物理量，即

I=dPdt冲击单位/时间导出量，状态量
冲击深度——H系统功能改变的最大值，即

H=|ΔQmax|功能单位读出量
冲击时间——t1从冲击开始到冲击极值所用时间秒/月/年读出量
恢复时间——t2从冲击极值到系统恢复所用时间秒/月/年读出量
韧性周期——T从冲击开始到系统恢复所用时间，称为一个韧性的周期秒/月/年读出量
冲击时间——d从冲击开始到冲击结束所用时间，称为冲击时间秒/月/年读出量
功能——Q衡量系统功能水平的物理量，即系统韧性指标的表征量功能单位基础量，状态量

功能变化率——k表征不同阶段功能变化速度的物理量，推广到一般模型时，基本定义式为

k=dQdt功能单位/时间导出量，状态量

3.1.1.3韧性方程
1. 韧性方程与韧性度量
韧性的作用体现在外界冲击使系统功能变化能力(大小、变化率)的量度上，韧性本身是系统的一种复杂性特征，靠单一指标量化很难全面表征，韧性在定义中所表达的抵抗、恢复等也是不同性质的动力学特征，无法通过一个参数度量。因此，我们提出从系统自身的复杂性出发，基于系统功能变化的机理分析，推导可能具有多维度特征的韧性表征模型和度量指标。
考虑一个复杂系统，在外界冲击下的功能变化不是单一过程，至少涉及三个方面的效应： 系统功能对外界冲击的响应； 系统对于功能瞬时状态改变趋势的抵抗； 系统向稳定功能状态的恢复。从系统动力学的角度，包含上述三方面效应的系统功能变化过程可用如下微分系统描述： 


dPdt-R1dQdt+R2Q=R0d2Qdt2(31)
其中，方程左端首项表示外界冲击的瞬时作用，在韧性过程中表征系统受到外界冲击的影响； 括弧中第一项表征韧性过程中包含的“抵抗”效应，与系统功能的瞬时变化率成正比； 第二项表征韧性过程中包括的“恢复”效应，与系统功能相对原稳定状态的偏离量成正比； 方程右端表征系统功能对外界冲击的响应，表现为系统功能对时间的二阶导数与外界冲击强度成正比。
由此可见，方程(31)表征了外界冲击作用于系统韧性引起的系统动力学响应过程，我们定义为韧性方程； 上述各项正比关系的比例系数代表了系统韧性从不同方面表现出来的多种含义，即我们认为韧性的度量应当是多维度的，定义为韧性分量系数，包括： 
R0(韧性惯量)： 外界冲击、自身抵抗力综合作用于系统功能并使其变化的作用因数，即系统功能变化二阶时间导数的比例系数。
R1(抵抗韧性)： 抵抗系统功能发生变化的系数。该抵抗是指某一时刻，系统自发产生的、相对于现在状态改变的抵抗，因此作用于系统功能的变化率项dQ/dt，与系统本身功能水平无关。
R2(恢复韧性)： 使系统趋于恢复到稳定状态的变化系数。该趋势是指当系统功能水平偏离稳定状态值时，系统内部产生的、使系统向稳定值靠拢的趋势，因此作用于系统功能的改变值Q。此公式下，可通过归一化等方法，令系统稳定值为零。
各韧性分量的定义方式列于表33。


表33韧性的分量表达



韧性分量机理大小方向韧性分量系数

抵抗分量系统对功能状态改变的瞬时趋势产生抵抗作用正比于功能状态改变率
(改变越快，抵抗越强)总是抵抗功能状态改变的瞬时趋势抵抗韧性R1

(Resisting resilience)

恢复分量系统对功能状态变化的积累趋势(即相对于期望水平的偏离量)具有恢复能力正比于功能状态偏离量
(偏离越大，恢复越强)向期望值改变(抵抗功能的偏离方向)恢复韧性R2

(Restoring resilience)

惯性分量外部冲击及系统抵抗和恢复的共同作用，导致了功能状态变化趋势的改变正比于功能状态变化率的变化率一致韧性惯量R0

(Inertial resilience)


2. 定性讨论
当系统遭受骤然冲击时，冲击的特点表现为时间Δt极短而强度dP/dt极大，高强冲击在极短时间内的累积量通常为有限值，因此在Δt内对方程(31)积分，得到


P-R1Q+R2∫Qdt=R0dQdt(32)

考虑到∫Qdt为无穷小量，可忽略。这样有


P-R1Q=R0dQdt(33)

若给出冲击极大的条件，则上式进一步简化为


P=R0dQdt(34)
上式表明，在外界冲击极大、骤然的条件下，系统韧性主要表现为惯量性质，或者说韧性中的惯性分量起主导作用，由此导出的韧性惯量表达式为


R0=P×t1H(35)
该式的意义为： 外界冲击一定的条件下，系统韧性惯量与系统在冲击下抵抗的时间成正比，与系统在冲击下改变的变化最大值成反比。需要指出的是，式(35)给出了系统韧性惯量的一种定义式，而非决定式，不能认为系统韧性惯量与外界冲击量成正比、与冲击深度成反比。系统韧性惯量R0为复杂系统的固有属性，不随外界冲击的状态变化。在某些冲击的极端条件下，甚至可能改变系统固有属性，即系统韧性惯量可能作为冲击条件的函数出现，这种高度耦合的情况非常复杂且极端，我们只讨论一般的情况。
更一般的情况下，采用上述化简过程中的一步结论，由方程(33)，骤然冲击下，Q与dQ的值均为H，有


P=HR1+R0t1(36)

由式(36)给出的结论，在外界冲击P的作用下，冲击深度H的大小结果取决于韧性惯量R0和抵抗韧性R1的线性组合，R0的系数为1/t1。也就是说，抵抗韧性的作用效果与时间参量无关。若抵抗韧性R1足够大，则一定外界冲击下，冲击深度H基本由抵抗韧性决定。
当系统遭受的冲击为一变化性的持续冲击时，考虑dP/dt比P更有意义。定义冲击强度I=dP/dt，I的函数形式将对Q(t)曲线的预测具有十分重要的指导意义。
在外界冲击结束、抵抗韧性R1很小的情况下，恢复韧性项R2Q对系统的韧性起到决定性的作用。此时的式(31)化简为


R2Q=-R0d2Qdt2(37)

式(37)为标准的谐振型系统关系式。此时，韧性系统被描述为一个偏离平衡位置的弹簧振子，弹簧振子的周期按照式(37)可给出


t=2πR0R2(38)

可预测，系统将在1/4周期之后回到平衡位置(稳定状态)，即


t2-t1=π2R0R2(39)
在t1t2的情况下，有


T=π2R0R2(310)
若满足一定条件： 外界冲击骤然、抵抗韧性极小，则韧性系统的韧性周期T将只由固有韧性与恢复韧性的比值决定，与外界冲击大小无关。反之，对于韧性惯量、恢复韧性比值相同的两个系统，韧性周期T的不同反映了系统抵抗韧性的不同； 对于抵抗韧性可忽略的两个系统，韧性周期T的相同反映了系统内固有韧性与恢复韧性存在的相同比例关系。
此简化模型还指出，在外界冲击结束后，系统功能水平超过平衡位置存在可能。
3.1.2韧性曲线建模与分型研究
3.1.2.1模拟韧性曲线分类

在韧性方程（3.1）的基础上，进行典型韧性图像的模拟和分类研究。
1. 骤然冲击型曲线
骤然冲击型曲线，指冲击的发生时间远小于韧性周期，冲击过程相对于整个韧性过程极短情况下的功能曲线。从曲线上看，骤然冲击型曲线具有以下显著特征： 功能下降迅速，系统功能在很短时间内降低到极值，抵抗期极短； 恢复过程相对平缓，恢复时间与抵抗时间有接近的数量级，但恢复期稍长于抵抗期。骤然冲击型曲线的形状特征是具有狭长的“尖峰”。
以公式(31)为基础，模拟生成骤然冲击型曲线。给予一个骤然性的外部冲击P，初始变化率(k)为零，功能(Q)稳定值为1，韧性惯量R0=10，抵抗韧性R1=0.1，恢复韧性R2=0.0001，图像时间范围为0~10000，外界冲击强度0.0005，持续时间5001~5200，冲击时间比e=d/T=200/20000=1%。得到韧性曲线见图33。


图33骤然冲击模拟


2. 持续冲击型曲线
持续冲击型曲线，指外界冲击持续时间与韧性周期具有相近数量级，冲击对于整个韧性过程都有显著影响情况下的功能曲线。
在模拟持续冲击型曲线时，将矩形脉冲的时间扩充，即可得到区别于图33的持续冲击型曲线，如图34所示。给予一个持续性的外部冲击P，初始变化率(k)为零，功能(Q)稳定值为1，韧性惯量R0=10，抵抗韧性R1=0.1，恢复韧性R2=0.00045，图像时间范围为0~1200，外界冲击强度0.0005，持续时间101~600，冲击时间比e=500/1200=42%。


图34持续冲击模拟


从外观上看，持续冲击型曲线与骤然冲击型曲线的区别在于，持续冲击型曲线具有更为平缓的抵抗期，这种“平缓”是相对于恢复期而言的。调节持续冲击型曲线的各分量，发现抵抗期的平缓是受到了恢复韧性R2的控制作用，当R2略微增大时，可以得到更平缓的抵抗期曲线。
3. 背景变化型曲线
背景变化型曲线，指无外界冲击时，系统功能本身存在变化的趋势情况下的曲线。需要指出的是，根据背景变化进行分类的分类方式，与根据冲击情况进行分类的分类方式不是并列的，是一种交叉分类的方式。在研究了最为简单的功能稳定值不变的情况后，需要考量系统功能本身存在的变化趋势在冲击过程中的影响。
考虑给前文中的骤然冲击添加一个强度为0.0000005的自然增长(约为外界冲击的0.1%)，增长方式为增长系数×时间，即


自然增长=0.0000005×t(311)

初始变化率k为零，功能Q稳定值为1，韧性惯量R0=10，抵抗韧性R1=0.1，恢复韧性R2=0.0001，图像时间范围为0~10000，外界冲击强度0.0005，持续时间5001~5500，冲击时间比e=500/20000=2.5%。得到图35。



图35背景变化型曲线模拟


分析背景变化对曲线形状的影响，由于背景增长与功能恢复同向，因此背景变化使曲线的抵抗期变得平缓，恢复期变得短促； 同时，背景增长也使得功能稳定值增大，这意味着功能在恢复到冲击前水平时，系统的恢复期尚未结束，功能仍然需要向更高的水平恢复，变相延长了恢复期。这让恢复期从视觉效果上反而变得更为“平缓”。
4. 波动型曲线
波动型曲线具有的特征为： 恢复期相对于抵抗期更加短促，恢复时间稍短于抵抗时间，存在明显的反复期曲线，并且有显著的“二次冲击峰”。微分公式(31)中的一阶导数项、二阶导数项可能对波动过程产生影响。改变各韧性分量的数量级关系，得到了呈现波动状的曲线结果，见图36。


图36波动型曲线模拟


图36的初始变化率k为零，功能Q稳定值为1，韧性惯量R0=100，抵抗韧性R1=0.01，恢复韧性R2=0.00001，图像时间范围为1~100000，外界冲击强度0.001，持续时间5001~10000，冲击时间比e=5000/100000=5%。可以看到，在此情况下，功能曲线出现了明显的恢复期波动。
3.1.2.2韧性曲线分型标准
在以上内容中，我们介绍了不同类型的韧性功能曲线，它们在形状外观上的区别是分类的重要依据。那么，如何将这种依据进一步量化为指标？在冲击归一化的条件下，功能曲线的形状实际是受到各韧性分量和冲击时间控制的，为了给出确切的量化标准，对冲击时间d、韧性惯量R0、抵抗韧性R1和恢复韧性R2的内在关系进行研究，定义

C1=R2d（312a）

C2=3R1d（312b）

C3=4R0d(312c)

研究发现，韧性曲线的形状实际受到C1、C2、C3的控制，我们将这三个定义数称为韧性曲线的控制数。控制数的意义表明，当冲击持续时间d扩大s倍时，为了保持图形的形状不变，R0、R1、R2应该分别扩大s2、s3、s4倍，即

R1/40∶R1/31∶R1/22=C3∶C2∶C1(313)
式(313)称为韧性曲线的形状控制式，该式的提出，对于韧性曲线的研究具有十分重要的指导意义，给予任意一条曲线准确的模拟基础。应用该公式分类前述曲线，对于骤然冲击型曲线，其形状特点为尖峰型，抵抗期短促，恢复期平缓，其控制系数的关系为


C3∶C2∶C1≈177∶46∶1(314)

对于持续冲击型曲线，其形状特点为对称型，抵抗期与恢复期具有对称的特点，其控制系数的关系为


C3∶C2∶C1≈177∶46∶2(315)

对于波动型曲线，其形状特点为反复期具有明显的波动，其控制系数的关系为


C3∶C2∶C1≈17.7∶4.6∶1(316)

对比几类韧性曲线分型标准如表34。


表34韧性曲线分型对比



曲线分型模 拟 曲 线冲 击 特 征韧 性 范 围

骤然冲击型
冲击骤然，冲击时间比e=1%C3∶C2∶C1≈177∶46∶1
持续冲击型
冲击持续，冲击时间比e=30%C3∶C2∶C1≈177∶46∶2
续表



曲线分型模 拟 曲 线冲 击 特 征韧 性 范 围

背景变化型
冲击骤然，冲击时间比e=2.5%C3∶C2∶C1≈177∶46∶1
波动型
冲击持续，冲击时间比e=12.5%C3∶C2∶C1≈17.7∶4.6∶1


3.1.2.3实际韧性曲线拟合
以式(31)为基础，通过多元回归拟合实例曲线的韧性分量(R0、R1、R2)，再以韧性分量作为参数，对实例曲线进行再现。
1. 拟合理论
在多元回归拟合韧性分量时，需要经历的过程为： 归一化——矩阵化——多元回归。以公式表达，即



dPdt-R1dQdt+R2Q=R0d2Qdt2

fQ，dQdt，d2Qdt2，dPdt=0

Q，dQdt，d2Qdt2R=dPdt(317)

R矩阵即为所求韧性的表示矩阵，矩阵元素即为韧性分量R0、R1、R2。
2. 拟合过程
以埃博拉病毒在塞拉利昂的感染人数为例，实际的感染曲线如图37所示。



图37埃博拉感染实例(数据来源： 世界卫生组织WHO周报)


拟合之前首先要对韧性曲线做分型控制的处理。韧性曲线拟合的瓶颈在于，外界冲击的时间函数是未知的。在模型研究中，我们将外界冲击都处理为有限区间的矩形函数，则由冲击强度和两个参数唯一确定。根据功能曲线，冲击强度采用归一化处理，持续时间参考冲击深度对应的时刻，并利用控制数，估算拟合参数的大致范围，确保冲击曲线的分型正确，再进行严格的定量拟合。
(1) 归一化处理。令纵轴取值范围为0~1，通过实际病例数/总病例数，得到归一化后的感染比例。
(2) 初步模拟。通过前文中所述的控制原理，对需拟合曲线进行尝试性模拟。根据原始图像的冲击时间d、韧性周期T，模拟出具有近似形状与周期的模拟图像。在形状吻合的情况下，根据模拟的结果，确定变化率k、惯量n的数量级。
(3) 数据扩充。系统样本数据为T=76，为了更好地拟合结果，将数据复制1000次，扩充为76000的样本，再通过Origin软件的soomth功能(每20000个数据进行回归)，得到平滑、渐变的新数组Q。
(4) 背景消除。对于简单的背景增长或衰减，需要在拟合之前进行背景消除工作，得到稳定状态为不变量的Q值曲线。
(5) 微分处理。对Q值曲线进行一次微分、二次微分，得到一阶导数数组k、二阶导数数组n。此时，需要评价步骤(3)中的soomth过程，如果k、n的数量级符合模拟结果，则认为平滑处理过程是合理的； 反之，根据数量级的不同，需要改变soomth的回归步长。
(6) 外界冲击假设。根据实际Q值曲线，确定韧性曲线分类为持续性冲击。从图37中读出t1=18天，故根据冲击开始与结束的时间，设置外界冲击矩阵I为


I=I1
I2
I3(318)

其中


I1=［111…1］T(5000项)

I2=［000…0］T(18657项)

I3=［111…1］T(52343项)(319)