,*,+,2), Michael Beer**FUNCTIONAL PERSPECTIVE OF UNCERTAINTY QUANTIFICATION FOR STOCHASTIC PARAMETRIC SYSTEMS AND GLOBAL SENSITIVITY ANALYSIS1)
Wan Zhiqiang*,+, Chen Jianbing,*,+,2), Michael Beer**通讯作者: 2) 陈建兵, 教授, 主要研究方向: 随机动力学与控制及结构可靠性. E-mail:chenjb@tongji.edu.cn
收稿日期:2020-09-21网络出版日期:2021-03-18
| 基金资助: |
Received:2020-09-21Online:2021-03-18
作者简介 About authors
摘要
不确定性因素广泛存在于实际工程分析与设计之中. 例如,土木工程结构的力学性能参数可能存在不可忽略的随机性.故而随机系统参数灵敏度分析, 对合理的工程设计与决策具有相当重要的意义.本文首先论述了随机系统中不确定性量化与传播的泛函空间分析观点. 在此基础上,给出了由泛函Fréchet导数所定义的整体灵敏度指标,讨论了该整体灵敏度指标的一些基本数学物理性质,定义了该整体灵敏度指标所对应的工程实用的重要性测度与方向灵敏度,并分别解释了二者的几何物理意义. 然后, 具体论述了在$\varepsilon$-等效分布定义下该整体灵敏度指标的参数化表达形式.结合概率密度演化-概率测度变换方法,给出了该指标的具体数值求解算法及整体灵敏度分析流程. 通过4个算例分析,包括正态随机变量线性组合解析函数分析、隧道锚杆支护系统稳定性分析、大坝下稳态受限渗流量分析以及钢筋混凝土平面框架结构随机地震响应分析,具体说明了该指标的实用性及重要性.
关键词：
Abstract
Uncertainty exists broadly in real engineering design and analysis. For instance, some mechanical parameters of structures in civil engineering may be of randomness and usually cannot be ignored. Therefore, the process of uncertainty quantification, e.g., the sensitivity analysis on parameters of stochastic systems is, of paramount significance to reasonable engineering design and decision-making. In the present paper, the perspective of functional space analysis on uncertainty quantification and propagation in stochastic systems is firstly stated. On this basis, the global sensitivity index (GSI) is introduced based on the functional Fréchet derivative, of which some basically mathematical and physical properties are studied. Besides, the correspondingly defined importance measure and direction sensitivity of the GSI are also discussed, in terms of their geometric and physical meanings. Moreover, based on the definition of $\varepsilon$-equivalent distribution, the parametric form of the proposed GSI is elaborated in detail. By incorporating the probability density evolution method (PDEM) and the change of probability measure (COM), the numerical algorithm of the GSI and the procedure of sensitivity analysis are illustrated. Four numerical examples, including the analytical function of the linear combination of normal random variables, stability analysis of the rock bolting system of tunnel, the analysis of steady-state confined seepage below the dam, and the stochastic structural analysis of the reinforced concrete frame, are analyzed to demonstrate the effectiveness and significance of the GSI.
Keywords：
PDF (2220KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
万志强, 陈建兵, Michael Beer. 随机参数系统不确定性量化的泛函观点与整体灵敏度分析1). 力学学报[J], 2021, 53(3): 837-854 DOI:10.6052/0459-1879-20-336
Wan Zhiqiang, Chen Jianbing, Michael Beer.
引言
现代化土木工程建设正逐步走向整体化、精细化及全概率化的全新分析与设计范式. 与之相应的是, 数值仿真计算模型已然成为工程结构系统设计与分析中不可或缺的一部分. 然而, 设计与分析过程中不可避免地会遇到诸多不确定性因素的影响, 如仿真计算模型的边界约束可能具有不确定性、初始条件可能难以精确给定、系统的某些力学性能参数可能具有不可忽略的随机性等[1-2].一般而言, 这些不确定性因素均可表征为某种参数化的形式,如混凝土材料的抗压强度与弹性模量[3]、岩土材料的有效黏聚力与有效内摩擦角[4]、随机风场模型的地面粗糙度与剪切速度[5]等.由于系统输入参数存在不确定性, 必然导致系统输出亦存在一定程度的不确定性,因而仿真预测的结果不可避免地与实际情况存在一定偏离. 研究表明,不同系统输入参数对系统输出不确定性的“贡献”一般是不同的. 因此,量化系统输入参数对系统输出不确定性的影响(或称量化系统输出关于系统输入参数的灵敏度),即灵敏度分析, 对工程设计、优化与决策具有重要意义[6].
灵敏度分析的目标是量化系统输入参数对系统输出不确定性的影响,即选择并计算灵敏度指标. 例如, 基于差商的基本效应指标[7] (亦称为Morris指标)及基于导数的灵敏度指标[8], 主要适用于确定性系统,在本质上属于局部灵敏度的范畴. 另一方面, 针对随机参数系统,人们发展了基于方差分解的Sobol'指标[9]以及基于条件概率密度函数的矩独立指标[10]等整体灵敏度指标.这些指标侧重于对灵敏度“大小”的计算, 因而往往其值恒为非负.这容易使人产生输入变量的不确定性总是使得输出变量的不确定性增大的错觉[11].事实上, 系统输入参数对输出不确定性的影响是有“方向”的,而这一“方向”对结构分析、特别是考虑不确定性影响下的结构优化设计至关重要[12].例如, 基于失效概率的灵敏度指标便是一类能同时反映灵敏度“大小”与“方向”的指标[13-14],被广泛用于基于可靠度的结构优化与设计之中[12-15].
在上述研究工作的基础上, Chen等[16]提出了一类全新的基于Fréchet导数的整体灵敏度指标,并给出了相应的高效算法. 其中, 随机系统中不确定性传播的泛函观点虽已初具雏形,但尚未得以深入论述. 同时, 如何针对具体的工程问题采用基于Fréchet导数的整体灵敏度指标,特别是对参数进行灵敏度排序等相关问题的研究尚不够充分.
基于此, 本文从随机系统中不确定性传播的泛函表达角度,进一步探讨了该整体灵敏度指标的数学定义、物理意义和基本性质. 进而,在$\varepsilon $-等效分布的定义下, 讨论了该整体灵敏度指标的参数化表达形式,给出了其对应的重要性测度和方向灵敏度的定义及其物理意义.建立了该整体灵敏度指标与矩独立指标之间的内在联系.结合概率密度演化-概率测度变换方法[17],提出了该整体灵敏度指标及其对应的重要性测度与方向灵敏度的数值计算方法和流程.最后, 通过4则算例分析详细讨论了该灵敏度指标的工程实用性.
1 随机系统中不确定性量化与传播的泛函观点
1.1 随机变量变换的泛函空间观点
在概率论框架下, 随机系统中输入变量的不确定性可以随机变量(向量)的数学形式进行表征[18-19].记概率空间为$( \varOmega ,\mathfrak{F}$, $\mathbb{P}$), 其中$\varOmega$为样本空间、$\mathfrak{F}$为$\sigma $-代数($\varOmega $的事件域)、$\mathbb{P}$为概率测度, $\varpi \in \varOmega $代表样本空间中的一个点(即基本事件).不失一般性, 考虑一随机系统其中, ${g}(\cdot ):\mathbb{R}^{n}\to \mathbb{R}^{m}$为一连续映射, 其分量为$( g_{1}$, $g_{2}$, $\cdots,g_{m})^{T}$; ${\boldsymbol \varTheta}(\varpi)=\left( {\varTheta_{1}(\varpi),\varTheta_{2} (\varpi),\cdots ,\varTheta_{n}(\varpi)} \right)^{T}$为$n$维输入随机向量, 其联合概率密度函数$p_{{\boldsymbol \varTheta }} (\boldsymbol \theta)$已知, 其中${\boldsymbol \theta }=\left( {\theta_{1} ,\theta_{2} ,\cdots ,\theta_{n} } \right)^{T}$为随机向量${\boldsymbol \varTheta}$的可实现值; ${X}(\varpi)=\left( {X_{1} (\varpi),X_{2} (\varpi),\cdots ,X_{m} (\varpi)}\right)^{T}$为$m$维系统输出向量, 记其联合概率密度函数为$p_{{X}}\left( {{x}} \right)$, 其中${x}=\left( {x_{1} ,x_{2} ,\cdots,x_{m} } \right)^{T}$为随机向量${X}$的可实现值.
由概率论可知, 若已知$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta }}\right)$, 则根据上式可确定输出随机向量的概率密度函数$p_{{X}} \left({{x}} \right)$. 事实上,根据概率守恒原理[20-21], 可得
其中$\varOmega_{{\boldsymbol \varTheta }} $为输入随机向量值域空间的任意子域, $\varOmega_{{X}} $为由式(1)中的映射决定的在输出随机向量值域空间中与$\varOmega_{{\boldsymbol \varTheta }} $对应的子域. 由此可得
为简单计, 先考察$n=m=1$的情况. 此时, 式(1)成为标量形式$X(\varpi)=g\left( {\varTheta (\varpi)} \right)$,而式(3)则相应地退化为$\int_{\varOmega_{X} } {p_{X} \left( x \right){d}x}=\int_{\varOmega_{\varTheta } } {p_{\varTheta } \left( \theta \right){d}\theta} $. 注意此时积分区域为$\varOmega_{X} =g\left( {\varOmega_{\varTheta } } \right)$.进一步, 考虑下列情况.
(1) 函数$g(\cdot )$为单调函数. 此时, 有
其中$\lambda (\cdot )$表示子域的Lebesgue测度,对一维情况即为其子域的长度; $\left| J \right|=\left| {{d}g^{-1}\left( x\right)/{d}x} \right|$为Jacobi行列式, 对一维情况即为其导数之绝对值,$g^{-1}(\cdot )$表示$g(\cdot )$的反函数.
(2) 函数$g(\cdot )$为非单调函数. 此时, 总可以将其定义域剖分为若干个不相交子域, 并使之在每个子域中为单调函数[21]. 记其相应的反函数为$g_{j}^{-1} (\cdot )(j=1,2,\cdots ,N)$, $N$为子域划分段数, 此时有
式中$\varOmega_{\varTheta ,j} $为$\varOmega_{\varTheta } $之第$j$个子域,且每个这样的子域都满足$\varOmega_{\varTheta ,j} =g_{j}^{-1} \left( {\varOmega_{X} }\right)$; $\left| {J_{j} } \right|=\left| {{d}g_{j}^{-1} \left( x\right)/{d}x} \right|$为第$j$个反函数对应的Jacobi行列式.
若$m=1$而$n>1$, 则式(1)退化为$X(\varpi)$ $=$ $g\left( {{\boldsymbol \varTheta}(\varpi)} \right)$, 而式(3)相应地成为$\int_{\varOmega_{X} }{p_{X} \left( x \right){d}x} =\int_{\varOmega_{\boldsymbol \varTheta } } {p_{\boldsymbol \varTheta }\left( \boldsymbol \theta \right){d}\boldsymbol \theta } $. 由此, 有
对其进一步化简可得
其中${\boldsymbol \theta }_{\sim 1} =\left( {\theta_{2} ,\theta_{3} ,\cdots,\theta_{n} } \right)^{T}$, $g_{j}^{-1} \left( {x,{\boldsymbol \theta }_{\sim 1} } \right)$是以${\boldsymbol \theta }_{\sim 1} $为参数时$g(\cdot )$的第$j$个反函数, $\left| {J_{j} } \right|=\left| {\partial g_{j}^{-1} \left( {x,{\boldsymbol \theta }_{\sim 1} } \right)/\partial x}\right|$.
另一方面, 也可以通过Dirac函数$\delta _{{D}} \left[ \cdot \right]$的筛选性质体现的映射关系获得输出变量的概率密度函数. 事实上, 对于上述情形, 有
其中$p_{X\vert {\boldsymbol \varTheta }} \left( {x\vert {\boldsymbol \varTheta}={\boldsymbol \theta }} \right)$为条件概率密度函数. 由于存在映射关系$X(\varpi)=g\left( {{\boldsymbol \varTheta }(\varpi)} \right)$,因而有$p_{X\vert {\boldsymbol \varTheta }} \left( {x\vert {\boldsymbol \varTheta}={\boldsymbol \theta }} \right)=\delta _{{D}} \left[ {x-g\left({{\boldsymbol \theta }} \right)} \right]$, 由此可得式(8)中的第二个等式.
显然, 可认为式(8)是输出概率密度函数$p_{X} \left( x\right)$关于输入概率密度函数$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta }}\right)$的隐式表达,而式(4)、式(5)和式(7)则是对式(8)进行积分获得的显式解析表达. 上述各式表明,输出概率密度函数$p_{X} \left( x \right)$由输入概率密度函数$p_{{\boldsymbol \varTheta}} \left( {{\boldsymbol \theta }} \right)$所决定,因而可以进一步写成如下算子形式的表达
其中下标$g$用来表示算子$\psi_{g} $决定于函数$g(\cdot )$.在实际工程中, 函数$g(\cdot )$描述了问题的物理机制, 因此,式(9)表明: 不确定性在系统中的传播, 即从输入概率信息到输出概率信息的变换,是由系统物理机制所决定的. 从函数空间来看,这一变换的实现是通过由物理机制决定的算子作用来实现的.
对于$m\geqslant 2$的情况, 上述基本思想是完全类似的. 限于篇幅, 兹不赘述.
1.2 随机动力系统中不确定性传播的泛函空间观点
这一泛函空间观点, 对于随机动力系统中的不确定性传播也是适用的. 为了明确起见,考虑如下随机动力系统[21]其中${X}=\left( {X_{1} ,X_{2} ,\cdots ,X_{m} } \right)^{T}$为状态向量, ${\dot{{X}}}=\partial {X}/\partial t$;${G}=\left( {G_{1} ,G_{2} ,\cdots ,G_{m} } \right)^{T}$为确定性函数向量; ${\boldsymbol \varTheta }=\left( {\varTheta_{1} ,\varTheta_{2},\cdots ,\varTheta_{n} } \right)^{T}$为刻画随机源信息的随机向量,其概率密度函数$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta }} \right)$已知,这里${\boldsymbol \theta }=\left( {\theta_{1} ,\theta_{2} ,\cdots ,\theta_{n}} \right)^{T}$为其对应的可实现值. 式(10)的初始条件为${X}\left( 0\right)={x}_{0} $, 其中${x}_{0} $为确定性初始值.
根据概率守恒原理[20], 记$\left( {{X}\left( t\right),{\boldsymbol \varTheta }} \right)$的联合概率密度函数为$p_{{X\boldsymbol \varTheta }}\left( {{x},{\boldsymbol \theta },t} \right)$,则其满足如下广义概率密度演化方程[21]
其初始条件为
求解上述方程即可得${X}\left( t \right)$的概率密度函数
由此可见, 若已知初始随机源概率信息(输入概率信息) $p_{{\boldsymbol \varTheta }}$ $\left({{\boldsymbol \theta }} \right)$,即可通过求解式(11)、式(12)和式(13)获得输出概率信息$p_{{X}} \left({{x},t} \right)$. 因此,上述3个方程决定了从输入概率密度函数到输出概率密度函数的变换.从函数空间来看, 即存在如下算子方程
这里用下标${G}$表示算子$\psi_{{G}} $由方程(10)所决定,而方程(10)是描述系统演化内在物理机制的方程. 因此, 从这一角度来看,概率密度演化理论[20-21]正是算子方程(14)的具体实现过程,而不确定性在系统中从输入概率信息到输出概率信息的传播,是由系统物理机制(通过确立传播算子)来决定的. 从这一意义上说, $\psi_{{G}} $可称为不确定性传播算子.
事实上, 式(9)和式(14)中的算子$\psi_{g} $和$\psi_{{G}} $,正是Frobenius-Perron算子[21-22]. 值得注意的是, 上述算子是线性算子,仅依赖于物理方程(例如以上的$g(\cdot )$和${G}(\cdot ))$, 而与输入概率信息(例如以上$p_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta }} \right))$的形式是无关的. 换言之,物理系统本身不因随机源概率信息的变化而变化. 这一随机系统客观不变性原理,正是上述不确定性传播算子具有不变性的内在原因[23].
2 基于Fréchet导数的整体灵敏度
2.1 基于Fréchet导数的整体灵敏度指标定义
基于上述泛函空间理解, Chen等[16]提出了基于Fréchet导数的整体灵敏度指标. 为清晰起见,下面以随机系统(10)与相应的不确定性传播算子方程(14)为基础简要阐述其基本思想.考虑如下问题: 若输入概率密度函数发生了变化, 例如从$p_{{\boldsymbol \varTheta }}\left( {{\boldsymbol \theta }} \right)$变为$\tilde{{p}}_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta }} \right)$,则系统响应(输出)的概率密度函数也将相应地发生改变, 且由算子方程(14)可知
为方便起见, 记输入概率密度函数的变化为$\delta p_{\boldsymbol \varTheta}=\tilde{p}_{\boldsymbol \varTheta}(\boldsymbol \theta)$ $-$ $p_{\boldsymbol \varTheta}(\boldsymbol \theta)$, 输出概率密度函数的变化为$\delta p_{X} =\tilde{p}_{X}(x,t)$ $-$ $p_{X}(x,t)$. 由式(15)两端分别减去式(14)可得
由此可见, 若输入概率信息具有变化,则输出概率信息的变化亦由不确定性传播之线性算子所决定.
事实上, 上述概率信息“变化”的传播算子, 就是不确定性传播算子的Fréchet导数, 其严格的数学定义为[24]:
令$V$和$W$为Banach空间, 定义映射$\psi :V\to W$.若存在一有界线性算子${\cal F}_{\psi } \in {\cal L}\left( {V,W} \right)$,满足
则称${\cal F}_{\psi } $为在$u_{0} \in U\subset V$处的算子$\psi$的Fréchet导数, 记为${\cal F}_{\psi } ={\psi }'\left[ {u_{0} } \right]$.式中, $\left\| \cdot \right\|_{V} $和$\left\| \cdot \right\|_{W}$分别对应Banach空间$V$和$W$上的范数.
在本文中, Banach空间$V$上的范数采用基于总变差距离的范数
其中, ${\cal L}^{1}\left( {\varOmega_{{\boldsymbol \varTheta }} }\right)$为定义在$\varOmega_{{\boldsymbol \varTheta }} $上的$\ell_{1} $范数, 即
而空间$W$上的范数同样可定义为
由于不确定性传播算子是线性算子,因此概率信息“变化”的传播算子就是不确定性传播算子本身.
在实际计算中, 可以利用Fréchet导数与Gateaux导数的等价性.Gateaux导数的定义为: 令$V$和$W$为Banach空间, 定义映射$\psi :V\to W$. 若存在一有界线性算子${\cal G}_{\psi } \in \left( {V,W} \right)$, 满足
则称${\cal G}_{\psi } $为在$u_{0} \in U\subset V$处的算子$\psi$的Gateaux导数, 记为${\cal G}_{\psi } ={\psi }'\left[ {u_{0} } \right]$.
显然, Fréchet导数即为Gateaux导数. 反之, 若Gateaux导数在点$u_{0}$处连续或式(21)对$h$一致成立且$\left\| h \right\|_{V} =1$, 则点$u_{0}$处的Gateaux导数亦为Fréchet导数[24].
无论是确定性系统还是随机系统,系统灵敏度的定义都可以阐述为[25]“系统灵敏度是输入量的扰动导致输出量的扰动这一效应的度量”.基于这一观点,在经典的确定性系统分析中往往采用函数的偏导数作为系统灵敏度的度量. 显然,在随机系统中采用算子的Fréchet导数作为系统的灵敏度度量,乃是一种自然的推广. 但是, 以函数的偏导数作为确定性系统的灵敏度,仅能反映系统在给定输入值一点处的变化情况, 因而是一种局部的度量;而对随机系统而言,Fréchet导数给出了在系统输入全部可能值(通过概率密度函数刻画)处的变化导致在系统输出全部可能值上的变化情况,因而其本质是一种整体灵敏度指标(global sensitivity index, GSI)[16].
基于Fréchet导数的整体灵敏度指标有4个重要性质:
(1) 无论物理系统本身是线性还是非线性的, 因为Fréchet导数是线性算子,因此基于Fréchet导数的整体灵敏度指标与$p_{{\boldsymbol \varTheta }}$的具体形式是无关的(类似于线性函数$x=k\theta$中的局部灵敏度${d}x/{d}\theta =k$是与$\theta $取值无关的).这意味着, 这一整体灵敏度指标不仅适用于输入概率密度信息具有小扰动的情况,也适用于输入概率密度信息具有大扰动、甚至若干输入变量从随机变量变为确定性变量的情况(例如概率密度函数坍塌为Dirac函数).
(2) 基于性质(1), 若考虑输入概率模型具有某种程度的认知不确定性[1-2,17-18], 例如, 输入概率密度函数可能为$p_{{\boldsymbol \varTheta }}^{\left( 1\right)} $或$p_{{\boldsymbol \varTheta }}^{\left( 2 \right)} $,则这种认知不确定性传播的量化可以直接由该整体灵敏度指标给出, 即${\cal F}_{\psi } \circ \left( {p_{{\boldsymbol \varTheta }}^{\left( 2 \right)}-p_{{\boldsymbol \varTheta }}^{\left( 1 \right)} } \right)$.
(3) 对于随机动力系统, 系统响应是时变的, 因此基于Fréchet导数的整体灵敏度指标为一时变量. 在此情况下,可根据问题的需要考虑系统响应的某代表性值, 如系统响应时程的极大值, 即${Z}_{{ext}} =\max \left\{ {{X}\left( t\right),\mbox{, }t\in \left[ {0,T} \right]} \right\}$, 简记为${Z}$.
(4) 注意到Fréchet导数是线性有界算子, 在式(17)、式(18)和式(20)定义下的范数为1(证明见第2.3节).由于Fréchet导数就是不确定性传播算子本身, 因此, 不确定性传播算子的范数为1. 这在本质上是概率守恒原理[20]在泛函空间观点下的数学表述之一.
2.2 整体灵敏度指标的参数化表达
在工程中, 特别具有实际意义的是参数化表达的整体灵敏度指标. 对绝大多数实际工程问题而言, 系统输入参数${\boldsymbol \varTheta}$一般为连续型随机向量, 记其概率密度函数为$p_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$, 其中${\boldsymbol \zeta }=\left({\zeta_{1} ,\zeta_{2} ,\cdots ,\zeta_{s} } \right)^{T}$为其对应的分布参数向量. 此时, $p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$的变分为在文献[16]中, 考虑各分布参数之间相互独立变化, 即当仅考虑$\delta \zeta_{j}\ne 0$而其他为$0$时, 给出了基于Fréchet导数的整体灵敏度指标的参数化表达
引入式(14), 式(23)可以进一步改写为
这表明基于Fréchet导数的整体灵敏度指标参数化表达是不确定性传播算子作用于单位化的输入函数变化上的结果.因此, 上述整体灵敏度指标参数化表达的几何意义为单位化的输入函数变化作为抽象空间中的点(向量)经系统作用后的“方向”及“长度”变化.
更具体地, 在实际工程问题中, 有时通过随机变量$\varTheta$的统计矩参数来近似构造其概率密度函数. 例如, 四阶矩正态变换方法假定随机变量$\varTheta $可表达为标准正态随机变量$U$的多项式函数[26]
其中, 参数$a_{0} $, $a_{1} $, $a_{2} $和$a_{3} $由$\varTheta$的前四阶矩即其均值$\mu_{\varTheta } $、标准差$\sigma_{\varTheta }$、偏度$\alpha_{\varTheta } $和峰度$\beta_{\varTheta } $计算得到
由随机变量的函数变换法则[19-21], 随机变量$\varTheta $的概率密度函数可表示为
其中, $\phi \left( u \right)$为标准正态分布的概率密度函数,反函数$u=u^{-1}\left( \theta \right)$的显式表达详见文献[27].
可见, 此时$\varTheta $的概率密度函数由其前四阶矩确定. 不妨考虑随机向量${\boldsymbol \varTheta }=\left( {\varTheta_{1} ,\varTheta_{2} ,\cdots ,\varTheta_{n} } \right)^{T}$相互独立, 而对于非独立情况可通过Nataf变换[28]、Copula函数[29]或随机函数模型[3]转换为独立情况.此时, 该随机向量的概率密度函数可表示为
式中, $p_{\Theta_{\ell}}$为第$\ell $个随机变量$\Theta_{\ell}$的概率密度函数, $\left(\mu_{\Theta}^{(\ell)}, \sigma_{\Theta}^{(\ell)}, \alpha_{\Theta}^{(\ell)}, \beta_{\Theta}^{(\ell)}\right)$分别对应$\Theta_{\ell}$的均值、标准差、偏度和峰度.由式(23)和式(28)可得随机系统(1)基于Fréchet导数的整体灵敏度指标为
矩表达形式中前四阶矩的物理意义是十分明确的,其对应的全局灵敏度指标反映了各个输入变量的前四阶矩信息变化对输出概率密度函数的影响.
更一般地, 若${\boldsymbol \varTheta }$的概率密度函数是已知的但是非参数化的,如$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta }}\right)$是通过核密度估计[1]或概率密度演化理论[21]得到的数值解,则可通过$\varepsilon $-等效分布的方法将其参数化. 为此,记随机变量${\boldsymbol \varTheta }$的概率密度函数为$p_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta }} \right)$. 若${\boldsymbol \varTheta }$可以由有限个参数进行表达,则称$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta}} \right)$为$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta }}\right)$的$\varepsilon $-等效分布, 满足
其中, ${\boldsymbol \zeta }$为有限维参数向量, 可由优化算法[12]求解式(30)给出; $\varepsilon $为给定的容许误差, 如0.01.注意, 这里$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$可以是任意形式的. 例如,有确定形式的概率密度函数[2-3]、由矩表达形式给出的概率密度函数[26-27]或由级数展开方法[30]计算得到的概率密度函数等. 在上述基础上将$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$替换为$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$, 则随机系统(1)基于Fréchet导数的整体灵敏度指标的参数化表达可统一为式(23).
2.3 重要性测度
如式(23)所示的整体灵敏度指标描述了系统响应${X}$的值域范围内每一点${x}$处概率密度的灵敏度. 进一步, 对整体灵敏度指标取$\left\| \cdot \right\|_{W} $范数可得其重要性测度现分别考察式(31)的几何意义和物理意义.
2.3.1 几何意义
考虑Fréchet导数的算子范数, 其定义为[24]
式中各符号同式(17)中定义. 任意取${\boldsymbol \varTheta }$的概率密度函数为$p_{2} $和$p_{1} $, 令$p=p_{2} -p_{1} $, 由均值定理有[24]
此处利用了不确定性传播算子是线性算子的性质.
注意$tp_{2} +\left( {1-t} \right)p_{1} $为概率密度函数. 因此, 由概率相容性和不确定性传播算子方程(9)有
整理式(32),$\sim \!$式(34)即可得$\left\| {{\cal F}_{\psi } } \right\|_{\ell _{1} } =1$. 显然, 此即由式(31)所定义的重要性测度的上限
因此, 重要性测度总是小于1. 结合前述式(24)的几何意义表明, 重要性测度的几何意义是单位化(即长度为1)的输入函数变化作为抽象空间中的点(向量)经系统作用后的长度,这一长度不可能超过1, 且该长度越大, 则其重要性越大.
2.3.2 物理意义
考虑式(31)的差分表达形式. 即, 取${\boldsymbol \zeta }$沿$\zeta_{j} $方向的增量为$\Delta \zeta_{j} $(一般取$\Delta \zeta_{j} =10^{-2}\zeta _{j} $或$10^{-3}\zeta_{j} $; 若$\zeta_{j} =0$则取$\Delta \zeta_{j} =10^{-2}$或$10^{-3})$, 并简记$p_{{X}} \left( {{x};{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} } \right)=p_{{X}}^{+} $, $p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} } \right)=p_{{\boldsymbol \varTheta }}^{+} $, $p_{{X}} \left( {{x};{\boldsymbol \zeta }} \right)=p_{{X}} $和$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)=p_{{\boldsymbol \varTheta }} $, 式(31)等价于
其中, ${e}_{j} $为$s$维单位向量, 且仅其第$j$位上元素为$1$而其它为$0$;$d_{{TV}} \left( {p,q} \right)$为概率密度函数$p$和$q$的总变差距离(totalvariation distance), 即[16]
类比于连续介质力学中的变形梯度并观察式(36)不难发现,由式(31)所定义的整体灵敏度对应的重要性测度反映了某种“体积变化率”(如图1所示):当第$j$个分布参数$\zeta_{j} $发生扰动,使得输入概率空间产生“体积变化”$d_{{TV}} \left( {p_{{\boldsymbol \varTheta}}^{+} ,p_{{\boldsymbol \varTheta }} } \right)$后,随机系统对应的输出概率空间亦相应地产生了“体积变化” , 即$d_{{TV}}\left( {p_{{X}}^{+} ,p_{{X}} } \right)$. 由式(36)定义可知,这两者的比值的物理意义便是由不确定性传播导致的“体积变化率”. 需要注意的是,由式(36)定义的“体积变化率”恒为正值,因此可视其为“体积变化率”的绝对值或将其称为“幅值”. 如此,由式(31)给出的整体灵敏度指标对应的重要性测度的物理意义是十分明确的.
图1
新窗口打开|下载原图ZIP|生成PPT图1基于Fréchet导数的整体灵敏度指标对应的重要性测度的物理意义
Fig. 1Physical meanings of the importance measure with respect to the Fréchet derivative based GSI
2.4 方向灵敏度
另一方面, 考虑如(23)所示的整体灵敏度指标在系统响应${X}$值域范围内的方向灵敏度式中, sgn$\left( A \right)$为符号函数: 当$A\geqslant 0$时取$+1$(代表“正”或“$\mbox{+}$”)、否则取$-1$ (代表“负”或“$-$”). 计算整体灵敏度指标在$\varOmega_{{X}} $上的积分, 易得
式(39)是由概率相容条件导致的必然结果. 定义$\varOmega_{{X},j}^{+}:{x}\in \varOmega_{{X}} $满足${\cal F}_{\psi ,j} \geqslant 0$和$\varOmega_{{X},j}^{-} :\mbox{, }\varOmega_{{X}} \backslash \varOmega_{{X},j}^{+} $, 对式(39)进行正负分解后有
事实上, 由式(31)可知
利用式(38)的定义, $\varOmega_{{X},j}^{+} $和$\varOmega_{{X},j}^{-}$的统一表达为
其中$I\left\{ A \right\}$为示性函数, 当$A$为真时取1, 否则取0.
由式(38)所定义的方向灵敏度的意义在于, 系统响应的值域$\varOmega_{{X}}$可以通过式(42)划分为正、负两个子值域(可能是多连通域). 具体而言: (1)输入分布参数的微小正增长, 如$\Delta \zeta_{j} >0$将导致正值域$\varOmega_{{X},j}^{+} $内的概率密度函数增大、负值域$\varOmega_{{X},j}^{-}$内的概率密度函数减小; 反之, 输入分布函数的微小负增长将导致正值域$\varOmega_{{X},j}^{+} $内概率密度函数减小, 而负值域$\varOmega_{{X},j}^{-}$内的概率密度函数增大. 因此, $\varOmega_{{X},j}^{-} $与$\varOmega_{{X},j}^{+} $分别表示了由于参数变化导致的概率净减和概率净增值域; (2)若按式(38)定义可得有限个正、负值域,则正或负值域的不连续子值域个数(如图2)可反映物理系统的复杂程度. 一般而言,系统非线性程度越复杂, 其不交叉子值域个数越多.在第4节将结合具体实例对其进行详细讨论.
图2
新窗口打开|下载原图ZIP|生成PPT图2基于Fréchet导数的整体灵敏度指标对应的方向灵敏度示意图
Fig. 2Schematic representation of direction sensitivity with respect to the Fréchet derivative based GSI
2.5 与矩独立灵敏度指标的联系
考虑随机系统(1), 基于矩独立的全局灵敏度指标定义为[10]其中, $\mathbb{E}_{\varTheta_{i} } \left[ \cdot \right]$为关于$\varTheta_{i}$的期望算子; $p_{{X}_{?i} } \left( {{x}\vert {\boldsymbol \varTheta}_{i} =\bar{{\theta }}_{i} } \right)$为条件概率密度函数.
式(43)可进一步写作
在式(37)定义下, $s_{i} \left( {\varTheta_{i} =\bar{{\theta }}_{i} } \right)$即为概率密度函数$p_{{X}_{i} } $与$p_{{X}} $的总变差距离, 不妨记为$d_{{TV}} \left( {p_{{X}_{i} } ,p_{{X}} } \right)$. 考虑各输入随机变量$\left( {\varTheta_{1} ,\varTheta _{2} ,\cdots ,\varTheta_{n} } \right)^{T}$是相互独立的, 则$p_{{X}_{i} } \left( {{x}\vert {\varTheta }_{i} =\bar{{\theta }}_{i} } \right)$对应的输入概率密度函数有
而此时$p_{{\boldsymbol \varTheta }_{i} } $与$p_{{\boldsymbol \varTheta }} $的总变差距离为
上式意味着当第$i$个随机变量$\varTheta_{i} $坍塌为确定性变量$\bar{{\theta}}_{i} $时, 输入概率空间产生了单位体积的变化. 此时,基于矩独立的全局灵敏度指标等价于
在第2.3节中, 式(36)解释了基于Fréchet导数的整体灵敏度指标的重要性测度的物理意义:它是随机系统在点$p_{{\boldsymbol \varTheta }} $处沿其分布参数$\zeta_{j}$方向的(绝对)体积变化率. 类比于此, 式(47)可理解为第$i$个随机变量$\varTheta_{i}$沿其所有可实现值$\bar{{\theta }}_{i} $方向(绝对)体积变化的期望值.
事实上, 若考虑参数化形式下的${\boldsymbol \varTheta}$的概率密度函数$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta}} \right)$, 其中${\boldsymbol \zeta }=\left( {{\boldsymbol \zeta }_{1} ,{\boldsymbol \zeta }_{2} ,\cdots ,{\boldsymbol \zeta }_{n} } \right)^{T}$分别对应$\left( {\varTheta _{1} ,\varTheta_{2} ,\cdots ,\varTheta_{n} } \right)^{T}$的分布参数向量.当第$i$个随机变量$\varTheta_{i} $坍塌为确定性变量$\bar{{\theta }}_{i} $时, 一方面, 其对应的概率密度函数由$p_{\varTheta_{i} } \left( {\theta_{i} ;\zeta _{i} } \right)$变为$\delta _{{D}} \left( {\theta_{i} -\bar{{\theta }}_{i} } \right)$; 另一方面, 若有${\boldsymbol \zeta }_{i} =\left( {\mu_{i},\sigma_{i} } \right)^{T}$, 其中$\mu_{i} $为均值、$\sigma_{i}$为标准差, 则有$\delta _{{D}} \left( {\theta_{i} -\bar{{\theta }}_{i}} \right)=p_{\varTheta_{i} } \left( {\theta_{i} ;\mbox{, }\mu_{i} \to\bar{{\theta }}_{i} ,\sigma_{i} \to 0} \right)$. 一个具体的例子是, 若$\varTheta_{i} $服从均值为$\mu_{i} $、标准差为$\sigma_{i} $的正态分布, 则当$\varTheta_{i} $坍塌为确定性变量$\bar{{\theta }}_{i} $时, 其概率密度函数是$\exp \left\{ {-\left( {\theta_{i} -\mu_{i} } \right)^{2}/2\sigma_{i}^{2} } \right\}/\sqrt {2\pi } \sigma_{i} $在$\mu _{i} \to \bar{{\theta }}_{i} $和$\sigma_{i} \to 0$时的极限形式.
在此理解的基础上, 记第$i$个随机变量$\varTheta_{i} $向其均值$\mu_{i}$坍塌前${\boldsymbol \varTheta }$的联合概率密度函数为$p_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta };\mu_{i} ,\sigma_{i} }\right)$、坍塌后为$p_{{\boldsymbol \varTheta }_{i} } \left( {{\boldsymbol \theta };\mu_{i} ,0} \right)$, 则由均值定理有[24]
这里${\cal F}_{\psi ,i} $为基于Fréchet导数的整体灵敏度指标的参数化表达. 注意到${\cal F}_{\psi ,i} $为一线性算子, 式(48)可进一步写为
式中不等式右边两项分别对应点$p_{{\boldsymbol \varTheta }_{i} } $和$p_{{\boldsymbol \varTheta}} $处的关于标准差参数的重要性测度. 在点$p_{{\boldsymbol \varTheta }} $处有
其中${\cal S}_{i} \left( {\sigma_{i} } \right)$为由式(31)给出的第$i$个随机变量$\varTheta_{i} $关于其标准差参数$\sigma_{i} $基于Fréchet导数的整体灵敏度指标的重要性测度.
3 数值方法
基于Fréchet导数的整体灵敏度计算分为两部分: (1) 系统响应的概率密度函数关于输入分布参数的导数; (2) 输入概率密度函数关于分布参数的导数的范数. 对于第(2)部分, 由于$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$是已知的, 因此计算该范数项是容易的; 而对于第(1)部分则可以通过中心差分法对其进行数值近似, 即式中, ${e}_{j} $和$\Delta \zeta_{j} $意义同式(36).
对$j=1,2,\cdots ,s$, 式(51)包含系统响应的概率密度函数$p_{{X}} \left({{x};{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} }\right)$和$p_{{X}} \left( {{x};{\boldsymbol \zeta }-{e}_{j}\Delta \zeta_{j} } \right)$,它们分别对应输入概率密度函数$\tilde{{p}}_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta };{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} }\right)$和$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta};{\boldsymbol \zeta }-{e}_{j} \Delta \zeta_{j} } \right)$.
对一些简单随机系统而言, 可通过显式解析表达的Frobenius-Perron算子[31]直接求解式(14), 从而获得系统响应的概率密度函数. 然而, 对绝大多数实际工程问题而言, Frobenius-Perron算子难以直接进行理论求解. 因此, 本文将基于概率密度演化理论(probability density evolution method, PDEM)[21]进行求解, 从而获得$p_{{X}} \left( {{x}} \right)$, 详见第3.1节.
一般而言, 对$j=1,2,\cdots,s$分别进行概率密度演化分析即可获得基于Fréchet导数的整体灵敏度指标. 这总计需要$2Ns$次确定性分析. 具体而言, 不妨考虑${\boldsymbol \varTheta }$为$10$维独立的随机向量、且各随机变量均为两参数分布. 一般取概率密度演化分析中的概率空间剖分数为$500$, 则总共需要求解$20000$次确定性物理系统！这毫无疑问地给实际工程问题分析带来了极大的挑战. 最近, Chen和Wan[17]从概率测度变换(change of probability measure, COM)的角度提出了降低确定性分析次数的方法, 详见第3.2节.
3.1 随机系统不确定性传播的概率密度演化理论
基于概率密度演化理论, 对系统(10)的数值求解算法的基本步骤包括[21]:(1) 概率空间剖分与点集优选: 记${\boldsymbol \varTheta }$的样本空间为$\varOmega_{{\boldsymbol \varTheta }} $, 将其剖分为$N$个子域$\varOmega_{1} ,\varOmega_{2} ,\cdots,\varOmega_{N} $, 并使其满足$\cup_{q=1}^{N} \varOmega_{q} =\varOmega_{{\boldsymbol \varTheta }} $且对任意$1\leqslant p\ne q\leqslant N$有$\varOmega_{p} \cap \varOmega_{q} =\emptyset $. 对$q=1,2,\cdots ,N$, 在每个子域中选取一个代表点,记为${\boldsymbol \theta }_{q} \in \varOmega_{q} $, 并计算其对应的赋得概率.赋得概率定义为子域$\varOmega_{q} $相对于全域$\varOmega_{{\boldsymbol \varTheta }}$的Voronoi体积, 即
此步骤所涉及的点集优选策略及其点集误差分析可进一步见文献[32,33].
(2) 对$q=1,2,\cdots ,N$, 给定代表点${\boldsymbol \theta }={\boldsymbol \theta }_{q} $并求解系统(10)得到${\dot{{X}}}\left( {{\boldsymbol \theta }_{q} ,t} \right)={G}\left( {{X}\left( {{\boldsymbol \theta }_{q} ,t} \right),{\boldsymbol \theta }_{q} ,t} \right)$.
(3) 对$q=1,2,\cdots ,N$, 求解半离散化的广义概率密度演化方程
其初始条件为$p_{{X\boldsymbol \varTheta }}^{\left( q \right)} \left({{x},{\boldsymbol \theta }_{q} ,0} \right)=\delta _{{D}} \left[{{x}-{x}_{0} } \right]P_{q} $、边界条件为$p_{{X\boldsymbol \varTheta}}^{\left( q \right)} \left( {+\infty ,{\boldsymbol \theta }_{q} ,t}\right)=p_{{X\boldsymbol \varTheta }}^{\left( q \right)} \left( {-\infty,{\boldsymbol \theta }_{q} ,t} \right)=0$.
(4) 对上述结果进行求和集成
需要说明的是, 若考虑系统(10)的极值分布$p_{{Z}} \left( {{z}} \right)$, 其中${Z}=\max \left\{ {{X}\left( t \right),\mbox{, }t\in \left[ {0,T} \right]} \right\}$. 此情况下, 基于等价极值事件的基本思想、通过构造虚拟随机过程即可获得该极值分布$p_{{Z}} \left( {{z}} \right)$, 详见文献[34].
3.2 随机系统不确定性传播的概率测度变换
考虑两组不同的输入概率密度函数$p_{{\boldsymbol \varTheta }}^{\left( 1 \right)} \left( {{\boldsymbol \theta }} \right)$和$p_{{\boldsymbol \varTheta }}^{\left( 2 \right)} \left( {{\boldsymbol \theta }} \right)$, 分别对应系统响应的概率密度函数$p_{{X}}^{\left( 1 \right)} \left({{x},t} \right)$和$p_{{X}}^{\left( 2 \right)} \left( {{x},t} \right)$. 其中, $p_{{X}}^{\left( 1 \right)} \left( {{x},t} \right)$可由第3.1节中的概率密度演化理论给出, 而$p_{{X}}^{\left( 2 \right)} \left( {{x},t}\right)$则由概率测度变换(COM)得到, 基本步骤如下[17]:(1) 由概率密度演化理论计算$p_{{\boldsymbol \varTheta }}^{\left( 1 \right)} \left( {{\boldsymbol \theta }} \right)$对应的系统响应概率密度函数$p_{{X}}^{\left( 1 \right)} \left( {{x},t} \right)$, 保存分析中的代表点集${\cal M}^{\left( 1 \right)}=\left\{ {{\boldsymbol \theta }_{q}^{\left( 1 \right)} ,P_{q}^{\left( 1 \right)} } \right\}_{q=1}^{N} $及其对应的确定性解$\left\{ {{\dot{{X}}}\left( {{\boldsymbol \theta }_{q}^{\left( 1 \right)} ,t} \right)} \right\}_{q=1}^{N} $.
(2) 当输入概率信息由$p_{{\boldsymbol \varTheta }}^{\left( 1 \right)} \left( {{\boldsymbol \theta }} \right)$变为$p_{{\boldsymbol \varTheta }}^{\left( 2 \right)} \left( {{\boldsymbol \theta }} \right)$, 保持代表性点不变、重新计算赋得概率有
其中, 子域$\varOmega_{q} $由$\left\{ {{\boldsymbol \theta }_{q}^{\left( 1 \right)}} \right\}_{q=1}^{N} $确定.
(3) 重新计算第3.1节中第3步, 将初始条件更换为$p_{{X\varTheta }}^{\left({q,2} \right)} \left( {{x},{\boldsymbol \theta }_{q}^{\left( 1 \right)},0} \right)=\delta _{{D}} \left[ {{x}-{x}_{0} }\right]P_{q}^{\left( 2 \right)} $.
(4) 对计算结果进行集成, 即
$\begin{eqnarray}\label{eq56}p_{{X}}^{\left( 2 \right)} \left( {{x},t}\right)=\sum\limits_{q=1}^N {p_{{X\varTheta }}^{\left( {q,2} \right)}\left( {{x},{\boldsymbol \theta }_{q}^{\left( 1 \right)} ,t} \right)}\end{eqnarray}$
3.3 基于Fréchet导数的整体灵敏度指标的概率密度演化概率测度变换算法
基于Fréchet导数的整体灵敏度指标的数值计算流程如下:(1) 由式(30)给出输入概率密度函数$p_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta}} \right)$的$\varepsilon $-等效分布$\tilde{{p}}_{{\boldsymbol \varTheta }} \left({{\boldsymbol \theta };{\boldsymbol \zeta }} \right)$, 其中${\boldsymbol \zeta }=\left({\zeta_{1} ,\zeta_{2} ,\cdots ,\zeta_{s} } \right)^{T}$为其对应的分布参数向量.
(2) 由概率密度演化理论计算输入概率密度函数$\tilde{{p}}_{{\boldsymbol \varTheta }}\left( {{\boldsymbol \theta };{\boldsymbol \zeta }}\right)$对应的系统输出概率密度函数$p_{{X}} \left( {{x}}\right)$.
(3) 对$j=1,2,\cdots ,s$,由第3.2节概率测度变换给出输入概率密度函数$\tilde{{p}}_{{\boldsymbol \varTheta }}\left( {{\boldsymbol \theta };{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} }\right)$对应的系统输出概率密度函数$p_{{X}} \left({{x};{\boldsymbol \zeta }+{e}_{j} \Delta \zeta_{j} } \right)$;以及输入概率密度函数$\tilde{{p}}_{{\boldsymbol \varTheta }} \left( {{\boldsymbol \theta};{\boldsymbol \zeta }-{e}_{j} \Delta \zeta_{j} }\right)$对应的系统输出概率密度函数$p_{{X}} \left({{x};{\boldsymbol \zeta }-{e}_{j} \Delta \zeta_{j} } \right)$.
注意这一步不需要额外的确定性物理系统分析.
(4) 对$j=1,2,\cdots ,s$计算各分布参数对应的范数$\left\|\partial \tilde{p}_{\Theta}(\theta ; \zeta) / \partial \zeta_{j}\right\|_{V}$.
(5) 计算式(51).
4 算例分析
4.1 算例一: 正态随机变量的线性组合函数
考虑一随机系统式中, $\varTheta_{1} $、$\varTheta_{2} $和$\varTheta_{3}$为相互独立的正态随机变量, 其均值为$\mu_{1} $, $\mu_{2} $和$\mu_{3} $,标准差为$\sigma_{1} $, $\sigma_{2} $和$\sigma_{3} $. 易知,$X$服从均值为$\mu_{1} -\mu_{2} +\sqrt 2 \mu_{3} $,标准差为$\sqrt {\sigma_{1}^{2} +\sigma_{2}^{2} +2\sigma_{3}^{2} } $的正态分布.式(57)所示随机系统用于阐述基于Fréchet导数的整体灵敏度指标的性质,并检验本文所采用的PDEM-COM方法的计算精度[16]. 本例中取$\mu_{i}=0$和$\sigma_{i} =1$ $(i=1,2,3)$.
4.1.1 基于Fréchet导数的整体灵敏度指标
记相对于分布参数$\mu_{1} $, $\mu_{2} $和$\mu_{3}$的整体灵敏度指标分别为${\cal F}\left( {\mu_{1} } \right)$, ${\cal F}\left({\mu_{2} } \right)$和${\cal F}\left( {\mu_{3} } \right)$;相对于分布参数$\sigma_{1} $、$\sigma_{2} $和$\sigma_{3}$的整体灵敏度指标分别为${\cal F}\left( {\sigma_{1} } \right)$, ${\cal F}\left( {\sigma_{2} } \right)$和${\cal F}\left( {\sigma_{3} } \right)$.本例中已知$\varTheta_{1} $, $\varTheta_{2} $和$\varTheta_{3} $为正态随机变量,因此其$\varepsilon $-等效分布即为正态分布.基于Fréchet导数的整体灵敏度指标绘于图3和图4之中,对应的重要性测度如表1所示, 按其大小排序有: $\mu_{3} >\mu_{2} =\mu_{1}$以及$\sigma_{3} >\sigma_{2} =\sigma_{1} $. 从表中可见,所有的重要性测度均满足式(35), 即其值不超过1.
图3
新窗口打开|下载原图ZIP|生成PPT图3相对于均值参数的基于Fréchet导数的整体灵敏度指标(算例一)
Fig. 3Fréchet derivative based GSI in terms of the mean values (Example 1)
图4
新窗口打开|下载原图ZIP|生成PPT图4相对于标准差参数的基于Fréchet导数的整体灵敏度指标(算例一)
Fig. 4Fréchet derivative based GSI in terms of the standard deviations (Example 1)
Table 1
表1
表1基于Fr′echet导数的整体灵敏度指标的重要性测度(算例一)
Table 1
 |
