可靠度分析在结构设计和评估等领域占有重要地位．可靠度分析已经发展多种分析方法，主要包括一次可靠度分析方法、二次可靠度分析方法[1]、蒙特卡罗法[1]、响应面法[2]和矩方法[3-5]等．相比较而言，矩方法由于原理简单、易于实现，且不依赖于验算点，因此受到了较多的关注，其中在基于立方正态变换[5]的矩方法中各参数可由统计矩解析确定，应用最为方便．不难发现，上述方法均是以变量的概率信息已知为前提，然而实际工程中获取变量的完备概率信息并不总是很方便，例如：有的变量仅已知低阶统计矩，有的变量仅已知样本数据，甚至是少量样本数据．对于此类概率信息不完全系统，相应的可靠度研究颇为缺乏．概率信息不完全系统可靠度分析的核心在于如何由不完备的概率信息近似得到完备的概率信息．对于已知变量统计矩的概率信息不完全系统，立方正态变换能够较为方便地建立变量的近似概率信息[6]．当统计矩已知时，只能采用矩匹配原则确定系数，即矩匹配法(PM)[3]；但当样本数据已知时，既可以先获得统计矩再由PM确定系数，又可以先从样本获得其他的概率信息，然后再根据误差最小原则确定系数，如最小二乘法(LS)[7]、线性矩法(LM)[8]等．对比上述方法可知LS与LM较PM更优[9-10]．然而基于立方正态变换近似已知样本数据变量的完全概率信息的研究大多针对单变量情形，双变量或多变量情况研究较少．文献[11]尽管考虑了双变量情形，但是当确定各变量的立方正态变换系数时采用了线性矩，等效相关系数的确定却以传统的统计矩(泊松相关系数)为基础；文献[12]分别采用单变量和双变量的线性矩相关系数，给出了立方正态变换系数及等效相关系数的确定方法，但是等效相关系数对于变量排序的依赖性不利于方法的推广应用．为此，本研究在文献[6]的概率不完全系统的基础上进一步拓展，针对已知少量样本数据或统计矩的变量，介绍了立方正态变换系数的求解方法，同时给出了不依赖于变量排序的双变量等效相关系数的确定方法．以此为基础，结合自适应降维近似模型和直接积分法，发展了概率信息不完全系统可靠度分析的一种简单易行的方法．1 概率不完全系统完全概率信息近似为不失一般性，结构的响应函数Z为Z=GΘ，(1)式中：Θ={Θ1，Θ2，…，Θs}为系统的基本随机向量，其中s为变量个数；G(∙)为反映Z和Θ关系的函数．根据概率信息是否完全，各变量可以分为3类，详见表1．10.13245/j.hust.211115.T001表1单变量的概率信息分类分类变量信息描述①已知少量样本数据的信息不完全变量②已知前几阶统计矩的信息不完全变量③已知边缘概率密度函数的信息完全变量对于双变量情形，可以分为6种子类，详见表2，对于已知联合分布的情况，实际上可转换为表2中的子类VI．对于多变量情形，可视为表2中6种子类的不同组合，不再赘述．10.13245/j.hust.211115.T002表2双变量的概率信息分类分类组合分类组合I①+①IV②+②+ρII①+②+ρV②+③+ρIII①+③+ρVI③+③+ρ注：ρ为变量之间的相关系数矩阵．2 概率信息不完全变量的完全概率信息近似2.1　单变量的完全概率信息近似2.1.1　已知样本数据若已知Θi的n个样本数据，即分类①，则根据Θi的次序统计量，可计算其r阶线性矩[8,12]为λr,i．对Θi引入立方正态变换，即Θi=a0,i+a1,iXi+a2,iXi2+a3,iXi3，(2)式中：Xi为标准正态变量：al,i(l=0，1，2，3)为待定系数，可由线性矩匹配原则确定[13]，即Ai=B-1M-1λi，其中，Ai=(a0,i，a1,i，a2,i，a3,i)T，λi=(λ1,i，λ2,i，λ3,i，λ4,i)T，B=1000-12001-660-112-3020，M=M0,0M0,1M0,2M0,3M1,0M1,1M1,2M1,3M2,0M2,1M2,2M2,3M3,0M3,1M3,2M3,3，Mr,s=∫xisΦr(xi)ϕ(xi)dxi，Φ(∙)和ϕ(∙)分别为标准正态累积分布函数和概率密度函数．2.1.2　已知样本统计矩若已知Θk的均值μk、标准差σk、偏度系数α3,k和峰度系数α4,k，即分类②，则可引入立方正态变换[5]，即Θk=μk+σk(b0,k+b1,kXk+b2,kXk2+b3,kXk3)，式中bl,k(l=0，1，2，3)由如下确定：lk=(6α4,k-8α3,k2-14-2)/36;b0,k=-b2,k=-α3,k/[6(1+6lk)];b1,k=(1-3lk)/(1+b0,k2-lk2);b3,k=lk/(1+b0,k2+12lk2). (3)2.2　双变量的完全概率信息近似2.2.1　子类I若Θi和Θj已知少量样本数据，则Θi与Θj的线性矩相关系数η¯ij可计算[14]为η¯ij=λ2,ij/λ2,i，(4)式中λ2,ij为变量Θi和Θj的联合线性矩，可通过对应变量的伴随量计算确定．若对于Θi和Θj均引入立方正态变换，即式(2)，则Θi与Θj线性矩相关系数η˜ij的理论表达式[12]为η˜ij=ρ[ij]*(c1,i+3c3,i-c3,iρ[ij]*2/2），(5)式中：ρ[ij]*为Θi与Θj对应的标准正态变量Xi和Xj的等效相关系数；c1,i和c3,i可通过文献[12]确定．由式(4)和式(5)确定的线性矩相关系数相等即可给出等效相关系数ρ[ij]*所须满足的方程，即λ2,ij/λ2,i=ρ[ij]*(c1,i+3c3,i-c3,iρ[ij]*2/2）．若调换Θi和Θj的排序，则对应的等效相关系数ρ[ji]*应满足λ2,ji/λ2,j=ρ[ji]*(c1,j+3c3,j-c3,jρ[ji]*2/2）．显然，由上述两式确定的ρ[ij]*与ρ[ji]*并不总是相等．换言之，以线性矩相关系数为基础确定的等效相关系数与变量的排序有关．为避免等效相关系数对变量排序的依赖性，本研究以ρ[ij]*与ρ[ji]*的平均值定义基于线性矩相关系数确定的等效相关系数，即ρij*=(ρ[ij]*+ρ[ji]*）/2．2.2.2　子类II若Θi已知少量样本数据，Θj已知前四阶矩，且两者的相关系数为ρij，则类似于文献[15]，基于Mehler公式可得到ρij与ρij*的关系为ρij=1σi∑k=131k!IP,k,iIP,k,jρij*k，式中：σi=a1,i2+6a1,ia3,i+2a2,i2+15a3,i2；IP,1,i=a1,i+3a3,i；IP,1,j=b1,j+3b1,j；IP,2,i=2a2,i；IP,2,j=2b2,j；IP,3,i=6a3,i；IP,3,j=6b3,j；bl,j可由式(3)确定．2.2.3　子类III若Θi已知少量样本数据，Θj已知累积分布函数Fj(∙)以及ρij，则ρij*可由如下公式确定，即ρij=1σi∑k=131k!IP,k,iIR,k,jρij*k，式中IR,l,j=∫-∞+∞[(Fj-1(Φ(u))-μj)/σj]Hl(u)ϕ(u)du，其中Hl(∙)为l阶埃尔米特多项式．2.2.4　其余子类关于子类IV、子类V和子类VI，等效相关系数ρij*的确定详见文献[15]．2.3　概率信息不完全系统的独立标准正态变换结合2.1节和2.2节，概率信息不完全系统的各变量可用标准正态变量表示，即Θ1,Θ2,⋯,ΘsT=GN(U)，(6)式中：U为独立的标准正态向量：GN(∙)为Θ和U的变换关系．显然，GN(∙)是文献[6]中广义Nataf变换的进一步拓展．3 概率信息不完全系统可靠度分析的矩方法将式(6)代入式(1)，有Z=G(Θ)=G(GN(U))≜gU．3.1　基于双变量降维近似模型的统计矩估计引入g(U)的自适应双变量降维近似[16]，Z的r阶统计矩可表示为   MZ,r=E[gU-C]r≈E{∑ij[Ih,ijg(Ui, Uj,  u̲ij,c)+(1-Ih,ij)(g(Ui,  u̲i,c)+g(Uj,  u̲j,c)-g(uc))]-s-2∑k=1sg(Uk,  u̲k,c)+s-1s-2∙g(uc)/2-C}r, (7)式中：当r=1时，C=0，MZ,1为Z的均值；当r1时，C=MZ,1；MZ,r为Z的r阶中心矩；uc={0，0，⸱⸱⸱，0}为参考点；ui1i2⸱⸱⸱il,c为除变量Ui1，Ui2，⸱⸱⸱，Uil外其余变量对应的参考点坐标；Ih,ij为Ui和Uj的交叉项示性函数，其判定方式可由文献[16]确定．采用直接积分法[17]，式(7)可改写为   MZ,r≈∑t1=1l∑t2=1l⋯∑tn=1l∏i=1swGH,tiπ{∑ij[Ih,ijh(2uGH,ti,2uGH,ti,  u̲ij,c)+(1-Ih,ij)(h(2uGH,ti,  u̲i,c)+h(2uGH,tj,  u̲j,c)-h(uc))]-n-2∑k=1sh(2uGH,tk,u̲k,c)+n-1n-2h(uc)/2-C}r,式中：uGH,tk和wGH,tk分别为高斯-埃尔米特求积公式的第tk个节点坐标和权系数；l为求积节点的数量，本研究中取l=7．3.2　基于立方正态变换的可靠度矩方法由MZ,3和MZ,4可确定Z的偏度系数α3,Z和峰度系数α4,Z，再结合式(3)可确定系数b0,Z，b1,Z，b2,Z和b3,Z，于是概率不完全系统的可靠指标β [3]为β=23p-q+Δ3--q+Δ323+l13k2，式中：Δ=q2+4p3;p=(3b1,Zb3,Z-b2,Z2)/(9b3,Z2); q=(2b2,Z3-9b1,Zb2,Zb3,Z)/(27b3,Z3)+(μZ/σZ-b2,Z)/b3,Z．4 算例分析算例1 考察已知N个样本数据的变量Θ．为便于验证，样本数据来自于均值为10、标准差为1的对数正态分布，其偏度系数和峰度系数分别为0.301和3.161 5．为便于阐述，本研究将数据来源变量称为对照变量．对于Θ的N个样本数据，分别采用PM，LM和LS计算立方正态变换中的系数，并根据文献[18]计算Θ的前四阶矩．将上述三种方法得到的前四阶矩与对照变量的前四阶矩进行对比，表3给出了各方法计算Θ的均值、标准差、偏度和峰度的相对差异，即εμ，εσ，εα3和εα4．10.13245/j.hust.211115.T003表3前四阶矩的结果比较N方法εμεσεα3εα410PM0.4259.92447.67433.248LM0.4250.1916.62314.461LS0.16216.61727.98312.30020PM0.1513.47141.25211.525LM0.1513.25117.31123.412LS0.03124.80217.617111.38730PM0.1870.90368.11314.257LM0.1873.59230.8787.648LS0.12915.20931.52742.10840PM0.0490.0527.0458.316LM0.0493.6647.37812.426LS0.05113.39546.25946.82150PM0.0731.14525.46915.399LM0.0741.5967.4751.571LS0.0037.34210.30311.797%从表3可以发现：a．总体而言，随着样本容量的增加，由PM，LM和LS确定的前四阶矩与对照变量的相对差异呈减小趋势；b．对于少样本情形，相比于PM和LS，由LM确定得到的统计矩与对照变量的相对差异更小．表3仅是一次样本数据的结果，为了考察样本随机性的影响，对本算例重复计算9次，所得结论与本算例类似．此外，本研究还分别以极值I型分布、韦布尔分布和卡方分布等常见分布为对照变量，比较了PM，LM和LS三种方法的性能，结论与本算例类似．由于篇幅原因，不再赘述．综上，在少量样本数据的情况下，建议采用LM计算立方正态变换的系数．算例2 考察图1所示结构地基A点处沉降的可靠度，其功能函数为Z=G(p0,Δp,Cc,e0,H)=2.5/H-[Cc/(1+  e0)]lg(1+Δp/p0)，式中：Cc为土壤压缩指数；e0为沉降之前的黏土层孔隙率；H为黏土层厚度；p0为B点加载前的原始有效应力；Δp为结构物建造之后在B点所增加的有效应力．Cc和H分别为正态分布和对数正态分布，其均值分别为0.396和120，其标准差分别为0.099和6，p0和Δp均已知50个样本数据，e0的均值、标准差、偏度和峰度分别为1.09，0.163 5，1.139 547和5.4．变量p0，Δp，Cc，e0和H的相关系数矩阵为ρ=1.0*0.10.10.1*1.00000.101.00.40.450.100.41.00.30.100.450.31.0  ，式中*表示相关系数未知．10.13245/j.hust.211115.F001图1土壤沉降模型显然，本算例属于多变量概率信息不完全的情况，其中：p0和Δp属于子类I；p0和e0属于子类II；p0和Cc，p0和H属于子类III；Cc和e0，e0和H属于子类V；Cc和H属于子类VI；子类Ⅳ已在文献[6]中讨论，这里不再赘述．现有方法无法求解此可靠度问题，可采用本研究方法求解．由于p0的样本数据来自于均值为3.72和标准差为1.116的对数正态分布，Δp的样本数据来自于均值为0.5和标准差为0.2的极值I型分布，p0和Δp的等效相关系数为0.3，e0的前四阶矩来自于极值I型分布，因此可将p0，Δp和e0具有完整概率信息的情况作为对照工况，以对照工况的蒙特卡罗模拟(MCS)结果为标准解验证建议方法的可行性．首先，确定p0，Δp和e0的立方正态变换系数，为验证系数的合理性，表4给出了p0和Δp的前四阶矩与对照变量前四阶的相对差异，结果表明：均值和标准差的相对差异较小，偏度系数和峰度系数的相对差异略有增加．10.13245/j.hust.211115.T004表4变量的统计矩相对差异变量εμεσεα3εα4p00.2660.17615.12612.786Δp0.2829.22929.22227.475%然后，按照建议方法求解等效相关系数ρ*，即ρ*=1.0000.4340.1000.1030.1030.4341.0000000.10001.0000.4130.4500.10300.4131.0000.3080.10300.4500.3081.000  ．(8)对照工况的等效相关系数ρ¯*为ρ¯*=1.0000.3110.1020.1050.1020.3111.0000000.10201.0000.4120.4500.10500.4121.0000.3080.10200.4500.3081.000 ．(9)对比式(8)和(9)可知：对于子类V和VI，即已知累积分布函数和已知统计矩，本研究提出的等效相关系数求解方法具有较高的精度，而对于子类II和III，即已知样本数据和累积分布函数或统计矩，通过采用本研究提出的等效相关系数求解方法，相对误差最大仅为1.96%，而对于子类I，即已知样本数据，相对误差为39.5%．同时，为了验证本研究提出的等效相关系数的求解方法的可行性，本算例对比了由式(5)确定的等效相关系数ρ[p0,Δp]*(记为EP-1)与ρ[Δp,p0]*(记为EP-2)、平均等效相关系数ρp0,Δp*(记为EP-3)、由p0与Δp的样本统计矩结合矩匹配方法计算的等效相关系数(记为EP-4)以及由文献[11]确定的p0与Δp的等效相关系数(记为EP-5)．由EP-1，EP-2，EP-3，EP-4和EP-5计算等效相关系数的相对差异分别为36.415%，42.593%，39.504%，46.122%和47.319%．通过等效相关系数的相对差异可以明显地看出：相对EP-4和EP-5，EP-3的相对差异较小，同时EP-3的相对差异介于EP-1和EP-2之间，说明变量的排序是对等效相关系数有影响．表5给出了由本文方法得到Z的前四阶矩的计算结果．从表5的结果可以看出：对于概率信息不完全系统，本研究建议的统计矩求解方法是可行的，并且统计矩的最大相对误差为0.574%．本研究建议方法的可靠指标为1.799，采用MCS求解的可靠指标为1.737，相对误差为3.570%．10.13245/j.hust.211115.T005表5算例2的前四阶矩计算结果方法μZσZα3,Zα4,Z分析次数MCS0.0100.005-1.0465.1141×106PPM0.0100.005-1.0405.125391为了考察样本随机性对可靠指标的影响，对本算例重复计算9次，可靠指标的相对误差的范围为0.211%~3.597%，可以看出：样本的随机性对可靠指标的误差有一定影响，但总体上是在可接受的范围内．5 结语针对概率信息不完全系统，本研究首先根据变量的不完全的概率信息进行分类，然后结合立方正态变换，对单变量和多变量的近似完全概率信息进行重建，同时根据变量的不完全信息推导出了等效相关系数的求解方法，并结合直接积分法和双变量自适应降维近似模型，提出了概率信息不完全系统的可靠度分析方法．算例分析表明：本研究方法能够通过变量的有限样本数据建立变量的近似完全概率信息，且平均等效相关系数求解方法是可行的；所建议的可靠度分析方法适用于处理概率信息不完全系统的可靠度分析问题．