一种风荷载效应极值不确定性的分析方法

摘要

本发明公开了一种风荷载效应极值不确定性的分析方法，对非高斯风荷载效应极值在任意分位点处的不确定性进行计算。计算过程中，采用HPM估计母分布；采用转换过程方法得到极值分布；归一化过程极值分布参数估计的经验公式；基于概率分布函数和极值任意分位点的高斯推理，提出了计算极值任意分位点的不确定性的两种分析方法。本发明在考虑气动随机性的情况下，估算非高斯风荷载效应极值在任意分位点处的不确定性，对结构抗风可靠度设计有重要意义。

说明书

技术领域

本发明属于结构抗风设计领域，具体而言是在考虑气动随机性的情况下，两种非高斯风荷载效应极值在任意分位点处的不确定性的分析方法。

背景技术

在结构抗风可靠度设计中，风压以及风致响应在给定时距内的极值确定是十分重要的工作。在过去的几十年里，大量的研究方法被提出以确定风荷载效应的峰值因子或者某一分位点的极值(e.g.,Davenport 1964；Sadek和Simiu 2002；Chen和Huang 2009；Kwon和Kareem 2011；Huang et al.2013；Yang et al.2013；Ding和Chen 2014；Huang etal.2016b；Ma et al.2016)。其中，基于Hermite多项式模型(Hermite polynomial model，HPM)的转换过程方法由于可以根据偏度和峰度对函数形状进行简便灵活的调整，且具有较高的准确度(e.g.,Peng et al.2014，Huang et al.2016b)，因此得到了广泛的应用。最近，Huang等(Huang et al.(2016a)和Luo和Huang(2016))研究了风压与其极值的相关性结构。

众所周知，由于众多随机因素(如气动力环境、模型选取、校准误差等)的存在，由短期样本(如10min或者1h)估计得到的风压极值和响应极值具有很大的变异性。一般地，在极值估计的不确定性研究中，这些随机因素可以分为两类：偶然因素(变量内在的随机性)和认知因素(缺乏认识或数据)(Kiureghian和Ditlevsen 2009)。Minciarelli等研究了多种认知因素(如风洞试验中的时间间隔、地面粗糙度及其他误差)的不确定性对风荷载效应极值的影响。如今，更多极值不确定性的研究运用风洞试验所测的超长样本。Yang和Tian利用多样本数据，对非高斯风压系数峰值因子的概率模型进行建模。Gavanski等采用Gumbel分布拟合方法，研究了取样时长和极值取样频率的不确定性对风压极值估计的影响。Wu等利用超长试验数据研究有限样本数据所带来的认知因素在各重现期风致响应极值中的影响。

一般地，关于极值不确定性的研究主要针对极值的峰值因子(即Gumbel分布的57％分位点)。然而，其他分位点处极值的研究在结构抗风设计中也是十分必要的。例如，78％分位点的风压系数极值经常用于设计规范中设计风荷载的确定，以考虑风载或风速的极值不确定性(e.g.,Cook 1990；Chen和Huang 2010)。

发明内容

本发明所要解决的技术问题是提供一种风荷载效应极值不确定性的分析方法，有效计算风荷载效应极值任意分位点的不确定性。

为解决上述技术问题，本发明采用的技术方案是：

一种风荷载效应极值不确定性的分析方法，包括以下步骤：

步骤1：采用HPM估计母分布

风压以及风致响应往往具有一定的非高斯性，利用基于HPM的转换过程方法估计风荷载效应的极值分布，即采用HPM估计母分布。

假设非高斯过程Y(t)的均值、标准差、偏度和峰度分别表示为r₁、r₂、r₃和r₄，对其归一化后有X(t)＝[Y(t)-r₁]/r₂。根据HPM，转换到标准高斯过程U(t)：

x＝k[u+h₃(u²-1)+h₄(u³-3u)]>