新窗口打开|下载CSV
在式(57)所示的随机系统中,各输入随机变量关于系统输出是线性且独立同分布的关系. 因此,线性系数反映了各个输入随机变量与系统响应的物理关系. 具体而言,分别就均值参数和标准差参数考察$\varTheta_{1} $, $\varTheta_{2} $和$\varTheta_{3}$的重要性测度(表1)如下.
(1) 均值参数反映了随机变量的位置关系. 注意到$\varTheta_{1} $, $\varTheta_{2}$和$\varTheta_{3} $在式(57)中的系数比例为$1:1:\sqrt 2 $, 这与${\cal F}\left({\mu_{1} } \right)$, ${\cal F}\left( {\mu_{2} } \right)$和${\cal F}\left({\mu_{3} } \right)$对应的重要性测度是一致的.这里所表现出的一致性是因为有关系$\mu_{X} =\mu_{1} -\mu_{2} +\sqrt 2 \mu_{3} $.
(2) 标准差参数反映了随机变量的离散程度. 注意到$\sigma_{X}^{2} =\sigma_{1}^{2} +\sigma_{2}^{2} +2\sigma_{3}^{2} $, 而${\cal F}\left( {\sigma_{1} } \right)$, ${\cal F}\left( {\sigma_{2} } \right)$和${\cal F}\left({\sigma_{3} } \right)$对应的重要性测度比例为$1:1:2$, 二者也是一致的.
除了通过表1计算所得的重要性测度外, 直接观察图3和图4中整体灵敏度指标的峰值点亦能得到上述相关结论, 即$\varTheta_{1}$, $\varTheta_{2} $和$\varTheta_{3} $的重要性测度关于均值参数有比例关系$1:1:\sqrt 2 $、而关于标准差参数有比例关系$1:1:2$.
由式(38)所定义的方向灵敏度绘于图5和图6中. 从图5可见, ${\cal F}\left( {\mu_{1} } \right)$与${\cal F}\left( {\mu_{3} }\right)$的方向是相同的、与${\cal F}\left( {\mu_{2} }\right)$的方向则是相反的. 而在图6中, ${\cal F}\left( {\sigma_{1} }\right)$、${\cal F}\left( {\sigma_{2} } \right)$和${\cal F}\left( {\sigma_{3} } \right)$的方向一致. 主要原因如下:
图5
新窗口打开|下载原图ZIP|生成PPT图5相对于均值参数的基于Fréchet导数的整体灵敏度指标的方向灵敏度(算例一)
Fig. 5Direction sensitivity of the Fréchet derivative based GSI in terms of the mean values (Example 1)
图6
新窗口打开|下载原图ZIP|生成PPT图6相对于标准差参数的基于Fréchet导数的整体灵敏度指标的方向灵敏度(算例一)
Fig. 6Direction sensitivity of the Fréchet derivative based GSI in terms of the standard deviations (Example 1)
(1) 均值参数反映了随机变量的总体趋势. 这意味着, 若$\mu_{1} $或$\mu_{3}$有一个正增量$\Delta \mu $, 则在$x<0$范围内的概率密度函数会相应地减小,而在$x\geqslant 0$范围内的概率密度函数则会相应地增大,最终导致系统响应的概率密度函数整体向右移动; 反之则反. 这一结论,正是由物理方程(57)中$\varTheta_{1} $, $\varTheta_{2} $和$\varTheta_{3}$前的符号“$+$”、“$-$”和“$+$”所决定的.
(2) 对于标准差参数而言, 由图6可知, 增大$\sigma_{1} $, $\sigma_{2}$或$\sigma_{3} $均会使得$x\in \left( {-2,+2}\right)$范围内的概率密度函数相应地减小, 而$x\in \left( {-\infty ,-2}\right]\cup \left[ {+2,+\infty } \right)$范围内的概率密度函数会相应地增大,最终导致系统响应的概率密度函数变得更为扁平. 换言之,本例中系统响应标准差将随着$\sigma_{1} $, $\sigma_{2} $或$\sigma_{3}$的增大而增大. 显然, 由$\sigma_{X}^{2} =\sigma_{1}^{2} +\sigma_{2}^{2}+2\sigma_{3}^{2} $可知这一结论是正确的.
4.1.2 与矩独立灵敏度指标的关系
考虑随机系统(57)的矩独立灵敏度指标[10], 对$i=1,2,3$, $s_{i} \left({\varTheta_{i} =\mu_{i} } \right)$分别为
而由式(50)计算可得$\hat{{s}}_{1} =0.121,0$, $\hat{{s}}_{2}=0.121,0$和$\hat{{s}}_{3} =0.242,0$, 满足不等式(49). 值得注意的是,这里的$\hat{{s}}_{j} $依然保留了比例关系$\hat{{s}}_{1} :\hat{{s}}_{2}:\hat{{s}}_{3} =1:1:2$, 而$s_{i} $及矩独立指标$\mathbb{E}_{\varTheta_{i} } \left[{s_{i} } \right]$(计算可得0.185,0, 0.185,0和0.305,7)则显然不具备这一比例特性.
需注意的是, 虽然本例中分布参数$\mu_{1} $和$\mu_{2} $的重要性测度相同,但它们的方向灵敏度是相反的. 可见, 单凭重要性测度提供的灵敏度信息是不够充分的,这侧面反映了基于Fréchet导数的整体灵敏度指标的优越性.
4.2 算例二: 隧道锚杆支护系统
考虑一隧道锚杆支护系统[27,36], 如图7所示, 假定系统是轴对称与各向同性的. 在岩体压力$P_{{rock}}$、喷浆混凝土衬砌承压力$P_{{shot}} $与锚杆承压力$P_{{bolt}}$的共同作用下, 该支护系统的功能函数可表达为[36]图7
新窗口打开|下载原图ZIP|生成PPT图7隧道锚杆支护系统
Fig. 7Rock bolting system of the tunnel
式中, 岩体压力$P_{{rock}} $可隐式表达为[36]
其中$P_{0} $为地应力, $\tilde{{c}}$为修正后岩体有效黏聚力, 有
喷浆混凝土衬砌承压力$P_{{shot}} $为径向变形的函数
其中, $K_{{s}} $为支护刚度模量(考虑$d\ll r_{0} )$
式(62)中$u(\cdot )$为径向变形场,它是离隧道中心距离$r$的反比例函数, 即
锚杆承压力$P_{{bolt}} $按下式计算
式(59),$\sim \!$式(65)中各确定性模型参数取值及其物理意义详见表2.
Table 2
表2
表2隧道锚杆支护系统确定性参数(算例二)
Table 2
 |
