随机响应分析

面向多敏感值的个性化随机响应机制设计与分析宋海娜;罗涛;韩新宇;李剑峰【摘要】在实际数据收集中,不同敏感值的敏感度有很大差异,隐私保护需求也不相同.然而,现有的基于随机响应的本地化隐私保护模型针对所有敏感值都执行同样程度的隐私保护,从而可能造成某些低敏感度的敏感值过度保护,而某些高敏感度的敏感值却保护不足.基于此,本文在常规随机响应(Conventional Randomized Response,CRR)模型的基础上,考虑个性化的隐私需求,引入敏感值权重,并将其引入到随机响应的决策中,提出一种面向多敏感值的个性化随机响应(Personalized Randomized Response,PRR)机制,该机制能够确保不同的敏感值群体均能达到各自期望的隐私保护程度,实现个性化的隐私保护.理论分析和仿真实验表明,在机制的主观隐私泄露程度一定时,相比于CRR模型,本文所提的PRR机制统计估计误差更小,即获得的统计数据的质量更高,同时又保证了个性化的隐私保护.【期刊名称】《电子学报》【年(卷),期】2019(047)006【总页数】8页(P1236-1243)【关键词】随机响应;敏感值权重;主观隐私泄露程度;数据质量;个性化隐私保护【作者】宋海娜;罗涛;韩新宇;李剑峰【作者单位】北京邮电大学北京先进信息网络实验室,北京100876;北京邮电大学网络体系构建与融合北京市重点实验室,北京100876;北京邮电大学北京先进信息网络实验室,北京100876;北京邮电大学网络体系构建与融合北京市重点实验室,北京100876;北京邮电大学北京先进信息网络实验室,北京100876;北京邮电大学网络体系构建与融合北京市重点实验室,北京100876;北京邮电大学北京先进信息网络实验室,北京100876;北京邮电大学网络体系构建与融合北京市重点实验室,北京100876【正文语种】中文【中图分类】TP309.21 引言随机响应(Randomized Response,RR)是基于数据失真的本地化隐私保护技术的主流扰动机制，模型简洁直观且易于实现，并且其扰动程度可直接量化，在统计特性方面性能优良，因此受到广泛的关注[1～4].RR采用依概率作答的方式来保护受访者的隐私[4]，它保证敏感问题作答具有很强的可否认性，已经在Google Chrome的隐私保护工具和Apple系统中应用.同时，RR充分考虑了数据采集过程中数据收集者窃取或泄露用户隐私的可能性，该模型中受访者能够独立地对个体数据进行隐私化处理，大大激发其参与数据收集的积极性，并且不再需要可信第三方的介入，也免除了不可信第三方数据收集者可能带来的隐私泄露与攻击[5]. Warner于1965年首次提出利用RR技术[1](简称W-RR)进行敏感数据的收集，借鉴统计学研究中的经典方法，模拟调查者在尽量不侵犯被访者隐私的前提下收集到有价值的统计数据的过程.基于W-RR模型，文献[6～8]使用离散无记忆二进制对称信道对其进行建模，并有学者将其推广到多元离散对称隐私信源提出K-RR[9]、O-RR[10]等，得到了广泛的应用.差分隐私[11]作为一种公认的强健的隐私保护模型被提出，受到了很多学者的青睐[12～15].顺应这一发展趋势，很多学者对RR模型下的本地化差分隐私机制进行研究，通过单一隐私预算参数ε来衡量机制的隐私保护程度，主要集中在最优的差分隐私机制的设计与分析，重点关注机制客观的隐私-效用的折中[14,15].以上随机响应模型(简称CRR模型)及其相关研究均默认不同隐私数据同等重要，隐私保护需求相同，进而对所有的敏感值都执行相同的隐私保护操作，主要关注的是整个体系的客观隐私保护程度.然而，这可能会造成某些敏感度低的敏感值过度保护，而某些敏感度高的敏感值欠缺保护.这是因为，实际中隐私信息并非同等重要，隐私保护的需求也不尽相同.换言之，受访者对不同敏感值的敏感度的主观感受存在很大差异，因而隐私保护需求也就不同，此时隐私信息的敏感度在一定程度上反映了隐私需求的高低[16].实际中，不同敏感值的敏感度有很大差异，隐私保护需求也不同.如果硬性地对所有敏感值进行相同等级的隐私保护，直观上会降低统计数据的质量.另一方面，对高敏感度的敏感值进行足够高的隐私保护，而对低敏感度的敏感值进行相对较小的隐私保护(达到其隐私需求即可)，并将各个敏感值的主观隐私泄露程度直观地反馈给受访者，能够激发受访者的积极性，且激励受访者更加客观地提交扰动后的数据，进而使得收集到的数据更加客观，数据可用性更高.根据课题组的调研来看，针对RR机制的相关研究，截止目前还未找到解决上述问题类似的研究.因此，本文在CRR模型的基础上，针对多元离散隐私信源，考虑个性化的隐私保护需求，引入敏感值权重的概念，并将其引入到随机响应决策中，提出一种面向多敏感值的个性化随机响应(Personalized Randomized Response,PRR)机制.分析表明，相比于CRR机制，本文所提的PRR机制收集到的统计数据的质量会更高，同时实现了个性化的隐私保护，具有很强的实际意义.2 系统模型2.1 RR模型对于m元离散有限隐私信源X∈X(其中，采用扩展的m元RR模型来进行本地化隐私保护数据的收集，其RR过程可建模为多元离散信道的传输过程，其输出为如图1所示.隐私保护机制PY|X可描述为X=x映射到Y=y的条件概率PY|X(Y=y|X=x)=PY|X(y|x).若多元离散信道是对称的，则X=Y，本文主要考虑该种对称隐私信道.此时，当i=j时，xi=yj(i,j = 1,2,…,m).隐私分布估计问题描述为：隐私信源的先验分布为PX=[P1,P2,…,Pm]，参与者的样本量为N(N≫m)，样本X1,X2,…,Xn,…,XN(n = 1,2,…,N)是来自于分布PX，不同个体间隐私数据相互独立且同分布，且假设每个个体仅拥有一种敏感值.样本Xn 经过信道PY|X的随机扰动处理，得到扰动数据Yn，如图1所示.其中，Yn的分布是根据PY=PXPY|X来确定的，由Xn到Yn的过程称为RR过程，信道转移概率PY|X称为隐私保护机制.本文基于RR提出个性化隐私保护的目的是：在本地化数据收集中，使个体对自己的隐私可控，高隐私需求的个体执行较高强度的扰动，低隐私需求的个体执行较低强度的扰动，从而实现个性化的隐私保护，而且能够激发受访者参与的积极性，激励受访者客观地提供扰动数据，提高统计数据的质量.另外，在一定的隐私泄露条件下，数据收集者由YN={Y1,Y2,…,Yn,…,YN}和PY|X进行统计分析，以得到有效的估计结果.2.2 CRR模型CRR模型中默认隐私信息同等重要，隐私保护需求相同.设置CRR模型的隐私保护机制为PCRR，可由m×m的行随机矩阵表示：(1)即对每种敏感值都以P的概率保持原值，以(1-P)/(m-1)的概率响应输出其他m-1种敏感值中的一种，称之为m-CRR模型，其中PCRR为m×m的对称矩阵，且满足|PCRR|≠0，即PCRR可逆.由式(1)可知，CRR模型对所有敏感值都执行相同的隐私保护操作，这样会造成某些敏感度低的敏感值过度保护，而某些敏感度高的敏感值欠缺保护.然而，在很多情况下隐私信息并非同等重要，隐私保护需求也不同.基于此，本文研究面向敏感值的PRR机制.2.3 PRR模型PRR模型考虑不同敏感值的敏感度不同，隐私保护需求亦不同.因此，根据敏感度的不同为各个敏感值赋予一定的权重因子，即敏感值权重，以便对不同敏感值进行个性化的设计与分析.定义1 敏感值权重wi：对于敏感值xi，通过权重wi (i = 1,2,…,m)来刻画用户对不同敏感值的主观敏感度，即敏感值权重，也称为主观敏感度.注：隐私信息的敏感度与具体的应用场景、个体的主观感受与偏好、经济效益等诸多因素有关，因此本文重点不在敏感值权重的设置上，重点关注PRR机制本身的设计与分析.但是，敏感值权重wi的设定应遵循以下两个原则：(1)原则一，敏感值的敏感度越高，相应的敏感值权重越大.比如：在“AIDS”和“Flu”两种敏感值中，“AIDS”更加敏感，其隐私泄露所带来的影响或损失更大，许多受访者会采取不合作的态度，不予回答，这给抽样结果的分析造成极大的困难，因此设定“AIDS”的敏感值权重相对更大是合理的.(2)原则二，敏感值的频率越低，相应的敏感值权重越大.借鉴信息论中信息熵概念，可知：敏感值xi所占的比例越小，其信息量越多.对于敏感值xi，其主观敏感度越高，相应的敏感值权重wi越大，隐私保护的需求也越高.PRR模型考虑个性化的隐私保护需求，将wi引入到敏感值xi随机响应的决策中，其隐私保护机制PPRR为(2)称之为m-PRR模型.其中，Pi,i(wi)可形象地描述为拥有敏感值xi的个体讲真话的概率，且Pi,i(wi)是关于wi的减函数.这是因为，Pi,i(wi)从侧面反映了敏感值xi的隐私保护程度：Pi,i(wi)越小，敏感值xi的隐私保护程度越高，相应的隐私泄露风险越小，反之亦然.因此，敏感值xi的敏感度越高，保持原值的概率Pi,i(wi)应该设置的越小.此时，PRR模型能够激励用户更加客观地提交数据，进而使得收集的数据更加客观，特别是高敏感度的敏感值.3 基于RR机制的隐私泄露程度与统计数据质量分析3.1 数据隐私泄露程度分析3.1.1 CMAP准则下CRR机制的客观隐私泄露程度若仅已知隐私信息的先验分布PX，则此时正确推测出隐私信息的平均概率为(3)其中，PX(xi)表示先验分布中X=xi的概率，Pc01反映的是攻击者仅已知先验信息条件下的最大隐私泄露程度，其值越大，说明此时的隐私泄露程度越大.定义2 CMAP准则：选择重构函数使之满足条件(4)则称该准则为常规最大后验概率(Conventional Maximum A Posteriori,CMAP)准则”或“最小错误概率准则”.已知后验信息，CMAP准则下平均错误概率为(5)其中，dH(kj,i)表示汉明距离，Pe也称为非加权错误率或贝叶斯风险，能够客观地从侧面反映数据的安全性.其中，(6)同理，令CMAP准则下平均正确概率为Pc，由Pe+Pc=1，则Pc为Pc=1-Pe(7)其中，Pc是关于PX和PY|X的函数，也称为非加权(平均)正确率或贝叶斯效用，直观衡量了数据客观上的安全性：Pc越大，数据客观上的安全性越低.在一次猜测中，已知后验信息下正确猜测出个体真实信息的平均最大概率为Pc.定义3 客观隐私泄露程度：CRR机制的客观隐私泄露程度定义为已知后验信息后的正确率与仅已知先验信息条件下的正确率的比值[17](8)其中，Lc≥1表示CMAP准则下由CRR机制带来的相对隐私泄露程度：Lc越大，相应的隐私保护程度越低.对于敏感值xi的群体而言，其更加关注的是敏感值xi的隐私保护程度，而并非整个隐私保护体系的客观隐私保护程度，CMAP准则下正确推测出敏感值xi的平均概率为(9)其中，是对dH(kj,i)进行取反操作.Pci反映的是CMAP准则下原始敏感值为xi且能够正确推测出xi的平均概率.此时，所有敏感值平均正确率之和为机制的平均正确率，即定义4 敏感值xi的客观隐私泄露程度：针对敏感值xi，其客观隐私泄露程度Lci定义为CMAP准则下原始隐私信息为xi且能够正确推测出xi的平均概率与xi的先验概率的比值：(10)其中，Lci满足：0Lci1.当Lci=0时，说明CRR机制对敏感值xi起到了完全的保护.3.1.2 PMAP准则下PRR机制的主观隐私泄露程度考虑到不同隐私信息并非同等重要，隐私保护需求亦不同，因此PRR机制将充分考虑敏感值权重带来的影响.PRR机制下，若仅已知PX，根据加权隐私效益最大化原则，此时正确推测出隐私信息的平均加权概率为(11)其中，Pcw0反映的是仅已知PX下，根据隐私加权概率最大化原则，所获得的最大加权效益，其值越大，说明隐私泄露程度越大.定义5 PMAP准则：选择重构函数使之满足(12)其中，则称该准则为“隐私加权最大后验概率(Privacy-weighted Maximum A Posteriori,PMAP)准则”或“隐私加权最小错误概率准则”.定义6 隐私加权正确率/贝叶斯效用：通过权重wi(i=1,2,…,m)刻画隐私谋取者获取后验信息的条件下，采用PMAP准则，令(13)则能够正确推测出带主观感受隐私信息的加权概率为(14)称该概率Pcw为隐私加权正确率或隐私加权贝叶斯效用.Pcw能够衡量机制主观数据安全性：Pcw越大，机制主观上的安全性越低.定义7 主观隐私泄露程度：PRR机制的主观隐私泄露程度定义为已知后验信息条件下的隐私加权正确率与仅已知先验信息条件下的隐私加权正确率的比值：(15)其中，Lcw≥1反映的是PMAP准则下由PRR机制带来的相对隐私泄露程度，其值越大，相应的隐私泄露程度越高.同理，针对PRR机制，PMAP准则下敏感值为xi且正确推测出xi的隐私加权概率为(16)其中，Pcwi满足：定义8 敏感值xi的主观隐私泄露程度：PRR机制下，敏感值xi的主观隐私泄露程度Lcwi定义为PMAP准则下敏感值为xi且正确推测出xi的隐私加权概率与敏感值xi的隐私加权先验概率的比值：(17)其中，Lcwi满足：0Lcwi1.Lcwi越大，说明敏感值xi主观上的隐私泄露程度越高，意味着PRR机制对敏感值xi保护程度越低.为了兼容CRR模型，敏感值权重wi的设置应满足“归m化”条件，即且0<wi<m.CRR模型在PRR模型中被视为敏感值权重相等的一种情况，即对所有的敏感值xi，其敏感值权重均为wi≡1(i=1,2,…,m)，此时基于PRR模型的所有理论推导与分析均适用于CRR模型.因此，采用归m化处理后，CRR模型为PRR模型的一种特例.3.2 统计数据质量分析基于RR的本地化隐私保护机制旨在将隐私泄露控制在源端，同时关注统计数据的质量，本小节使用估计误差来评估两种模型的隐私分布估计问题[10].将敏感值权重按照从大到小的顺序排列：wi1≥wi2≥…≥wim，相应敏感值的隐私保护需求依次减小，由wi与PY|X的关系，相应经过RR模型输出真实值的概率应该满足：Pi1,i1Pi2,i2…Pim,im，其中敏感值权重wi1对应的敏感值xi1群体隐私保护需求最大.3.2.1 CRR模型下的隐私分布的估计误差分析CRR模型下，不同敏感值同等重要，隐私需求相同.此时，对各个敏感值均执行相同的隐私保护操作，设置Pi1,i1=Pi2,i2=…=Pim,im=P.敏感值xi的先验概率PX(xi)简记为Pi，其经验估计值为根据最大似然估计准则：(18)其中，表示CRR模型下输出数据YN={Y1,Y2,…,Yn,…,YN}中满足Yn=yi的数量.是其真实值Pi的无偏估计，该过程中Pi的估计误差为(19)CRR模型下PX的经验估计记为由可描述为(20)其中，为PY的经验估计，为PCRR的逆矩阵.根据谢尔曼-莫里森公式，对于任意m≥2，N≫m和P≠1/m，则为(21)其中，(·)i,j表示矩阵的第i行，第j列的元素.采用来评估CRR模型的隐私分布估计问题[10]，称为估计误差，则(22)其中，表示CRR模型下输出数据YN中yj所占的实际比例，Pyj是的期望值.由式(19)和式(21)，式(22)亦可描述为(23)越小，说明估计误差越小，意味着隐私分布估计的准确性越高.考虑到敏感值xi1的敏感值权重wi1最大，其隐私保护需求最高，则Pi1,i1=P最小.本文借鉴K-RR模型的参数设计，引入参数λ>0，为了满足所有敏感值的隐私需求，同时要求隐私泄露风险最低，不妨构造式(1)中PCRR为(24)即对所有敏感值xi均以P=eλ/wi1/(m-1+eλ/wi1)的概率响应输出其真实值，这是满足CRR模型隐私需求的参数设置中的一个特例.实际中，λ可由敏感值xi1具体的隐私保护需求进行设置.另外，Lc和Lci均与λ有关，实际中可通过调节λ的大小来实现相应的隐私保护程度.3.2.2 PRR模型下的隐私分布的估计误差分析PRR模型中考虑到不同敏感值的敏感度不同，造成其隐私保护需求不同，因此进行个性化的隐私保护.PRR模型下Pi的经验估计值记为的经验估计记为(25)其中，为PPRR的逆矩阵.同理，PRR模型下隐私分布的估计误差为(26)其中，表示PRR模型下输出数据YN中yj所占的实际比例，Pyj表示的期望值. 同理，对于敏感值xi，根据wi与PPRR的关系，不妨构造式(2)中的PPRR为(27)即敏感值xi群体以eμi/wi/(m-1+eμi/wi)的概率响应其真实值，同时需满足：Pi1,i1Pi2,i2…Pim,im.其中，μi可由敏感值xi群体的具体隐私需求进行设置，实现了面向多敏感值的个性化隐私保护.当所有敏感值的敏感度均相同时，即对于任意i,j∈{1,2,…,m},μi=μj=λ，wi=wj=1时，式(27)与CRR模型的表达形式相同，因此式(27)是一种更为广泛的定义形式.对于敏感值xi而言，其wi越大，意味着其隐私保护需求越高，倾向响应输出其真实值的概率相对较低，而式(27)的约束规则刚好符合这一趋势.PRR模型的隐私分布估计值相对于CRR模型的隐私分布估计值的效率定义为相对效率Re[1]，描述为(28)在一定的主观隐私泄露程度下，最大化数据的统计质量，满足：Re≥1.这是因为CRR模型对所有的敏感值都采用相同的隐私保护操作，为了满足各自的隐私需求同时最小化隐私泄露程度，一般采用隐私保护需求最大的敏感值所需要的隐私保护参数设置，如此会造成其他敏感度低的敏感值过度保护，而PRR模型是根据具体敏感值的个性化隐私需求设置其隐私参数，实现了个性化的隐私保护.另一方面，低敏感度的敏感值群体对隐私保护的需求较低，倾向响应输出其真实值的概率相对较大，如此为统计数据的可用性做出较大贡献.显然，本文PRR机制下的估计值相比于CRR模型更加有效.4 实例分析为了验证本文所提PRR方案的有效性和可行性，本小节通过MATLAB仿真模拟调查者使用图1的RR模型来进行隐私分布估计的过程.设置参与样本数N=105，样本取值集合X={x1,x2,x3}={AIDS,Cancer,Flu}，即m=3，每种疾病代表一种敏感值，样本来源于分布PX=[P1,P2,P3]=[0.2,0.3,0.5](注：实际中PX是未知的，且本小节设置m=3，PX=[0.2,0.3,0.5]只是一个特例，本文PRR机制适用于任意m≥2元离散有限敏感值下基于RR的个性化隐私保护).其中，根据wi的设置原则，三种敏感值权重如表1所示.可知，“AIDS”的敏感度最高，其敏感值权重最大；“Flu”的敏感度最低，其敏感值权重最小.同时，3-CRR和3-PRR模型分别采用式(24)和(27)的隐私保护机制.为便于分析，PRR模型中设置μi=μj=μ(i,j=1,2,3).表1 三种敏感值权重及其先验概率分布敏感值Xx1x2x3疾病类型AIDSCancerFlu 敏感值权重wi(×3)0.60.30.1概率0.20.30.5图2分别展示了3-CRR模型在CMAP准则下机制的客观隐私泄露程度Lc和3-PRR模型在PMAP准则下机制的主观隐私泄露程度Lcw.可知，CRR(或PRR)模型下机制的客观(或主观)隐私泄露程度随着参数λ(或μ)的增大而增大，意味着机制的隐私保护程度越来越低，即隐私保护程度与λ(或μ)呈负相关的关系，这与3.1小节的分析一致.在实际应用中，应该综合考虑隐私保护和数据质量的需求来设置合适的λ和μ值，进而将机制的隐私泄露程度控制在一定的安全范围，同时又能保证一定的数据可用性.图3给出了3-CRR模型在CMAP准则下各个敏感值的客观隐私泄露程度Lci.可知，当λ≤1.65时(该临界值与具体的先验分布有关)，隐私泄露程度由底到高依次为x1、x2、x3，表明敏感度最高的敏感值x1的隐私保护程度最高，基本对其实现了完全的保护.但是，当λ>1.65时，三条曲线基本重合，说明此时CRR机制对三种敏感值客观上的隐私保护程度相同.实际中，各个敏感值的隐私需求并不同，但CRR模型对所有敏感值的隐私保护程度相同，这是其不足之处.图4给出了3-PRR模型在PMAP准则下各个敏感值的主观隐私泄露程度Lcwi.可知，三个敏感值的主观隐私泄露程度均随μ的增大而增大.当μ一定时，x1、x2、x3的主观隐私保护泄露程度依次增大，意味着相应的主观隐私保护程度依次增大，这也是本文PRR模型相比于CRR模型的优势所在.换言之，PRR机制对敏感度最高的敏感值x1的隐私保护程度最高，对敏感度最低的敏感值x3的隐私保护程度相对较低，这是因为其相应的隐私保护需求较低.为了验证所提PRR机制在统计质量方面优于CRR机制，本小节设置两种模型下的主观隐私泄露程度一定.图5给出3-CRR和3-PRR模型隐私分布估计误差的仿真值，并给出其理论值作为参考(蒙特卡洛仿真设置104次).可知，两种模型下隐私分布估计误差的仿真值与理论值基本一致，验证了本文理论分析的正确性.另外，当主观隐私泄露程度Lcw一定时，CRR模型隐私分布的估计误差均高于PRR模型，即PRR模型统计估计的准确性更高，性能更佳.另一方面，两种模型下的估计误差均随着主观隐私泄露程度Lcw的增大而减小.但是，实际中，受访者一般希望其主观隐私泄露程度越小越好，但如此会造成数据可用性降低，说明数据可用性与隐私之间是相互矛盾的.实际中，应该根据数据质量和隐私需求进行良好的折中.表2展示了五种主观隐私泄露程度下两种模型各自因此参数设置以及隐私分布估计的相对效率Re的仿真值，并给出理论值作为参考(蒙特卡洛仿真设置104次).可知，当两种模型的主观隐私泄露程度相同时，Re的理论值和仿真值基本一致，验证了本文理论分析的正确性和有效性(与图5结论一致).另外，Re均满足：Re≥1，意味着主观隐私泄露程度一定时，PRR模型在统计数据可用性方面优于CRR模型，并且能够针对不同敏感值灵活调控参数μ，实现个性化的隐私保护，这与3.2小节的理论分析一致.再者，两种模型的主观隐私泄露程度均随着λ和μ的增大而增大，此时隐私数据的安全性在不断降低.因此，实际中需要根据具体隐私保护需求设置合适的λ和μ值，进而将隐私泄露程度控制在合理的范围.表2 主观隐私泄露程度一定的条件下两种模型的隐私参数设置及其相对效率Re的仿真值与理论值主观隐私泄露程度LcwCRR参数:λPRR参数:μ相对效率Re理论值仿真值1.20001.63690.76881.16121.18211.40002.33161.23671.41431.30021.60003. 11611.86851.44571.44401.80004.11172.75361.31971.32012.00005.72054.3 1901.16121.14705 小结与展望实际中，隐私信息并非同等重要，隐私需求也不尽相同.而现有的CRR模型对所有的敏感值都执行相同的隐私保护操作，如此会造成某些低敏感信息过度保护，而某些高敏感信息欠缺保护.基于此，本文考虑个性化的隐私需求，引入敏感值权重的概念，提出了一种面向多敏感值的PRR机制，实现了个性化隐私保护，具有重要的实际意义.本文针对CRR和PRR两种机制下的隐私泄露程度和数据质量进行了充分的阐述与分析.分析结果表明，在主观隐私泄露程度一定时，相比于CRR机制，PRR机制所获得的统计估计的准确性较高.另外，为了便于讨论，本文仿真将PRR 模型的隐私参数设置为μi = μj =μ，这样只能对各敏感值的隐私保护程度进行粗粒度的调节.因此，接下来的研究可对各敏感值设置合理的μi值，旨在对各敏感值进行灵活且细粒度的隐私保护调节，实现个性化精准的隐私保护，并在一定主观隐私泄露程度的限制下，通过优化设计，最大化统计数据的可用性.参考文献【相关文献】[1] Warner S L.Randomized response:A survey technique for eliminating evasive answer bias[J].Journal of the American Statistical Association,1965,60(309):63-69.[2] 罗永龙,黄刘生,荆巍巍,等.一个保护私有信息的布尔关联规则挖掘算法[J].电子学报,2005,33(5):133-136.Luo Y L,Huang L S,Jing W W,et al.An algorithm for privacy-preserving boolean association rule mining[J].Acta Electronica Sinica,2005,33(5):133-136.(in Chinese)[3] Hsieh S H,Lee S M,Tu S H.Randomized response techniques for a multi-level attribute using a single sensitive question[J].Statistical Papers,2018,59(1):291-306.[4] Tian X,Taylor J.Selective inference with a randomized response[J].The Annals of Statistics,2018,46(2):679-710.[5] 叶青青,孟小峰,朱敏杰,等.本地化差分隐私研究综述[J].软件学报,2018,29(7):159-183.Ye Q Q,Meng X F,Zhu M J,et al.Survey on local differential privacy[J].Journal of Software,2018,29(7):159-183.(in Chinese)[6] Lin B C,Wu S H,Tsou Y T,et al.PPDCA:Privacy-preserving crowdsensing data collection。

随机振动响应分析技术研究一、引言随机振动响应分析是结构工程领域中一个非常重要的课题。

结构物的振动响应具有随机性、复杂性和非线性等特点，因此，能够对结构物在随机激励下的振动响应进行研究和分析，对于提高结构物的可靠性、耐久性和安全性非常关键。

二、随机振动响应分析的方法随机振动响应分析技术主要包括两种方法：频域分析和时域分析。

1. 频域分析频域分析是指将随机振动信号分解成一系列特定频率的正弦波分量，然后对这些正弦波分量进行分析、计算和处理。

这种方法一般使用离散傅里叶变换（DFT）或快速傅里叶变换（FFT）进行处理，可以方便地进行频率分析和频率响应。

2. 时域分析时域分析是指基于时间序列的方法，通过对随机振动信号的时间序列进行分析，得到结构物的响应特性。

这种方法可以使用自相关函数、互相关函数、功率谱密度和相干函数等分析工具。

三、随机振动响应分析的应用随机振动响应分析技术在各个领域都有广泛的应用。

1. 土木工程在土木工程中，随机振动响应分析技术可以用来评估建筑物、桥梁、隧道等结构物在地震或风荷载下的响应情况，以及评估疲劳损伤的程度。

2. 航空航天工程在航空航天工程中，随机振动响应分析技术可以用来评估航天器在发射过程中的响应情况，以及评估机体结构在飞行过程中的疲劳损伤程度。

3. 机械工程在机械工程中，随机振动响应分析技术可以用来评估机械系统在振动环境下的可靠性和安全性，以及寻找和消除机械系统的振动问题。

四、随机振动响应分析技术的发展趋势随着科学技术和计算机技术的快速发展，随机振动响应分析技术也得到了极大发展和应用。

未来，随机振动响应分析技术的发展主要将呈现以下几个趋势：1. 多物理场耦合建模针对涉及多种物理场同时作用的振动问题，将机械、声学、热学、流体力学等多种物理场有机结合起来，建立更加全面且真实的多物理场耦合模型，以便更好地分析和解决复杂振动问题。

2. 精细化建模分析建立尽可能精细的结构物和振动环境的建模，以更加准确地反映实际情况，预测结构物的振动响应和疲劳损伤情况，从而提高结构物的可靠性和安全性。

收稿日期:2003-10-26;　修回日期:2003-11-22 基金项目:国家电力公司资助项目(KJ 00-03-26-01) 作者简介:朱以文(1945-),男,教授,主要从事计算力学和结构防灾减灾研究文章编号:1000-1301(2003)06-0174-05TM D 多点控制体系随机地震响应分析的虚拟激励法朱以文,吴春秋(武汉大学土木建筑工程学院,湖北武汉430072)摘要:对于频率分布密集或受频带较宽的地震激励的结构,其响应不再以某一单一振型为主,须考虑采用多点控制。

本文对受T M D 多点控制的结构进行了研究。

文中建立了带有多个子结构系统的以模态坐标和子结构自由度为未知量的统一运动方程。

针对所得方程为非对称质量、非对称刚度、非经典阻尼的情况,本文给出了使用直接法求解的格式。

地震随机响应分析采用了虚拟激励法,可以考虑各振型之间的耦合项,计算量小且精度高。

本文的方法适用于带有多个子结构的系统的一般性问题,具有广泛的应用价值。

关键词:多点控制;主结构;子结构;随机地震响应中图分类号:P315.96 文献标识码:A Pseudo -excitation method for random earthquake response analysis of control system with MTMDZH U Yi -wen ,WU Chun -qiu(Civil and structural engineering school ,W uhan university ,Wuhan 430072,China )A bstract :The response of the structure is no t constituted with one sing le mode shape w hen the frequency distri -bution is dense o r the earthquake excitation 's frequency band is w ide .At this time ,it is necessary to adopt the multi -point control sy stem .The study on the structures w ith M TMD is carried out in this paper .The uniform dynamic equation w ith mode coordinate and slave system 's DOF as variables is established fo r the system w ith multi slave sy stem .The equatio n has asy mmetric mass m atrix ,asymmetric stiffness matrix and nonclassicaldamping m atrix ,and the direct solving format is given in this paper .The random earthquake response is studied by using pseudo -excitation method ,thus the coupling items between modes can be considered .The calculation is cheap and precision is high .The method in this paper is adaptable to the general case of the sy stem with multi -slave structures and has broad application wo rth .Key words :multi -point control ;master structure ;slave structure ;random earthquake response1　引言 对于高层建筑、大跨桥梁、高耸塔架等高柔结构采用TMD (Tuned Mass Damper )减小风振及地震响应是有效的,这一点得到了人们的普遍认同。

简支厚板的功率谱密度随机及简谐响应分析/PREP7 !进入前处理模块/TITLE, EX 8.4(4) by Zeng P, Lei L P, Fang GET,1,SHELL93,,,,,1 !设定1号单元(shell93)L=10 $t=1 !设定几何参数MP,EX,1,2E11 $MP,NUXY,1,0.3 $MP,DENS,1,8000 !设定弹性模量, 泊松比, 密度R,1,t !设定实常数N,1,0,0,0 $N,9,0,L,0 ! 生成节点FILL !填充生成节点NGEN,5,40,1,9,1,L/4 !复制生成节点N,21,1.25,0,0 $N,29,1.25,L,0 ! 定义节点FILL,21,29,3 !填充生成节点NGEN,4,40,21,29,2,L/4 !复制生成节点EN,1,1,41,43,3,21,42,23,2 !由连接节点来生成指定单元编号的单元EGEN,4,2,1 !复制生成单元EGEN,4,40,1,4 !复制生成单元D,ALL,UX,0,,,,UY,ROTZ !对所有节点施加约束UX=UY=ROTZ=0D,1,UZ,0,0,9,1,ROTX !对1-9节点施加约束UZ=ROTX=0D,161,UZ,0,0,169,1,ROTX !对161-169节点施加约束UZ=ROTX=0D,1,UZ,0,0,161,20,ROTY !对1-161(以20为增量)节点施加约束UZ=ROTY=0D,9,UZ,0,0,169,20,ROTY !对9-169(以20为增量)节点施加约束UZ=ROTY=0NSEL,S,LOC,X,.1,9.9 !选择位置0.1<x<9.9之间的节点NSEL,R,LOC,Y,.1,9.9 !再选择位置0.1<y<9.9之间的节点M,ALL,UZ !对所选择的节点设置主自由度UZNSEL,ALL !选择所有节点FINISH !结束前处理模块/SOLU !进入求解模块ANTYPE,MODAL !设置模态分析方式MODOPT,REDUC !设置缩减算法MXPAND,16,,,YES ! 设定模态扩展的阶数为16SFE,ALL,,PRES,,-1E6 ! 在单元上施加表面载荷SOLVE !进行模态求解*GET,F1,MODE,1,FREQ !提取1阶模态的频率值，赋给参数F1FINISH !退出模态分析/SOLU !进入求解模块ANTYPE,SPECTR !设定谱分析方式SPOPT,PSD,2,YES !设定2阶模态的PSD响应，并计算应力PSDUNIT,1,PRES !设定PSD谱分析为压力谱DMPRAT,0.02 !设定阻尼比为0.02PSDFRQ,1,1,1.0,80.0 ! 输入PSD数据表的频率值，前两位数值为数据表的编号PSDV AL,1,1.0,1.0 ! 输入PSD响应值，前一位数为数据表的编号LVSCALE,1 ! 模态叠加时的载荷乘子PFACT,1,NODE !计算针对节点激振的参与系数PSDRES,DISP,REL !设定针对位移的输出PSDCOM !设定功率谱密度的模态合并方法SOLVE $FINISH !进行求解，结束求解模块/POST1 !进入一般性后处理模块SET,3,1 !调入第3模态第1子步的结果/VIEW,1,2,3,4 !设置视图的显示方向PLNSOL,U,Z !图形显示Z方向的位移云图P1=NODE(L/2,L/2,0) !获取中间位置处的节点编号*GET,P_EQV,NODE,P1,S,EQV !提取节点P1的等效应力值，赋予P_EQV NSEL,,NODE,,P1 !选择节点P1PRNSOL,S,COMP !打印节点P1的应力值NSEL,ALL !选择所有节点FINISH !退出后处理/SOLU !进入求解模块ANTYPE,HARMIC ! 设定简谐分析方式HROPT,MSUP ! 设定模态叠加算法HROUT,OFF,ON !设置输出结果为幅值和相角方式KBC,1 !设置阶梯式加载方式HARFRQ,1,80 !定义简谐响应中的频率范围为1~80 DMPRAT,0.02 !设定阻尼比为0.02NSUBSTEP,10 !设定子步数SOLVE !进行求解FINISH !退出/POST26 !进入时间历程后处理FILE,,rfrq !设定输出结果文件，以便从简谐分析中获得位移结果PRCPLX,1 !设置幅值和位相角的输出方式NSOL,2,P1,U,Z !将节点P1(板的中心位置)的位移UZ设置为2号变量PSDDAT,6,1,1.0,80,1.0 !设置PSD计算的频率和PSD值，赋给6号变量PSDTYP,2PSDCAL,7,2 !对2号变量进行PSD计算，再设定为7号变量PSDPRT !打印变量值PRV AR,2,7 !列出2号及7号变量随频率变化的数值*GET,P1_UZ,V ARI,7,EXTREM,VMAX !提取7号变量中的最大值，赋给P1_UZ P1_UZ= P1_UZ*1000000 !将m^2单位转换为mm^2*STATUS !列显所有参数内容/VIEW/AXLAB,Y,PSD (M^2/HZ) !定义Y轴标签为PSD (M^2/HZ)PLV AR,7 !图形给出7号变量的曲线FINISH。

随机振动分析随机振动是只能从统计的角度描述的振动。

在任何给定的时间内，瞬态幅值都是未知的，它们用其统计特性（如平均值、标准方差和超过某个值的可能性）来表示。

随机振动的示例包括地震运动、海浪的高度和频率、飞机和高层建筑上的风压波动以及因火箭和喷气式发动机噪声引起的声音激励。

这些随机的激励通常用功率频谱密度 (PSD) 函数来描述。

NX Nastran 在频率响应分析之后的后处理步骤中执行随机响应分析。

频率响应分析用于生成传递函数（即输出输入比）。

将输入 PSD 乘以传递函数可得到响应PSD。

输入 PSD 可采用自动谱密度或交叉谱密度的形式。

随机响应输出由如下值组成：响应 PSD、ATOC（自相关函数）、每单位时间中具有正斜率的零交叉的数量以及响应的 RMS（均方根）值和 CRMS（累积均方根）。

可使用参数 RMSSF 对 RMS 和 CRMS 进行按比例调整。

默认情况下，频率响应输出会在随机震动分析中被抑制。

要获取频率响应输出，请指定 SYSTEM(524)=1。

NX随机响应计算用来定义随机过程的函数功率谱密度 (PSD) 函数功率谱密度 (PSD) 函数又称作自动谱密度函数，软件使用该工具来定义和计算随机过程（激励或响应）。

PSD 函数是实数函数，它用均方值的单面光谱密度来描述随机激励 (1)，其定义如下所示：方程 1其中 ( )* 是复共轭的转置矩阵多个激励之间的关联使用创建关联命令，可以将多个随机激励关联在一起。

关于更多信息，请参见PSD 相关。

PSD 相关随机事件是基于统计上的分析，这种分析允许您应用一次或多次功率谱密度(PSD) 函数激励。

PSD 激励可以表示对您并不知道其准确大小的力的取样。

默认情况下，PSD 函数是独立的（不相关）。

但是，新建相关性对话框允许您使用相位角或时间延迟将两个 PSD 激励相关。

注释您仅可以在相同类型激励之间定义相关性。

例如，可在某分布式载荷和另一分布式载荷之间定义相关性，在节点力和节点力之间或在强迫运动和强迫运动之间定义相关性。

• 37•纤维增强复合材料在结构上具有多尺度特性与空间随机性，其尺度结构、组份材料性能参数均会影响到材料的力学性能。

本文建立了一种基于PCE与Vine Copula方法的多尺度随机力学性能预测方法，能够为CFRP材料的力学性能预测与受力、变形状态评估提供参考价值。

1 材料特性与方法选择1.1 碳纤维复合材料碳纤维复合材料又称为碳纤维增强聚合物基复合材料（CFRP），是一种密度低、比模数大、比强度高的轻质复合材料，具备良好的力学性能，在当前电子产品轻量化趋势下被广泛应用于微型电路芯片、锂电池电极等电子产品的制造生产领域。

CFRP材料因其制备工艺、存储条件、组成相成分等均具有不确定性特征，这种特征反映在材料性质上主要体现为多尺度力学性能的随机性，最终将作用于材料的随机性能，因此本文拟针对CFRP材料的随机力学性能进行测定，并分析影响材料宏观力学性能预测结果的主要因素。

1.2 多尺度分析方法当前国内外学者在针对复合材料随机力学性能预测的研究方面取得了一系列进展：一方面从研究纤维束的尺度入手，现有研究成果主要通过调节纤维的角度、位移等参数，通过改变其约束条件生成所需的材料结构。

例如有学者建立了一种序列随机扰动算法，结合有限元分析方法判断改变纤维的随机分布结构后，纤维束的力学性能将发生哪些变化；有学者采用随机序列展开方法，以介观尺度作为研究切入点，运用图像分析方法与数学统计学方法建立具有随机性RVE结构，并利用仿真软件实现对结构特征的直观分析；有学者针对影响材料结构排列特征的参数进行相关性分析，运用混合高斯随机序列进行算法重构，重新生成符合随机性特征的RVE模型。

另一方面以解析细观力学方法作为切入点，结合计算细观力学存在的计算代价高等缺陷，将解析细观力学方法运用在不确定性预测研究领域，用于提高计算效率。

例如有学者选取复合材料层合板作为研究对象，利用多项式与函数进行材料随机自由振动分析，并运用随机有限元方法进行该材料微观结构的预测；有学者运用Copula函数表示出材料参数对于时复合材料结构、性能的影响，采用摄动法进行材料微观结构的不确定性分析；有学者提出基于PCE 的层级传递方法，针对材料微观结构的分布形态进行分析，进而实现对宏观材料力学性能的预测。

第11章随机响应分析11.1 动力学环境分类11.2 概述1）随机振动是统计意义下描述的振动，在任何瞬时大小未知，但其大小的概率超过一给定的值。

2）常见的例子如地震引起的地基运动、海洋波浪高度和频率、航天器和高耸建筑物受到的风压力、由于火箭与喷气发动机噪音引起的声波等。

3）MSC/NASTRAN 对随机响应分析是作为频率响应后处理进行的。

输入包括频率响应的输出、用户给定的载荷条件（形式为自相关的谱密度）。

输出为响应功率谱密度、自相关函数、响应的均方值。

4）MSC/NASTRAN 随机分析假设历经性随机过程5）随机动态环境例子11.3 自相关与自谱1)自相关函数注：R j(0)为均方值2)自谱函数Fourier变换为3）均方响应值4）外观频率为N06）例子11.4 各态历经性随机激励下线性系统响应计算1) 线性系统单输入输出关系由频率响应分析得到其中，H ja(ω)为频率响应或输入到输出的传递函数对多输入单输出其矩阵形式为输出自相关谱为其单个输入谱为2）线性系统的多输入输出关系多输入输出谱关系其中，输入互谱矩阵为其谱特性为3）常用特殊情况（1）单输入分析（完全相关输入）（2）不相关多输入11.5 MSC/NASTRAN中随机分析的实现1)如果由频率响应计算结果为H ja(ω),但并不直接计算2)如需要H ja(ω),令F a(ω)=111.5.1RANDPS卡片1)定义随机分析中使用的功率谱密度因子，频率相关形式为2) 格式3) 由情况控制卡RANDOM = SID选取4)自谱密度，J=K, X为大于0的整数，Y为05)TID=0, G(F)=011.5.2 TABRND1卡片1)用表格函数定义功率谱密度函数2)格式3）11.5.3 随机响应输入要求1）执行控制2）情况控制3）模型数据11.5.4 随机响应例子1）例1：单输入随机响应分析（1）问题描述：(a)矩形板如图；(b)基座运动（z方向）功率谱(PSD)表中给出；(c)整个频率范围的常临界阻尼比为0.03；(d)用log-log输入PSD;(e)使用模态求解法(2) 使用具有大质量的模态法（在边界处用REB2单元）确定a)9999点处的位移和加速度功率谱(PSD)b)确定结点33和55的位移功率谱(PSD)(3)输入文件ID SEMINAR, PROB10SOL 111TIME 30CENDTITLE= RANDOM ANALYSIS - BASE EXCITATIONSUBTITLE= USING THE MODAL METHOD WITH LANCZOSECHO= UNSORTEDSPC= 101SET 111= 33, 55, 9999ACCELERATION(SORT2, PHASE)= 111METHOD= 100FREQUENCY= 100SDAMPING= 100RANDOM= 100DLOAD= 100$OUTPUT(XYPLOT)XTGRID= YESYTGRID= YESXBGRID= YESYBGRID= YESYTLOG= YESXTITLE= FREQUENCYYTTITLE= ACCEL RESPONSE BASE, MAGNITUDEYBTITLE= ACCEL RESPONSE AT BASE, PHASEXYPLOT ACCEL RESPONSE / 9999 (T3RM, T3IP)YTTITLE= ACCEL RESPONSE AT TIP CENTER, MAGNITUDEYBTITLE= ACCEL RESPONSE AT TIP CENTER, PHASEXYPLOT ACCEL RESPONSE / 33 (T3RM, T3IP)YTTITLE= ACCEL RESPONSE AT OPPOSITE CORNER, MAGNITUDEYBTITLE= ACCEL RESPONSE AT OPPOSETE CORNER, PHASEXYPLOT ACCEL RESPONSE / 55 (T3RM, T3IP)$$ PLOT OUTPUT IS ONLY MEANS OF VIEWING PSD DATA$XGRID= YESYGRID= YESXLOG= YESYLOG= YESYTITLE= ACCEL P S D AT LOADED CORNER XYPLOT ACCEL PSDF / 9999(T3)YTITLE= ACCEL P S D AT TIP CENTERXYPLOT ACCEL PSDF / 33(T3)YTITLE= ACCEL P S D AT OPPOSITE CORNER XYPLOT ACCEL PSDF / 55(T3)$BEGIN BULKPARAM,COUPMASS,1PARAM,WTMASS,0.00259$INCLUDE ’plate.bdf’$GRID, 9999, , 0., 0., 0.$RBE2, 101, 9999, 12345, 1, 12, 23, 34, 45 $SPC1, 101, 12456, 9999$CONM2, 6000, 9999, , 1.0E8$$MAT1, 1, .1, , .1, .286$$ EIGENVALUE EXTRACTION PARAMETERS$EIGRL, 100 , , 2000.$$ SPECIFY MODAL DAMPING$TABDMP1, 100, CRIT,+, 0., .03, 10., .03, ENDT$$ POINT LOADING AT TIP CENTER$RLOAD2, 100, 600, , , 310$TABLED1, 310,+, 10., 1., 1000., 1., END T$DAREA, 600, 9999, 3, 1.E8$$ SPECIFY FREQUENCY STEPS$FREQ,100,30.FREQ1,100,20.,20.,50FREQ4,100,20.,1000.,.03,5$$ SPECIFY SPECTRAL DENSITY$RANDPS, 100, 1, 1, 1., 0., 111$TABRND1, 111,LOG,LOG+, 20., 0.1, 30., 1., 100., 1., 500., .1,+, 1000., .1, ENDT$ENDDATA（4）部分结果基座的加速度PSD(大小、相位)与频率关系结点55的加速度PSD(大小、相位)与频率关系结点33的加速度PSD(大小、相位)与频率关系例2：多输入随机响应分析问题：输入文件ID SEMINAR, PROB11SOL 111TIME 30CENDTITLE= FREQUENCY RESPONSE WITH PRESSURE AND POINT LOADSSUBTITLE= USING THE MODAL METHOD WITH LANCZOSECHO= UNSORTEDSPC= 1SET 111= 11, 33, 55DISPLACEMENT(PLOT, PHASE)= 111METHOD= 100FREQUENCY= 100SDAMPING= 100RANDOM= 100SUBCASE 1LABEL= PRESSURE LOADDLOAD= 100LOADSET= 100SUBCASE 2LABEL CORNER LOADDLOAD= 200LOADSET= 100$OUTPUT (XYPLOT)$XTGRID= YESYTGRID= YESXBGRID= YESYBGRID= YESYTLOG= YESYBLOG= NOXTITLE= FREQUENCY (HZ)YTTITLE= DISPLACEMENT RESPONSE AT LOADED CORNER, MAGNITUDE YBTITLE= DISPLACEMENT RESPONSE AT LOADED CORNER, PHASE XYPLOT DISP RESPONSE / 11 (T3RM, T3IP)YTTITLE= DISPLACEMENT RESPONSE AT TIP CENTER, MAGNITUDE YBTITLE= DISPLACEMENT RESPONSE AT TIP CENTER, PHASEXYPLOT DISP RESPONSE / 33 (T3RM, T3IP)YTTITLE= DISPLACEMENT RESPONSE AT OPPOSITE CORNER, MAGNITUDE YBTITLE= DISPLACEMENT RESPONSE AT OPPOSITE CORNER, PHASE XYPLOT DISP RESPONSE / 55 (T3RM, T3IP)$$ PLOT OUTPUT IS ONLY MEANS OF VIEWING PSD DATA$XGRID= YESYGRID= YESXLOG= YESYLOG= YESYTITLE= DISP P S D AT LOADED CORNERXYPLOT DISP PSDF / 11(T3)YTITLE= DISP P S D AT TIP CENTERXYPLOT DISP PSDF / 33(T3)YTITLE= DISP P S D AT OPPOSITE CORNERXYPLOT DISP PSDF / 55(T3)$BEGIN BULKPARAM,COUPMASS,1PARAM,WTMASS,0.00259$$ MODEL DESCRIBED IN NORMAL MODES EXAMPLE$INCLUDE ’plate.bdf’$$ EIGENVALUE EXTRACTION PARAMETERS $EIGRL, 100, 10., 2000.$$ SPECIFY MODAL DAMPING$TABDMP1, 100, CRIT,+, 0., .03, 10., .03, ENDT$$ FIRST LOADING$RLOAD2, 100, 300, , , 310$TABLED1, 310,+, 10., 1., 1000., 1., ENDT$$ UNIT PRESSURE LOAD TO PLATE$LSEQ, 100, 300, 400$PLOAD2, 400, 1., 1, THRU, 40$$ SECOND LOADING$RLOAD2, 200, 600, , , 310$$ POINT LOAD AT TIP CENTER$DAREA, 600, 11, 3, 1.$$ SPECIFY FREQUENCY STEPS$FREQ1, 100, 20., 20., 49$$ SPECIFY SPECTRAL DENSITY$RANDPS, 100, 1, 1, 1., 0., 100 RANDPS, 100, 2, 2, 1., 0., 200 RANDPS, 100, 1, 2, 1., 0., 300 RANDPS, 100, 1, 2, 0., 1., 400$TABRND1, 100,+, 20., 0.1, 30., 1., 100., 1., 500., .1,+, 1000., .1, ENDT$TABRND1, 200,+, 20., 0.5, 30., 2.5, 500., 2.5, 1000., 0.,+, ENDT$TABRND1, 300,+, 20., -.099619, 100., -.498097, 500., .070711, 1000., 0., +, ENDT$TABRND1, 400,+, 20., .0078158, 100., .0435791, 500., -.70711, 1000., 0., +, ENDT$ENDDATA。