式中，k，h₃及h₄为该HPM曲线的形状控制参数，通过Newton-Raphson迭代求解非线性方程组得到(e.g.,Ditlevsen>

式中，为标准高斯过程的概率密度函数。由X(t)的概率分布进一步得到Y(t)的概率分布。

对于偏度和峰度位于有效区域的过程，HPM能够很好地拟合其概率密度函数；对于偏度峰度位于有效区域外的过程，可以对偏度或者峰度进行调整，使其落在区域边缘近似估计其极值。

步骤2：采用转换过程方法得到极值分布

当母分布得到后，采用转换过程方法得到其极值(e.g.,Sadek和Simiu 2002)。

已知标准高斯过程U(t)的极值累计分布函数(时距T)为：

式中，λ₀＝ν_0,uT为高斯过程U(t)的零超越次数；ν_0,u为过程U(t)的零超越率，由下式计算：

式中，f为频率(Hz)，S_U(f)为高斯过程U(t)的功率谱密度函数，由于极值对零超越率不是很敏感(e.g.,Sadek和Simiu>U(f)采用非高斯过程Y(t)的功率谱密度函数S_Y(f)代替。得到高斯过程U(t)的极值分布后，通过等概率关系可将其转换得到非高斯过程X(t)的极值X_pk与相应的概率值，详细介绍可参考相关文献(e.g.,Sadek和Simiu>

研究表明，非高斯过程X(t)的极值分布接近于Gumbel分布(e.g.,Holmes和Cochran 2003；Huang et al.2016b)，表示为：

式中，δ_x和ψ_x分别表示Gumbel分布的位置参数和尺度参数。由此，进一步得到原过程Y(t)的极值分布表达式，其分布参数为：

δ_y＝r₁+r₂δ_x>

ψ_y＝r₂ψ_x>

式中，δ_y和ψ_y分别表示过程Y(t)的极值分布的位置参数和尺度参数。注意：如果采用广义极值分布对极值进行拟合，形状参数保持不变。

步骤3：归一化过程极值分布参数估计的经验公式

建立两组变量间(第一组：偏度、峰度和零超越次数；第二组：归一化过程极值分布的位置参数和尺度参数)的直接关系式，主要步骤如下：

首先，确定偏度、峰度和零超越次数的变化范围。根据Huang等(Huang etal.2016a和b)研究，典型建筑屋面风压和屋面板上升力的偏度峰度散点分布如图1所示。图中，虚线表示HPM的有效区域(Winterstein和MacKenzie(2013))。实线所围区域表示被选择的偏度峰度变化范围，其基本上包含实际中所有可能的偏度峰度数据组。基于两栋UF房屋屋面的风压时程(Yang et al.2013)，可以发现零超越率的变化范围为1.44～9.37。在工程实际中，零超越率的变化范围可认为是1～10。一般地，确定极值的时距为10min或者1hour。相应地，零超越次数的变化范围为600～36000。

其次，针对上述第一组变量范围内的所有取值组合，基于HPM转换过程方法求得对应的归一化非高斯过程的极值分布(Gumbel)参数值，作为第二组变量。在特定零超越次数下估计得到的偏度峰度值如图2、图3所示。与尺度参数相比，位置参数也许对较小的零超越次数比较敏感。

最后，通过多元线性回归分析对两组变量进行关系拟合，得到如下经验公式：

式中，δ_x和ψ_x的系数在表1中给出。可见，在已知风荷载效应过程的偏度、峰度和零超越次数的情况下，所得的经验公式能够直接用来估计极值分布。

表1 经验公式中δ_x和ψ_x的系数

为了验证该经验公式的准确性，采用被估计参数的相对误差(RE)进行评估。以δ_x为例，对于给定的偏度、峰度和零超越次数，其RE表示为：

式中，为通过式(8)求得的参数值，δ_x则直接基于HPM转换过程求得。对于给定变化区域内的特定偏度峰度，在不同的λ₀(变化范围：600-36000)下，δ_x的RE最大值如图4所示。同样地，ψ_x的RE最大值如图5所示。可以看出，对于δ_x和ψ_x，RE基本上都在3.5以下。另外，为了进一步评估经验公式的准确性，可采用式(10)定义极值任一分位点的RE。同理，在不同的λ₀下，极值57％和78％分位点的RE最大值如图6和图7所示。可以看出，其RE都小于4％。因此，该经验公式对归一化过程的极值分布参数的估计结果令人满意。

步骤4：极值的不确定性估计

Yang和Tian等(e.g.,Yang和Tian 2015)研究表明，风压过程前四阶矩(如均值、标准差、偏度和峰度)具有显著的变异性或者不确定性。相应地，极值分布也会存在不确定性。尽管Yang和Tian等研究表明零超越次数的不确定性对极值的影响很小，其影响需要进一步的讨论，因为图2显示位置参数对零超越次数的变化比较敏感。

前四阶矩和零超越次数将被看做随机变量以考虑其不确定性。设R₁,R₂,R₃和R₄分别用来表示均值、标准差、偏度和峰度四个随机变量，Λ₀表示零超越次数这个随机变量。

第一种方法：基于概率分布函数的分析方法(方法1)

为了描述前四阶矩和零超越次数的不确定性，通常采用概率分布。假设R₃,R₄和Λ₀的边缘累计分布函数分别为和其相关系数矩阵Σ从大量的样本中估计得到或直接指定。那么，非高斯变量R₃,R₄和Λ₀可以分别通过等概率关系与标准高斯变量联系：

式中，Φ为标准高斯分布。基于雅克比变换，R₃,R₄和Λ₀的联合概率密度函数表示为：

式中，表示三变量联合高斯概率密度函数，其相关系数矩阵为Σ_z，矩阵中元素通过Σ中元素求得，如Z₁和Z₂的相关系数通过下式进行估计：

式中，为R₃和R₄的相关系数；为相应的高斯变量Z₁和Z₂的相关系数；E[·]和D[·]分别表示随机变量的均值和标准差。

经验公式建立了偏度、峰度和零超越次数与位置、尺度参数间的直接联系，通过雅克比变换可以从R₃,R₄和Λ₀的联合概率密度函数中得到极值两参数Δ_x和Ψ_x的联合分布函数。然而，经验公式中高阶项的存在使得反函数求解十分复杂且需要大量的迭代计算，耗费大量时间，降低了不确定分析方法的效率。

为了有效地获取极值分布参数的联合概率函数，借助蒙特卡洛模拟(MCS)，具体方法为：首先，基于Nataf变换并通过蒙特卡洛模拟，得到一组R₃,R₄和Λ₀的模拟样本，其维持变量之间的相关性结构(e.g.,Huang>x和Ψ_x的样本；然后，对模拟得到极值分布中两参数的样本进行概率密度函数拟合，分别表示为和其相关系数从这些样本中计算得到。最后，通过下面的雅克比变换得到Δ_x和Ψ_x的联合概率密度函数

如果Λ₀的随机性对极值的影响较小，采用零超越次数的样本均值λ_0,m进行近似分析。从而，式(12)可简化为：

同样地，通过MCS，Δ_x和Ψ_x的联合概率密度函数可由R₃和R₄的联合概率密度函数求得。

当Δ_x和Ψ_x的联合概率密度函数求得后，Δ_y和Ψ_y也可进一步求得。假设在R₁＝r₁和R₂＝r₂的条件下，Δ_y和Ψ_y的条件联合概率密度函数为根据条件概率可得：

式中，表示Δ_y,Ψ_y,R₁和R₂的联合概率密度函数，表示R₁和R₂的联合概率密度函数，可根据式(12)求得。根据雅克比变换得下式：

式中，雅克比矩阵的秩可由式(6)和(7)推得：

将式(15)到(17)带入下式Δ_y和Ψ_y的联合概率密度函数：

得到：

过程Y(t)在任一分位点q(0<q<1)的极值Y_q表示为：

Y_q＝Δ_y-Ψ_yln(-lnq)>

当得到Δ_y和Ψ_y的联合概率密度函数后，Y_q的累计分布函数表示为：

通过雅克比变换，Y_q的概率密度函数则表示为：

第二种方法：基于极值任意分位点的高斯推理的分析方法(方法2)

虽然方法1具有满意的准确度，当样本数量有限时，此方法也许是不可行的。一般地，母过程的极值分布参数Δ_y和Ψ_y被认为服从独立同一分布。可以证明得到Δ_y和Ψ_y渐进趋于高斯分布。因此，式(20)中的Y_q将服从高斯分布。下面将讨论其均值和标准差的计算方法。

为了简要介绍方法2，这里将不考虑Λ₀的不确定性，也就是，用λ_0,m代替式(8)中的λ₀。在r₃＝E[R₃]和r₄＝E[R₄]时，δ_x(r₃,r₄,λ_0,m)的Taylor展开式表示为：

式中，k＝r₃-E[R₃]和l＝r₄-E[R₄]；Rn为残余项.

在式(23)中取n＝2和1，则Δ_x的均值和方差近似表示为：

式中，Cov(R₃,R₄)为R₃和R₄的协方差。同样地，得到E[Ψ_x]和D[Ψ_x]。

令η(r₃,r₄,λ_0,m)＝δ_x(r₃,r₄,λ_0,m)ψ_x(r₃,r₄,λ_0,m)，在式(24)中，以η⁰＝η⁰(E[R₃],E[R₄],λ_0,m)代替得到其均值。因此，Δ_x和Ψ_x的协方差为：

其中

和

R₁和R₂与母过程相关，Δ_x和Ψ_x则与归一化母过程的极值相关。因此，这两组参数见的相关性通常很弱，R₁和R₂与Δ_x和Ψ_x相互独立。从而，Δ_y和Ψ_y的均值直接由式(6)和(7)求得：

E[Δ_y]＝E[R₁]+E[R₂]E[Δ_x]>

E[Ψ_y]＝E[R₂]E[Ψ_x]>

通过Taylor展开，方差及协方差近似表示为：

D[Δ_y]≈D[R₁]+D[R₂]E[Δ_x]²+D[Δ_x]E[R₂]²>

D[Ψ_y]≈D[R₂]E[Ψ_x]²+D[Ψ_x]E[R₂]²>

Cov(Δ_y,Ψ_y)≈E[Δ_x]E[Ψ_x]D[R₂]+E[Ψ_x]Cov(R₁,R₂)+E[R₂]²Cov(Δ_x,Ψ_x)(31)

通过式(20)的关系式，估计得到Y_q的均值和方差：

E[Y_q]＝E[Δ_y]-ln(-lnq)E[Ψ_y]>

D[Y_q]≈D[Δ_y]+ln²(-lnq)D[Ψ_y]-2ln(-lnq)Cov(Δ_y,Ψ_y)>

与方法1相比，方法2具有简便且有效的优点，其适用于少量样本的情况。另外，在这个方法中，仅仅需要获得多元样本的前四阶矩的均值、方差和协方差。这些统计量可从样本中估计得到，或者从它们的概率密度函数中进行推导得到，或者直接给定。而且，计算Y_q的均值和方差的步骤可被延伸到进一步考虑Λ₀的不确定性的情况。然而，需要注意的是：如果一个或者几个风荷载效应样本由于其他因素(如突然的天气变化和没有检测到的仪器故障)而偏离总体很多，方法2中的矩估计会受到严重的影响，从而导致极值不确定性估计的不准确性。相反地，方法1对这些矩的全概率进行建模，避免了上述缺陷，在极值不确定性估计中具有稳定性。

与现有技术相比，本发明的有益效果是：本发明在考虑气动随机性的情况下，估算非高斯风荷载效应极值在任意分位点处的不确定性，对结构抗风可靠度设计有重要意义。

附图说明：

图1为风荷载效应的偏度峰度分布散点图。

图2为在5种零超越次数下估计得到的位置参数。

图3为在5种零超越次数下估计得到的尺度参数。

图4为不同λ₀下的相对误差最大值(位置参数)。

图5为不同λ₀下的相对误差最大值(尺度参数)。

图6为不同λ₀下的相对误差最大值(57％分位点)。

图7为不同λ₀下的相对误差最大值(78％分位点)。

图8为房屋模型FL30。

图9为屋面及测点布置。

图10为R₁和R₂，R₃和R₄这两组变量的相关系数

图11为R₁和R₂，R₃和R₄这两组变量的相关系数

图12为R₁和R₂，R₃和R₄这两组变量的相关系数(的最大值)。

图13为R₁和R₂，R₃和R₄这两组变量的相关系数(的最大值)。

图14为零超越次数的均值和变化幅度散点图。

图15为均值的直方图。

图16为方差的直方图。

图17为偏度的直方图。

图18为峰度的直方图。

图19为联合概率密度函数

图20为联合概率密度函数

图21为分布位置参数的直方图。

图22为分布尺度参数的直方图。

图23为分布位置参数和尺度参数的联合概率密度函数。

图24为布位置参数和尺度参数的联合概率密度函数。

图25为分布位置参数的边缘概率密度函数。

图26为分布尺度参数的边缘概率密度函数。

图27为极值分位点57％的概率密度函数及其直方图。

图28为极值分位点78％的概率密度函数及其直方图。

图29为风压系数在78％分位点处的极值变异系数。

图30为房屋的框架结构模型。

图31为在分位点57％处，底柱B的极值概率密度函数。

图32为在分位点78％处，底柱B的极值概率密度函数。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明，具体如下：

大量的风压样本将被用来评估所提方法的计算性能，这里，采用该方法计算风压系数和风致结构响应的不确定性。风压数据来源于加拿大西安大略大学(University ofWestern Ontario，UWO)的边界层风洞实验室。采用缩尺比为1:50的试验模型(FL30)，如图8所示(没有周边建筑；在B类地貌下)。屋面共布置474个测压点，其分布情况如图9所示。所用数据规定为120°风向下，采样频率400Hz，采样时长为3h(实验模型尺度上)。如果在实际尺度中10m高度处的风速为31.7m/s，那么速度缩尺比为1:5。相应地，实际尺度中的采样频率和时长分别为40Hz和30h。本发明将30h的超长风压数据分成180段10min数据。该模型的详细信息可参考文献Peng et al.(2014)。为了更好地进行解释，风压数据将乘以-1。

1、风压系数

首先，研究每个测点的风压系数前四阶矩以及零超越次数的相关性；其次，讨论零超越次数的不确定性对极值的影响；最后，评估两种分析方法的计算性能。

(1)前四阶矩以及零超越次数的相关性

对于每个屋面测点，R₁和R₂，R₃和R₄这两组变量的相关系数分别如图10和图11所示。可以看出，均值和标准差之间的相关性不是很强，偏度和峰度之间的相关性在许多情况下都大于0.8。

R₁和R₃,R₁和R₄,R₂和R₃,和R₂和R₄这四组变量的相关性也可计算得到。这些相关系数的最大值如图12所示。可以看出，均值与偏度(峰度)之间以及标准差与偏度(峰度)的相关性很弱，可以忽略不计，证实了Δ_x和Ψ_x可被认为独立于R₁和R₂。其解释如下：因为Δ_x和Ψ_x仅仅是R₃和R₄的函数，而R₁和R₂基本上独立于R₃和R₄，从而可以得出结论：Δ_x和Ψ_x独立于R₁和R₂。

同样地，Λ₀和R₃,和Λ₀和R₄这两组变量的相关系数也可通过计算得到。这些相关系数的最大值如图13所示。可以看出，零超越次数基本上独立于偏度和峰度。

(2)零超越次数的不确定性对极值的影响

在探讨零超越次数的不确定性对极值的影响之前，首先讨论零超越次数的变异特征。零超越次数的均值、最小值和最大值分别表示为λ_0,m,λ_0,min和λ_0,max，可从180段风压数据估计得到。零超越次数的变化幅度λ_0,r＝λ_0,max-λ_0,min。图14展示了474个测压点的零超越次数的均值和变化幅度散点图，发现这两者存在这近似线性关系。计算所有测压点的180段样本的零超越次数的变异系数，结果显示其变异系数的范围为1.3％～5.3％。由此可知，零超越次数的变异性是相对较弱的。

为了研究零超越次数的不确定性对极值的影响，分别在λ_0,m,λ_0,min和λ_0,max下，计算X(t)某一分位点的极值，可分别表示为x_pk,m,x_pk,min和x_pk,max。在极值估计中，采用180段样本的偏度和峰度均值进行计算。结果显示，对于所有测压点而言，在分位点57％和78％下，(x_pk,max-x_pk,m)/x_pk,m和(x_pk,m-x_pk,min)/x_pk,m的较大者都小于1.5％。可以得出结论：零超越次数的不确定对极值的影响不大。因此，在后续分析中将忽略零超越次数的不确定性因素，采用180段数据中零超越次数的平均值进行分析。

(3)两种分析方法的性能评估

图9中的测点A的风压数据将用来说明所提方法的计算性能。从180段风压数据中可计算得到均值、标准差、偏度和峰度。均值和标准差的直方图和其概率密度函数的高斯拟合曲线分别如图15和图16所示。可以看出，均值和标准差基本上服从高斯分布。根据中心极值定理，随着样本大小的增加，均值和标准差会接近高斯分布。偏度和峰度的直方图分别如图图17和图18所示，由于两者都存在着明显的非高斯性，采用HPM拟合概率密度函数，分别如图图17和图18所示。

为了求得均值和标准差以及偏度和峰度的联合概率密度函数，其相关系数应该分别直接由180段风压数据估计得到，其值为：和由于偏度和峰度的概率密度函数是采用HPM进行拟合的，其在式(13)中相应的相关系数可采用简化公式得到(Blaise et al 2016；Luo和Huang 2016)：

式中，k_i,h_3,i和h_4,i为R_i(i＝3,4)的模型参数。通过计算得到相应地，可求得联合概率密度函数和其分别如图19和图20所示。

归一化过程的极值分布位置参数和尺度参数的直方图如图21和图22所示。通过MCS方法，可以得到这两个参数的模拟样本，其概率密度如图21和图22所示。从图中可以看出，模拟样本接近原始数据直方图。另外，模拟得到偏度和峰度的相关系数为0.9932，与由原始数据估得的系数0.9955一致。发现模拟样本非常接近高斯分布，其解释如下：Δ_x和Ψ_x为Δ_y和Ψ_y的线性转换，R₁和R₂作为转换中的变量在概率中将收敛到固定值。因此，Δ_x和Ψ_x也将接近高斯分布。根据拟合得到的高斯边缘分布以及估计得到的相关系数，可以得到归一化过程极值分布位置参数和尺度参数的联合概率密度函数，如图23所示。

根据图23的联合概率密度函数，原始过程Y(t)的极值Gumbel分布两参数的联合概率密度函数可通过式(19)计算得到，如图24所示。两个参数的边缘概率密度函数则如图25和图26所示。可以看出，方法1推导得到的分析性模型能很好地对两个参数的直方图进行拟合。另外，也可从方法2中推导得到两个参数的高斯分析模型，如图25、26所示，拟合效果同样很好。从式(22)估计得到的极值分位点57％和78％的概率密度函数与其相应的直方图分别如图27和图28所示。可以发现其拟合效果很好。另外，从R₁,R₂,R₃和R₄的样本中，在两个分位点处的极值均值和标准差可通过方法2得到，其相应的概率密度函数如图27和图28所示，发现它们能很好地对计算极值的不确定性。

在分位点57％和78％处，测压点A的风压系数的极值变异系数大约为8％。所有测压点的风压系数在78％分位点处的极值变异系数如图29所示，其都大于6.5％，一些达到10％甚至20％。

2、风致响应

基于房屋模型FL30，设计框架结构模型。根据钢结构设计相关规范(GB50017-2003；CECS 102-2002)，进行结构布局，并对梁、柱以及檩条等的材料和截面进行确定。在Sap2000中建立有限元模型，如图30所示。柱底采用固定支座，梁与边柱采用刚接节点，梁与中柱则采用铰接节点。梁和柱采用钢材Q345，檩条则采用钢材Q235。将30h风压(风速为31.7m/s)输入该有限元模型中，得到180段结构响应时程数据。

由于柱底弯矩(base bending moments，BBM)在结构设计中非常重要，采用底柱B的柱底弯矩进行分析计算。与风压系数一样，零超越次数的不确定性将被忽略不计。在分位点57％和78％处，由方法1和2可计算得到极值概率密度函数，分别如图31和图32所示。可以发现，两种方法都能很好对极值响应不确定性进行估计。在两个分位点下，其极值的变异系数都大约为7％。

注：Δ_y和Ψ_y的高斯联合分布

风荷载效应的极值为Y_pk。假设：Y_pk的概率密度函数为这里服从Gumbel分布；θ＝[Δ_y,Ψ_y]；n表示风荷载效应极值的个数，其值足够大。参数的对数似然函数为：

参数的最大似然估计可在以下约束中得到：

式中，在Y_pk,i，i＝1,2,…,n服从独立同一分布的前提条件下，是存在的并收敛于θ₀。

在θ₀处，式(36)的一阶Taylor展开为：

式中，为当θ＝θ₀时的值；Rn'为残余项.令根据弱大数定理可知,在概率上收敛于B(Chung 1974)。

从而可得：

注意和根据中心极值定理，当n→∞时，

式中，表示具有零均值，方差的高斯分布，由此可得：

也就是说，θ＝[Δ_y,Ψ_y]服从两变量高斯分布。