新窗口打开|下载CSV
此外, 考虑式(59),$\sim \!$式(65)中的地应力$P_{0}$ (MPa)由四阶矩正态变换方法进行表征[27], 其前四阶矩分别为$\mu_{P_{0} } =9.975$, $\sigma_{P_{0} } =2.300$, $\alpha_{P_{0} } =-0.249,1$和$\beta_{P_{0} } =2.874,8$;
岩体弹性模量$E$ (GPa)服从对数正态分布, 其均值为$\mu_{E} =4.370$, 标准差为$\sigma_{E} =0.524,4$; 初始径向位移$u_{0}$ (cm)服从对数正态分布, 其均值为$\mu_{u_{0} } =3.200$,标准差为$\sigma _{u_{0} } =0.400$.
通常, 对数正态分布的概率密度函数为
因此, 实际分析中对数正态分布函数的参数$a$和$b$按下式进行换算[18]
其中$\mu $为均值、$\sigma $为标准差. 记$E$的分布参数为$a_{E} $和$b_{E} $,$u_{0} $的分布参数为$a_{u_{0} } $和$b_{u_{0} } $, 由式(67)计算得到.本例灵敏度分析结果绘于图8和图9之中.
图8
新窗口打开|下载原图ZIP|生成PPT图8相对于参数$P_{0} $前四阶矩的基于Fréchet导数的整体灵敏度指标(算例二)
Fig. 8Fréchet derivative based GSI in terms of the first fourth-moments of the parameter $P_{0} $ (Example 2)
图9
新窗口打开|下载原图ZIP|生成PPT图9相对于参数$E$和$u_{0} $的分布参数$a$和$b$的基于Fréchet导数的整体灵敏度指标(算例二)
Fig. 9Fréchet derivative based GSI in terms of the distribution parameters $a$ and $b$ of the parameters $E$ and $u_{0} $ (Example 2)
本例中的8个分布参数的重要性测度排序为
对应的重要性测度分别为0.910,5, 0.793,8, 0.713,6, 0.635,3, 0.275,9, 0.161,1,0.122,5和0.086,2, 且均满足不超过1的性质. 可见, 地压力$P_{0}$对功能函数的影响是最大的, 次之为岩体弹性模量$E$和初始径向位移$u_{0} $. 另一方面, 从图8和图9可以进一步观察到:
(1) 相较于低阶矩而言, 高阶矩参数对响应概率密度函数的影响更为复杂, 即其方向灵敏度有更多的不相交正值域或负值域.
(2) 若仅考虑均值分布参数, 可以看到参数$P_{0} $将使得$g$的概率密度函数向正向移动, 而参数$E$则使$g$的概率密度函数向负向移动.这一点可以在锚杆承压力计算式(65)与径向变形场式(64)中得以理解: 当$P_{0}$增大或$E$减小时, 径向变形增大, 从而使得锚杆承压力$P_{{bolt}} $增大, 最终导致$g$增大, 因而其概率密度函数向正向移动.
4.3 算例三: 大坝下稳态受限渗流模型
考虑一个大坝下稳态受限渗流分析模型[14,37], 如图10所示. 在分析中, 假定渗流是各向同性和均匀的、且渗流仅发生在$AB$段和$CD$段.记粉质碎石层的透水率为$T_{{g},x} (x$方向)和$T_{{g},y}$ $(y$方向); 粉质砂土层的透水率为$T_{{s},x}$ $(x$方向)和$T_{{s},y}$ $(y$方向). 记$AB$段的水头差为$h_{{W}} =h_{{D}} +20$, 则$AB$段到$CD$段的渗流控制方程为[14]图10
新窗口打开|下载原图ZIP|生成PPT图10大坝下稳态受限渗流模型的原理图[14]
Fig. 10Schematic representation of the steady-state confined seepage below the dam [14]
本文用有限元法(共计3413个节点, 1628个单元)求解式(69)所示控制方程[14]得到分析区域的水头差分布函数$h_{{W}}\left( {x,y} \right)$, 最后计算渗流量$q$
其中$q$的单位为$L$/(h$\cdot$m). 式(69)各参数取值列于表3中, 其中${\cal U}\left( {7,10} \right)$为区间$\left[{7,10} \right]$上均匀分布, ${\cal L}{\cal N}\left( {\mu ,\sigma }\right)$为均值为$\mu $、标准差为$\sigma $的对数正态分布.
Table 3
表3
表3大坝下稳态受限渗流模型参数(算例三)
Table 3
 |
新窗口打开|下载CSV
分析时, 记$T_{{g},x} $, $T_{{g},y} $, $T_{{s},x} $和$T_{{s},y} $的分布参数分别为$\left( {a_{{g},x} ,b_{{g},x} } \right)$, $\left( {a_{{g},y} ,b_{{g},y} } \right)$, $\left( {a_{{s},x} ,b_{{s},x} } \right)$和$\left( {a_{{s},y} ,b_{{s},y} } \right)$; 对于蓄水高度参数$h_{{D}} $, 考虑它的$\varepsilon $-等效分布为Beta分布$β$$\left( {a_{h} ,b_{h} } \right)$, 其中$a_{h} =1$和$b_{h} =1$. 由于Beta分布定义于区间$\left[ {0,1} \right]$中, 不妨令
式中${B}\left( {a_{h} ,b_{h} } \right)=\Gamma \left( {a_{h} }\right)\Gamma \left( {b_{h} } \right)/\Gamma \left( {a_{h} +b_{h} }\right)$, $\Gamma (\cdot )$为Gamma函数.
灵敏度分析结果绘于图11与图12中. 由表4可知, 按参数$a$进行重要性测度排序有
图11
新窗口打开|下载原图ZIP|生成PPT图11相对于参数$a$的基于Fréchet导数的整体灵敏度指标(算例三)
Fig. 11Fréchet derivative based GSI in terms of the parameter $a$ (Example 3)
图12
新窗口打开|下载原图ZIP|生成PPT图12相对于参数$b$的基于Fréchet导数的整体灵敏度指标(算例三)
Fig. 12Fréchet derivative based GSI in terms of the parameter $b$ (Example 3)
Table 4
表4
表4基于Fr′echet 导数的整体灵敏度指标的重要性测度(算例三)
Table 4
 |
新窗口打开|下载CSV
而按参数$b$进行重要性测度排序有
这意味着: (1)粉质砂土层透水率是非常重要的参数、而粉质碎石层透水率则是相对而言重要性较低;(2) 除在$h_{{D}} $中二者接近外, 参数$a$相较于参数$b$更为重要.但同样所有变量的重要性测度也均不超过1.
4.4 算例四: 钢筋混凝土平面框架结构
考虑一榀钢筋混凝土平面框架结构, 见图13. 结构底层柱截面为600 mm$\times$800 mm、二层柱截面为500 mm$\times$700 mm, 其他层柱截面为500 mm$\times$600 mm;中跨梁截面统一设计为300 mm$\times$450 mm、边跨梁截面均采用300 mm$\times$600 mm.截面配筋信息详见文献[17]. 采用OpenSEES软件对其进行有限元建模, 具体而言: (1) 单元类型设置为纤维梁单元; (2)采用与混凝土结构设计规范[38]对应的混凝土单轴损伤本构模型(“ConcreteD”命令)模拟混凝土材料的应力-应变关系;(3) 考虑各向同性强化的Giuffe-Menegotto-Pinto模型[17]模拟钢筋材料的应力-应变关系(“Steel02”命令).图13
新窗口打开|下载原图ZIP|生成PPT图13钢筋混凝土平面框架结构[17]
Fig. 13Reinforced concrete frame structure[17]
分析中, 考虑第1层, 第2层, 第3$\sim$6层及第7$\sim$10层的混凝土材料的抗压强度$f_{{c}1} $, $f_{{c}2} $, $f_{{c}3}$和$f_{{c}4} $为独立同分布的随机变量, 其概率信息详见表5. 注意到混凝土弹性模量$E_{{c}}$ ($10^{4}$ MPa)与其抗压强度$f_{{c}}$ (MPa)是相关的. 本例中, 采用混凝土结构设计规范中建议的经验关系[38], 即
Table 5
表5
表5混凝土抗压强度参数(算例四)
Table 5
 |
新窗口打开|下载CSV
考察该结构在南北向El Centro地震波作用下的随机动力响应, 地震动幅值调整为中震作用下对应的0.2 g.记感兴趣量为顶层位移响应$d_{{top}} \left( t \right)$的极值, 即$d_{{ext}} =\left\| {d_{{top}} \left( t \right)}\right\|_{\infty } $, 其中$\left\| \cdot \right\|_{\infty } $为$\infty $范数.
表中, ${\cal N}\left( {\mu ,\sigma } \right)$代表正态分布, $\mu $为均值,$\sigma $为标准差, 并对$i=1,2,3,4$, 记$f_{{c}i} $的均值参数为$\mu_{i}$, 标准差参数为$\sigma_{i} $.
本例分析结果如图14和图15所示, 所得整体灵敏度指标的重要性测度列于表6之中.同样, 计算所得的各重要性测度均不超过1. 可见, 与$f_{{c}2}$和$f_{{c}3} $相比, $f_{{c}1} $和$f_{{c}4}$对顶层位移响应极值的影响更大, 且尤以$f_{{c}4} $影响最大. 有意思的是,对底层位移做灵敏度分析, 即取$d_{{ext}} =\left\| {d_{{bottom}}\left( t \right)} \right\|_{\infty } $,此时相对于均值与标准差参数的重要性测度如表7所示. 此时,底层位移响应极值几乎完全由$f_{{c}1} $所影响.表8给出了在0.1$g$峰值加速度下底层位移与顶层位移的重要性测度(此时结构响应基本处于线性状态).值得注意的是, 此时除$f_{{c}1} $外, $f_{{c}3} $和$f_{{c}4}$均对底层位移响应极值的重要性测度亦有不可忽略的影响.
图14
新窗口打开|下载原图ZIP|生成PPT图14相对于均值的基于Fréchet导数的整体灵敏度指标(算例四)
Fig. 14Fréchet derivative based GSI in terms of the mean value (Example 4)
图15
新窗口打开|下载原图ZIP|生成PPT图15相对于标准差的基于Fréchet导数的整体灵敏度指标(算例四)
Fig. 15Fréchet derivative based GSI in terms of the standard deviation (Example 4)
Table 6
表6
表6基于Fréchet 导数的整体灵敏度指标的重要性测度(算例四/顶部位移/0.2g)
Table 6
 |
新窗口打开|下载CSV
Table 7
表7
表7基于Fréchet 导数的整体灵敏度指标的重要性测度(算例四/底部位移/0.2g)
Table 7
 |
新窗口打开|下载CSV
Table 8
表8
表8基于Fréchet 导数的整体灵敏度指标的重要性测度(算例四/0.1g)
Table 8
 |
新窗口打开|下载CSV
从定性上看, 这可能是由于结构在中震作用下进入了一定程度的非线性[11],且非线性主要表现在混凝土材料中, 如图16所示.这也再次表明随机系统的整体灵敏度与物理系统的复杂程度密切相关,且在非线性与随机性耦合情况下, 灵敏度指标排序可能出现与线性情况下不同的结果.
图16
新窗口打开|下载原图ZIP|生成PPT图16某底层柱混凝土材料的应力-应变曲线 (算例四)
Fig. 16Stress-strain curves of concrete materials of one bottom column (Example 4)
5 结论
基于随机系统中不确定性传播的泛函观点及其泛函表达,引入了基于Fréchet导数的整体灵敏度指标. 进一步, 通过定义$\varepsilon$-等效分布, 讨论了该整体灵敏度的参数化表达形式,定义了其对应的重要性测度与方向灵敏度, 并论述了其物理意义.结合概率密度演化理论(PDEM)和概率测度变换(COM),具体给出了该整体灵敏度对应的数值计算方法.通过4个算例对该整体灵敏度的工程实用性进行了分析与论证. 本文主要结论如下:(1) 与传统整体灵敏度不同的是,本文所提整体灵敏度可看作是灵敏度分析从确定性系统到随机性系统的一个自然延拓,因而更具有物理指导意义.
(2) 特别地, 本文所提整体灵敏度不仅能反映系统响应空间中任一点处的灵敏度,亦能通过其重要性测度给出宏观的重要性排序、或通过其方向灵敏度标定特定的设计区域.而且, 通过理论证明重要性测度不超过1. 因此,该整体灵敏度所包含的不确定性量化信息十分丰富.
(3) 概率密度演化理论与概率测度变换具有良好的计算精度与分析效率,可用于解决实际工程的灵敏度分析问题.
此外,如何根据实测数据进行灵敏度分析、如何确保高维概率空间下PDEM-COM方法的计算精度、以及如何考虑非精确概率模型下的灵敏度分析等一系列问题,还有待进一步深入研究.
参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
[本文引用: 3]
[本文引用: 3]
[本文引用: 3]
[本文引用: 3]
[本文引用: 3]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 3]
[本文引用: 2]
[本文引用: 2]
[本文引用: 3]
[本文引用: 1]
[本文引用: 6]
[本文引用: 1]
[本文引用: 6]
[本文引用: 8]
陈建兵, 彭勇波, 刘威等, 译.
[本文引用: 3]
Transl. Chen Jianbing, Peng Yongbo, Liu Wei, et al.
[本文引用: 3]
[本文引用: 2]
[本文引用: 4]
[本文引用: 10]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 5]
[本文引用: 1]
[本文引用: 2]
[本文引用: 4]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
DOIURL [本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 4]
[本文引用: 2]
[本文引用: 2]
[本文引用: 2]
