第２０卷第１期
２０１１年３月
计算机辅助工程
ＣｏｍｐｕｔｅｒＡｉｄｅｄＥｎｇｉｎｅｅｒｉｎｇ
Ｖ０１．２０Ｎｏ．１Ｍａｒ．２０１１
文章编号：１００６—０８７１（２０１１）０１．０１０１．０５
用于不确定性分析的高斯过程响应面模型
设计点选择方法
刘信恩，
肖世富，莫军
（中国工程物理研究院总体工程研究所，四川绵阳６２１９００）
摘要：为推动高斯过程响应面模型在复杂耗时数值模拟不确定性分析中的应用，提出一种可自动实现位置优化的高效设计点选择方法，即先在标准超立方体上生成拉丁超立方设计点，然后利用输入变量的已知概率分布将其映射回原始设计空间．将该方法与基于假设均匀分布的传统拉丁超立方设计方法进行比较，探讨它们对所建立的高斯过程响应面模型和不确定性分析结果的影响．算例表明新方法具有一定优势．
关键词：不确定性分析；蒙特卡罗方法；高斯过程；响应面模型中图分类号：０２１２．８；０２１１．９
文献标志码：Ａ
ＤｅｓｉｇｎｐｏｉｎｔｃｈｏｉｃｅｏｆＧａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌａｐｐｌｉｅｄｔｏｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓ
ＬＩＵＸｉｎ’ｅｎ，ＸＩＡＯＳｈｉｆｕ，ＭＯＪｕｎ
（Ｉｎｓｔｉｔｕｔｅ
ｏｆＳｙｓｔｅｍｓＥｎｇｉｎｅｅｒｉｎｇ，ＣｈｉｎａＡｃａｄｅｍｙｏｆ
Ｅｎｇｉｎｅｅｒｉｎｇ
Ｐｈｙｓｉｃｓ，Ｍｉａｎｙａｎｇ
６２１９００，Ｓｉｃｈｕａｎ，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ａｎｅｆｆｉｃｉｅｎｔｍｅｔｈｏｄｏｆｄｅｓｉｇｎｐｏｉｎｔｃｈｏｉｃｅｉｓｐｒｅｓｅｎｔｅｄ，ｗｈｉｃｈｌｏｃａｔｉｏｎ
ｏｆｔｈｅ
ｐｏｉｎｔｓ，ｔｏ
ａｄｖａｎｃｅｃｏｍｐｌｅｘ
ｔｈｅ
ｃａｎ
ａｕｔｏｍａｔｉｃａｌｌｙｏｐｔｉｍｉｚｅｔｈｅ
ｓｕｒｆａｃｅ
ｍｏｄｅｌ
ｔｏ
ａｐｐｌｉｃａｔｉｏｎｏｆＧａｕｓｓｉａｎｐｒｏｃｅｓｓ
ｒｅｓｐｏｎｓｅ
ｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓｏｆｔｈｅｄｅｓｉｇｎｐｏｉｎｔｓｉｎ
ａ
ｔｉｍｅ—ｃｏｎｓｕｍｉｎｇｎｕｍｅｒｉｃａｌ
ｓｉｍｕｌａｔｉｏｎ．Ｔｈｅｍｅｔｈｏｄ
ｐｒｏｄｕｃｅｓ
ｔｈｅ
ｓｔａｎｄａｒｄＬａｔｉｎｈｙｐｅｒｃｕｂｅ，ａｎｄｔｈｅｎｍａｐｓｔｈｅｍｂａｃｋｔｏｔｈｅｏｒｉｇｉｎａｌｄｅｓｉｇｎ
ｓｐａｃｅ
ａｃｃｏｒｄｉｎｇｔｏｔｈｅｋｎｏｗｎｐｒｏｂａｂｉｌｉｔｙｄｉｓｔｒｉｂｕｔｉｏｎｏｆｉｎｐｕｔｖａｒｉａｂｌｅｓ．ＴｈｅＨｙｐｅｒｃｕｂｅ
ｍｅｔｈｏｄ
ａｎｄｔｈｅｔｒａｄｉｔｉｏｎａｌＬａｔｉｎ
Ｄｅｓｉｇｎ
ｍｅｔｈｏｄｂａｓｅｄ
ｏｎａｌｌ
ａｓｓｕｍｅｄｕｎｉｆｏｒｍｄｉｓｔｒｉｂｕｔｉｏｎ
ａｒｅ
ｃｏｍｐａｒｅｄ
ａｎｄ
ｔｈｅｉｒｅｆｆｅｃｔ
Ｏｉｌ
ｔｈｅ
ｅｓｔａｂｌｉｓｈｅｄＧａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌｅｘａｍｐｌｅｉｎｄｉｃａｔｅｓｔｈｅａｄｖａｎｔａｇｅｏｆｔｈｅｍｅｔｈｏｄ．Ｋｅｙｗｏｒｄｓ：ｕｎｃｅｒｔａｉｎｔｙ
ａｎｄｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓｒｅｓｕｌｔｉｓｄｉｓｃｕｓｓｅｄ．Ａｎ
ａｎａｌｙｓｉｓ；ＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ；Ｇａｕｓｓｉａｎｐｒｏｃｅｓｓ；ｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌ
ｏ引言
工程中大多数复杂数值模拟虽然模型本身具有确定性，但某些输入可能存在不确定性（随机性或模糊性），由此产生的输出不确定性量化问题称为
磊篙票篱釜誓篇麓嬲≥罂羞喜
ｈｔｔｐ：／／ｗｗｗ．ｃｈｉｎａｃａｅ．ｃｎ
单，普适性强，但对于非常耗时的复杂数值模拟来说，直接大量抽样和模拟不现实，必须先建立可快速计算的代理模型，如回归模型、神经网络模型以及径
收稿日期：２０１０．０８．１７修回日期：２０１１・０１．１３
基金项目：国家自然科学基金委员会一中国工程物理研究院联合基金（１０８７６１００）
作者简介：刘信恩（１９７６一），男，四川南充人。副研究员，博士，研究方向为结构动力学，（Ｅ－ｍａｉｌ）ｌｉｕ—ｘｉｎｅｎ＠ｙａｈｏｏ．ｃｏｍ．ｃｎ
万方数据
１０２
计算机辅助工程
２０１１生
向基模型等．
ｍ（工）＝ｈ（工）１ｐ
（２）近年来，一种基于贝叶斯原理的新型代理模
ｃ（ｘ，工’）＝ｅｘｐ｛一（工一ｚ’）１．ｆ２（工一工’）｝
（３）
型——高斯过程响应面模型正日益引起重视．它有
式中：ｈ（ｘ）为需指定的回归函数基矢量；卢为回归三大优点：（１）灵活性好，属于非参数／半参数方法，系数矢量，均有ｇ个分量；盯２为方差；ｎ＝ｄｉａｇ（∞。，易逼近高度非线性函数；（２）精度高，属于插值型方∞：，…，∞。）为对角矩阵，其中的正参数∞＝［∞，，法，能通过所有已知数据点；（３）能量化自身不确定∞：，…，∞。］１’被称为粗糙度参数．通常，∞被处理为性，以后验分布形式描述．该模型虽然已在灵敏度分已知参数，由人为指定或者根据数据估计．
析…、不确定性分析［２］、参数校准或模型确认【３４１方根据定义，对于概率分布已知的高斯过程，任意面得到应用，但其在先验期望函数形式选取、粗糙度ｎ个输入点Ｄ＝｛工。，工：，…，工。｝（称为设计点）处的
参数取值等技术细节方面由于处理方式不同而可能输出ｄ＝［７１（工。），叩（Ｘ：），…，７１（ｚ。）］Ｔ（称为数据）
导致产生不同结果，这仍时常困扰着使用者．刘信恩服从多元正态分布，即
等【５３对上述不同处理方法产生的影响进行理论分ｄｌ卢，矿２一Ⅳ（邵，盯２Ａ）
（４）
析，并结合算例从计算简便性、数值稳定性和结果保式中：Ｈ＝ｈ（Ｄ）’＝［ｈ（工，），ｈ（工２），…，ｈ（ｘ。）］’，守性等角度给出推荐处理方法，初步扫清该方法向Ａ＝ｃ（Ｄ，Ｄ）＝［ｃ（ｘ；，ｘ，）］．当ｄ已知时，函数，，（・）实际工程应用推广的障碍．在此基础上，刘信恩的条件分布仍为一个高斯过程（卡尔曼滤波器），即等【６０进一步研究基于高斯过程响应面模型的贝叶叩（・）ＩＪＢ，盯２，ｄ—ｃＰ（ｍ＋（・），盯２ｃ＋（・，・））（５）
斯参数校准或模型确认问题，提出一种更彻底的简式中：
化方法，不仅让工程分析人员更易理解和实现，且效ｌｉｆｔ‘（ｘ）＝ｈ（ｘ）啊＋ｔ（ｘ）ＴＡ。（ｄ一邱）（６）
果几乎不变，同时，数值更稳定、效率更高，进一步推ｃ‘（，，ｘ’）＝ｃ（ｘ，工７）一ｔ（ｘ）ＴＡ．１ｔ（工’）
（７）
动该方法走向实际工程应用．
式中：ｔ（工）＝ｃ（ｚ，Ｄ）１‘＝［ｃ（工，Ｘ１），ｃ（ｚ，工２），…，
本文简要介绍基于高斯过程响应面模型进行不ｃ（ｘ，Ｘ。）］Ｔ．
确定性分析的基本理论，针对设计点选择问题，提出超参数ＪＢ和盯２的无信息先验分布常取为一种基于输入变量已知概率分布的拉丁超立方设计ｐ（ＪＢ，盯２）ｏｃ盯～，其后验分布根据贝叶斯公式方法，并结合算例对该方法与基于假设均匀分布的ｐ（Ｐ，矿２Ｉｄ）ｏｃＰ（ｄＩＪＢ，矿２）Ｐ（ｐ，盯２）计算，为一个
传统方法展开比较、研究，探讨其对所建立的高斯过正态一逆ｒ分布（自由度为ｎ—ｇ），即
程响应面模型和不确定性分析结果的影响．ｐ
Ｉ，，ｄ一Ⅳ（后，盯２（ＨＴＡ‘１日）。）
（８）盯２
Ｉ
１
高斯过程响应面法基本理论
ｄ一（，Ｉ—ｑ一２）＾２，…－２。
（９）
式中：
确定性的复杂数值模拟模型可视为隐式函数，
卢＝（ＨＴＡ＿日）一百ｒＡ～ｄ
即具体表达式未知的确定性函数ｙ＝田（ｚ）．通常，输入工为Ｐ维矢量，输出Ｙ为标量（或矢量的任意分舀ｚ：虫￡型塑堕窆掣坳（１１）
（１０）
ｎ—ｑ一二
量）．假设叼（・）的先验分布是个平稳高斯过程，即
将式（５），（８）和（９）相乘并对超参数卢和矿２积分，，７（・）ＩＪ６ｌ，矿２一ＧＰ（ｍ（・），盯２ｃ（・，・））（１）可求出函数田（・）的后验分布，这是一个自由度为式中：先验期望函数ｍ（・）和相关函数ｃ（・，・）常ｎ—ｇ的ｔ过程（类似于高斯过程，但需用多元ｔ分布取为
叼（・）Ｉ
ｄ—ＳＰ（ｎ一口，ｍ”（・），生型彦２ｃ”（・，・））
描述），即
（１２）
ｎ—ｇ
式中：
ｍ“（工）＝ｈ（ｘ）１ｐ＋ｆ（ｚ）ＴＡ－１（ｄ一正够）
（１３）ｃ”（Ｘ，Ｘ’）＝ｃ’（Ｘ，Ｘ’）＋（ｈ（ｘ）一百ｒＡ。ｔ（ｘ））’・（日ｒＡ。１日）。１（ｈ（ｘ’）一日’Ａ。１ｔ（ｘ’））
（１４）
通常，后验期望ｍ”（・）用作未知函数叼（・）的预文献［５］对高斯过程响应面法若干技术细节的测（点估计），而后验（协）方差６－２ｃ”（・，・）用于不同处理方法进行研究，并给出先验期望函数选用误差估计（后验方差就是给定输入下的均方误差），低阶多项式、粗糙度参数使用边缘后验众数法估计、也常用某高置信水平（如９５％）下置信区间形式描模型有效性必须经过验算点（新数据）检验等处理述预测值的可能范围．
方法．粗糙度参数∞的边缘后验密度函数为
万方数据
ｈｔｔｐ．／／ｗｗｗ．ｅｈｉｎａｅａｅ．ｅｌｌ
第ｌ期
ｐ（∞Ｉｄ）优Ｐ（∞）Ｉ
刘信恩，等：用于不确定性分析的高斯过程响应面模型设计点选择方法
ＡＩ—Ｔ
１０３
１日‘Ａ－１日Ｉ—Ｔ（存２）一７
（１５）
实际上不可能，因此常用“模拟设计点”技术［２１近似：预先选择不同于原始设计点Ｄ＝｛ｚ，，工：，…，工。｝的若干“模拟设计点”Ｄ＝｛ｉ，，ｉ：，…，ｉ。｝，先抽取出“模拟设计点”上的函数实现ｄｆ＝（叼ｆ（ｉ。），叼ｆ（ｊ：），…，叩ｆ（ｊ。））’，然后将其加入原始设计点集，只要“模拟设计点”足够多且分布合理，其后验方差总能小到足以忽略，此时就可用其后验期望近似任何其他输入Ｊ处的函数实现叼ｉ（工）．该过程不必重建高斯过程响应面模型，可直接通过简单迭代公式计算，即
其最大值所对应的∞值就是边缘后验众数估计值．
∞的无信息先验分布【纠可简单取为Ｐ（∞）ｏｃ
１／（∞ｌ×∞２×…×∞。），经验表明，其先验分布的不同取法一般对结果影响不大．２
基于高斯过程响应面模型的不确定性分析
对于已知确定性函数Ｙ＝叼（工），如果将具有不
确定性的输人工记为随机矢量ｘ（设概率分布已知），则输出Ｙ就成为随机变量Ｙ＝叼（Ｘ）．蒙特卡罗方法（简单随机抽样）首先从ｘ的概率分布中随机
，，ｌｆ”（工）＝ｍ”（工）＋，（工）１Ａ。（ｄｊ一丘）
（１６）
式中：矢量ｔ（工）的第ｉ个元素为ｃ”（ｘ，ｊｉ）；正的第ｉ个元素为ｍ”（ｉ；）；矩阵Ａ的第ｉ行第．『列为ｃ”（ｊｉ，ｊＡ
３
抽取输入样本｛Ｘ１＂，工；，…，石；｝（Ⅳ足够大），然后计
算相应的输出样本｛），？＝叼（ＸＩ’），），ｆ＝叼（卫ｆ），…，，，；＝田（ｘ毒）｝，最后对上述输出样本进行统计分析．
高斯过程响应面法将未知确定性函数田（・）视
为一个随机函数（它的每个实现都是对叼（・）的一个近似），因此输出ｙ的任何总体数字特征（如期望、方差等）都为随机变量，而概率分布和密度函数等为随机函数．可在２个层次上嵌套使用蒙特卡罗方法，即首先从田（・）的后验分布中随机抽取Ｍ个实现｛叼。（・），叼：（・），…，叼村（・）｝（Ｍ足够大），然后对每个实现叼ｊ（・）（歹＝１，２，…，Ｍ）分别使用蒙特卡罗方法进行不确定性分析，最后对不同实现下的分析结果进行统计分析．不仅能获得输出ｙ的期望、方差或概率分布／密度函数等传统的不确定性分析结果，而且还能获得对这些分析结果自身不确定性的估计．
获得随机函数叼（・）的一个精确实现叼，（・）
１．０
高斯过程响应面模型设计点选择
建立高斯过程响应面模型的前提为获得设计点
Ｄ＝｛ｘ。，工：，…，工。｝处的数据ｄ＝［叼（工１），叼（工：），…，叼（ｚ。）］Ｔ．一般使用“空间填充”的试验
设计方法选择设计点（如拉丁超立方设计方法［７１）．该方法将Ｐ维输人中的每一维等概率地分割为ｎ个子区间，每个子区间内依概率随机抽取一个样本作为分点（为简单起见，也可直接使用中点或端点，本文采用该处理方式），然后将不同维上的分点不重复地随机组合形成设计点．该方法设计点分布比较均匀，且投影到任何一维都不重复，代表性强，加密后仍为拉丁超立方设计，很受欢迎．图ｌ为一个二维拉丁超立方设计的例子．
Ｏ．８
０．６犁
０．４
点
密密
Ｏ．２
０
ｏ．２ｏ．４
Ｘｌ
ｏ．６ｏ．８１．０
（ａ）随机取点
（ｂ）取中点（ｃ）取端点
图１一个二维拉丁超立方设计的例子。，ｌ＝５
Ｆｉｇ．１
Ａｎ
ｅｘａｍｐｌｅｏｆ２ＤＬａｔｉｎｈｙｐｅｒｃｕｂｅｄｅｓｉｇｎ。，ｌ＝５
但在实践中，通常是先在标准超立方体（ｊｊ∈［０，１］Ｊ＝１，２，…，Ｐ）中生成拉丁超立方设计点，然
后映射回原始设计空间．具体是：标准超立方体中的
每个设计点ｉｉ＝（ｉｎ＇牙垃，…，ｉｌ；Ｐ）７（ｉ＝１，２，…，ｎ）
万方数据
ｈｔｔｐ：／／ｗｗｗ．ｃｈｉｎａｃａｅ．ｃｎ
１０４
计算机辅助工程
４．１高斯过程响应面模型
２０１１年
的第＿『个分量互ｉ通过其第Ｊ维边缘分布函数‘（・）的逆函数映射回原始设计空间，即并。＝Ｆ／１（ｉ。）．这里所用的概率分布可任意假定，一般使用均匀分布．考虑到开展不确定性分析时输入变量的概率分布已知，本文提出一种新方法，直接基于输入变量的已知概率分布从标准超立方体映射回原始设计空间，即Ｅ（・）使用已知的输入变量边缘分布函数．
为简单起见，本文将基于均匀分布的方法称为均匀取点方法，将基于已知概率分布的方法称为依概率取点方法．均匀取点方法的难点在于如何确定感兴趣的区间，因为不确定性分析常假定输出变量服从正态分布、对数正态分布等无限或半无限区间上的概率分布，基于均匀分布进行映射必须将其截断为有限区间，这是一个两难问题：如果使绝大多数样本点落在不确定性较小的内插区，就必须选择较
选取８个设计点（ｒｔ＝８），在均匀取点时取子区间端点，依概率取点时取子区间中点．高斯过程响应面模型先验期望函数选用ｌ阶多项式，出于演示目的，直接指定粗糙度参数ｔｏ＝０．５（如果使用推荐的边缘后验众数法，估计出的∞值要小许多，函数预测的置信区间太窄，不便于图示，因此本文作了保守处理）．
图２为在２种设计点选择方法下建立的高斯过程响应面模型．在均匀取点时，感兴趣区间取为［一４，４］，即４０＂范围（概率大于９９．９９％），在进行几万次蒙特卡罗抽样时几乎能覆盖所有样本点．可知，此时真实曲线被９５％的置信区间覆盖，但在整个感兴趣区间内，函数叼（・）的不确定性仍比较明显．如果选择３０－或２０－范围（概率分别为９９．７３％和９５．４４％），那么随着设计点间距减小，内插区的不确定性更小而外插区的不确定性明显增大，同时也会有更多样本点落到外插区．与均匀取点相比，依概率取点时９５％的置信区间同样覆盖了真实曲线，且内插区精度极高，置信区间宽度几乎为０，但外插区不确定性明显大得多，同时由于内插区范围比较窄（ｎ＝８时还不足２０＂范围），会有较多的样本点落到外插区．当然，随着设计点增多，其内插区范围会自动拓展从而缓解该问题，但同时高概率区的设计点也可能过密，不仅造成浪费，而且会带来数值上的困难（相关矩阵Ａ接近奇异，难以求逆），因此必要时须舍弃某些设计点．
１０８６４２＾０—２—４—６—８
大的范围，但设计点数量较少又可能使设计点间距
过大，内插区不确定性增大，反而使最终的不确定性分析结果精度不高，反之亦然．与均匀取点方法相比，依概率取点方法不需要人为指定感兴趣的区间，且自动在概率密度较大的区域安排较密的设计点以体现其相对重要性，实现设计点位置的某种“优化”，具有一定优势．
４算例和讨论
假设未知函数Ｙ＝叼（髫）＝１＋省＋ＣＯＳ算（代表实际应用中的复杂数值模拟），输入ｘ—Ｎ（０，１），建立高斯过程响应面模型并用蒙特卡罗方法估计输出ｙ＝ｎ（Ｘ）的概率分布．
１０８６４２＾０
—２—４—６—８
（ａ）均匀取点（ｂ）依概率取点
图２高斯过程响应面模型
Ｆｉｇ．２
Ｇａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌｓ
４．２不确定性分析的结果
图３为２种设计点选择方法下输出ｙ的概率分布函数（通过Ｋｅｒｎｅｌ－Ｓｍｏｏｔｈｉｎｇ方法拟合）．输入Ｘ的样本容量Ⅳ和随机函数叼（・）的样本容量肘均取１００００，经重复计算表明，在该样本容量下抽样所
带来的不确定性已很小，可忽略．由图３可知，在２种设计点选择方法下，输出ｌ，的概率分布函数在９５％的置信区间均覆盖其真实值，而且其期望值也与真实值非常接近（该期望值一般作为对概率分布函数的点估计）。相比而言，在均匀取点时概率分布
ｈｔｔｐ：／／ｗｗｗ．ＣｈＪｌ忸Ｃｇｌｅ．ｃ．ｎ万方数据
第１期刘信恩，等：用于不确定性分析的高斯过程响应面模型设计点选择方法
１０５
函数的不确定性更大一些，依概率取点时仅在两端低概率区表现出轻微的不确定性．尝试在均匀取点时将感兴趣的区间缩小为３０＂或２０－范围，结果发现所获得的概率分布函数不确定性减小．其中，当取２０＂范围时，与依概率取点时结果差不多（感兴趣区间取多大合适还与设计点数量、输入维数等有关，不
能轻易下结论）．此外，直接根据高斯过程响应面模型后验期望通过蒙特卡罗方法获得的概率分布函数也与真实值吻合很好，因为高斯过程响应面模型后验期望函数ｍ”（・）常常为未知函数田（・）的一个好的点估计（见图２）．
（ａ）均匀取点（ｂ）依概率取点
图３输出ｙ的概率分布函数
Ｆｉｇ．３
Ｐｍｂａ砌坶ｄｉｓｔｒｉｂｕｔｉｏｎ
ｆｕｎｃｔｉｏｎｓｏｆｏｕｔｐｕｔｖａｒｉａｂｌｅＹ
值得说明的是，Ｋｅｒｎｅｌ・Ｓｍｏｏｔｈｉｎｇ方法拟合出的是光滑的概率分布／密度函数，而本例中函数田（・）在Ｙ＝２．５７处有一个拐点，在该处输出ｙ的概率密度函数奇异，即其概率分布函数斜率为无穷大，拟合结果在该处有明显误差．５
期望、方差或概率分布／密度函数等分析结果自身不确定性的独特优势．
（２）针对高斯过程响应面模型设计点选择这个重要问题，提出一种基于输入变量已知概率分布的拉丁超立方设计方法，与基于假设均匀分布的传统方法相比，新方法不仅可以避开人为指定感兴趣区间的难题，而且能自动实现设计点位置优化．在相同数量设计点情况下，不确定性分析结果精度更高，更具优势．
结论
（１）基于高斯过程响应面模型和蒙特卡罗方法
展开复杂数值模拟不确定性分析，具有可量化输出
参考文献：
［１］ＯＡＫＬＥＹＪＥ，０’ＨＡＧＡＮＡ．Ｐｒｏｂａｂｉｌｉｓｔｉｃ
Ｍｅｔｈｏｄｏｌｏｇｙ。２００４，６６（３）：７５１－７６９．［２］ＯＡＫＬＥＹＪ
７６９．７８４．
Ｅ，Ｏ’ＨＡＧＡＮＡ．Ｂａｙｅｓｉａｎ
ｉｎｆｅｒｅｎｃｅｆｏｒｔｈｅｕｎｃｅｒｔａｉｎｔｙ
ｓｅｎｓｉｔｉｖｉｔｙ
ａｎａｌｙｓｉｓ
ｏｆ
ｃｏｍｐｌｅｘｍｏｄｅｌ：ａ
Ｂａｙｅｓｉａｎ
ａｐｐｒｏａｃｈ［Ｊ］．Ｊ
ＲＳｔａｔ
Ｓｏｅ：Ｓｅｒ
Ｂ：Ｓｔａｒ
ｄｉｓｔｒｉｂｕｔｉｏｎｏｆｃｏｍｐｕｔｅｒｍｏｄｅｌｏｕｔｐｕｔｓ［Ｊ］．Ｂｉｏｍｅｔｒｉｋａ，２００２，８９（４）：
［３］ＫＥＮＮＥＤＹ［４］ＢＡＹＡＲＲＩ
Ｒｅｐ，２００５．
ＭＭ
Ｃ，Ｏ’ＨＡＧＡＮＡ．ＢａｙｅｓｉａｎＪ，ＢＥＲＧＥＲＪ
Ｏ，ＰＡＵＬＯ
ｃａｌｉｂｒａｔｉｏｎｏｆｃｏｍｐｕｔｅｒ
ｆ／ｉｔｍｅｗｏｌ＇ｋ
ｍｏｄｅｌ［Ｊ］．Ｊ
ｖａｌｉｄａｔｉｏｎ
ｏｆ
ＲＳｔａｔＳｏｃ：ＳｅｘＢ：Ｓｔａｔ
Ｍｅｔｈｏｄｏｌｏｇｙ，２００１，６３（３）：４２５－４６４．
ＩｎｓｔＳｔａｒ
Ｒ，ｅｔａ／．Ａ
ｆｏｒ
ｃｏｍｐｕｔｅｒ
ｍｏｄｅｌｓ［Ｒ］．ＮｏｒｔｈＣａｒｏｌｉｎａ：Ｎａｔｌ
ＳｃｉＴｅｃｈ
［５］刘信恩，肖世富，莫军．高斯过程响应面法研究［Ｊ］．应用力学学报，２０１０，２７（１）：１９０—１９５．
ＬＩＵＸｉｎ’ｅｎ，ＸＩＡＯＳｈｉｆｕ，ＭＯＪｕｎ．Ａ
ｎｅｗ
ｒｅｓｐｏｎｓｅ
ｓｕｒｆａｃｅｍｅｔｈｏｄｂａｓｅｄｏｎＧａｕｓｓｉａｎ
ｐｒｏｃｅｓｓ［Ｊ］．ＣｈｉｎＪＡｐｐＩＭｅｃｈ，２０１０，２７（１）：１９０—１９５．
［６］
刘信恩，肖世富，莫军．复杂数值模拟的贝叶斯模型确认框架及其简化［ｃ］∥中国力学学会学术大会’２００９，ｃｃＴＡ№００９瑚３７８９。
郑州，２００９．
Ｍ
Ｄ，ＭＩＴＣＨＥＬＬＴ
［７］ＭＯＲＲＩＳＪ．Ｅｘｐｌｏｒａｔｏｒｙｄｅｓｉｇｎｓ
ｆｏｒｃｏｍｐｕｔａｔｉｏｎａｌｅｘｐｅｒｉｍｅｎｔｓ［Ｊ］．Ｊ
Ｓｔａｒ
Ｐｌａｎｎｉｎｇ＆Ｉｎｆｅｒｅｎｃｅ。１９９５，４３（３）：３８１－４０２．
（编辑陈锋杰）
万方数据
ｈｔｔｐ：／／ｗｗｗ．ｅＭｎａｃａｅ．ｃｎ
第２０卷第１期
２０１１年３月
计算机辅助工程
ＣｏｍｐｕｔｅｒＡｉｄｅｄＥｎｇｉｎｅｅｒｉｎｇ
Ｖ０１．２０Ｎｏ．１Ｍａｒ．２０１１
文章编号：１００６—０８７１（２０１１）０１．０１０１．０５
用于不确定性分析的高斯过程响应面模型
设计点选择方法
刘信恩，
肖世富，莫军
（中国工程物理研究院总体工程研究所，四川绵阳６２１９００）
摘要：为推动高斯过程响应面模型在复杂耗时数值模拟不确定性分析中的应用，提出一种可自动实现位置优化的高效设计点选择方法，即先在标准超立方体上生成拉丁超立方设计点，然后利用输入变量的已知概率分布将其映射回原始设计空间．将该方法与基于假设均匀分布的传统拉丁超立方设计方法进行比较，探讨它们对所建立的高斯过程响应面模型和不确定性分析结果的影响．算例表明新方法具有一定优势．
关键词：不确定性分析；蒙特卡罗方法；高斯过程；响应面模型中图分类号：０２１２．８；０２１１．９
文献标志码：Ａ
ＤｅｓｉｇｎｐｏｉｎｔｃｈｏｉｃｅｏｆＧａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌａｐｐｌｉｅｄｔｏｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓ
ＬＩＵＸｉｎ’ｅｎ，ＸＩＡＯＳｈｉｆｕ，ＭＯＪｕｎ
（Ｉｎｓｔｉｔｕｔｅ
ｏｆＳｙｓｔｅｍｓＥｎｇｉｎｅｅｒｉｎｇ，ＣｈｉｎａＡｃａｄｅｍｙｏｆ
Ｅｎｇｉｎｅｅｒｉｎｇ
Ｐｈｙｓｉｃｓ，Ｍｉａｎｙａｎｇ
６２１９００，Ｓｉｃｈｕａｎ，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ａｎｅｆｆｉｃｉｅｎｔｍｅｔｈｏｄｏｆｄｅｓｉｇｎｐｏｉｎｔｃｈｏｉｃｅｉｓｐｒｅｓｅｎｔｅｄ，ｗｈｉｃｈｌｏｃａｔｉｏｎ
ｏｆｔｈｅ
ｐｏｉｎｔｓ，ｔｏ
ａｄｖａｎｃｅｃｏｍｐｌｅｘ
ｔｈｅ
ｃａｎ
ａｕｔｏｍａｔｉｃａｌｌｙｏｐｔｉｍｉｚｅｔｈｅ
ｓｕｒｆａｃｅ
ｍｏｄｅｌ
ｔｏ
ａｐｐｌｉｃａｔｉｏｎｏｆＧａｕｓｓｉａｎｐｒｏｃｅｓｓ
ｒｅｓｐｏｎｓｅ
ｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓｏｆｔｈｅｄｅｓｉｇｎｐｏｉｎｔｓｉｎ
ａ
ｔｉｍｅ—ｃｏｎｓｕｍｉｎｇｎｕｍｅｒｉｃａｌ
ｓｉｍｕｌａｔｉｏｎ．Ｔｈｅｍｅｔｈｏｄ
ｐｒｏｄｕｃｅｓ
ｔｈｅ
ｓｔａｎｄａｒｄＬａｔｉｎｈｙｐｅｒｃｕｂｅ，ａｎｄｔｈｅｎｍａｐｓｔｈｅｍｂａｃｋｔｏｔｈｅｏｒｉｇｉｎａｌｄｅｓｉｇｎ
ｓｐａｃｅ
ａｃｃｏｒｄｉｎｇｔｏｔｈｅｋｎｏｗｎｐｒｏｂａｂｉｌｉｔｙｄｉｓｔｒｉｂｕｔｉｏｎｏｆｉｎｐｕｔｖａｒｉａｂｌｅｓ．ＴｈｅＨｙｐｅｒｃｕｂｅ
ｍｅｔｈｏｄ
ａｎｄｔｈｅｔｒａｄｉｔｉｏｎａｌＬａｔｉｎ
Ｄｅｓｉｇｎ
ｍｅｔｈｏｄｂａｓｅｄ
ｏｎａｌｌ
ａｓｓｕｍｅｄｕｎｉｆｏｒｍｄｉｓｔｒｉｂｕｔｉｏｎ
ａｒｅ
ｃｏｍｐａｒｅｄ
ａｎｄ
ｔｈｅｉｒｅｆｆｅｃｔ
Ｏｉｌ
ｔｈｅ
ｅｓｔａｂｌｉｓｈｅｄＧａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌｅｘａｍｐｌｅｉｎｄｉｃａｔｅｓｔｈｅａｄｖａｎｔａｇｅｏｆｔｈｅｍｅｔｈｏｄ．Ｋｅｙｗｏｒｄｓ：ｕｎｃｅｒｔａｉｎｔｙ
ａｎｄｕｎｃｅｒｔａｉｎｔｙａｎａｌｙｓｉｓｒｅｓｕｌｔｉｓｄｉｓｃｕｓｓｅｄ．Ａｎ
ａｎａｌｙｓｉｓ；ＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ；Ｇａｕｓｓｉａｎｐｒｏｃｅｓｓ；ｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌ
ｏ引言
工程中大多数复杂数值模拟虽然模型本身具有确定性，但某些输入可能存在不确定性（随机性或模糊性），由此产生的输出不确定性量化问题称为
磊篙票篱釜誓篇麓嬲≥罂羞喜
ｈｔｔｐ：／／ｗｗｗ．ｃｈｉｎａｃａｅ．ｃｎ
单，普适性强，但对于非常耗时的复杂数值模拟来说，直接大量抽样和模拟不现实，必须先建立可快速计算的代理模型，如回归模型、神经网络模型以及径
收稿日期：２０１０．０８．１７修回日期：２０１１・０１．１３
基金项目：国家自然科学基金委员会一中国工程物理研究院联合基金（１０８７６１００）
作者简介：刘信恩（１９７６一），男，四川南充人。副研究员，博士，研究方向为结构动力学，（Ｅ－ｍａｉｌ）ｌｉｕ—ｘｉｎｅｎ＠ｙａｈｏｏ．ｃｏｍ．ｃｎ
万方数据
１０２
计算机辅助工程
２０１１生
向基模型等．
ｍ（工）＝ｈ（工）１ｐ
（２）近年来，一种基于贝叶斯原理的新型代理模
ｃ（ｘ，工’）＝ｅｘｐ｛一（工一ｚ’）１．ｆ２（工一工’）｝
（３）
型——高斯过程响应面模型正日益引起重视．它有
式中：ｈ（ｘ）为需指定的回归函数基矢量；卢为回归三大优点：（１）灵活性好，属于非参数／半参数方法，系数矢量，均有ｇ个分量；盯２为方差；ｎ＝ｄｉａｇ（∞。，易逼近高度非线性函数；（２）精度高，属于插值型方∞：，…，∞。）为对角矩阵，其中的正参数∞＝［∞，，法，能通过所有已知数据点；（３）能量化自身不确定∞：，…，∞。］１’被称为粗糙度参数．通常，∞被处理为性，以后验分布形式描述．该模型虽然已在灵敏度分已知参数，由人为指定或者根据数据估计．
析…、不确定性分析［２］、参数校准或模型确认【３４１方根据定义，对于概率分布已知的高斯过程，任意面得到应用，但其在先验期望函数形式选取、粗糙度ｎ个输入点Ｄ＝｛工。，工：，…，工。｝（称为设计点）处的
参数取值等技术细节方面由于处理方式不同而可能输出ｄ＝［７１（工。），叩（Ｘ：），…，７１（ｚ。）］Ｔ（称为数据）
导致产生不同结果，这仍时常困扰着使用者．刘信恩服从多元正态分布，即
等【５３对上述不同处理方法产生的影响进行理论分ｄｌ卢，矿２一Ⅳ（邵，盯２Ａ）
（４）
析，并结合算例从计算简便性、数值稳定性和结果保式中：Ｈ＝ｈ（Ｄ）’＝［ｈ（工，），ｈ（工２），…，ｈ（ｘ。）］’，守性等角度给出推荐处理方法，初步扫清该方法向Ａ＝ｃ（Ｄ，Ｄ）＝［ｃ（ｘ；，ｘ，）］．当ｄ已知时，函数，，（・）实际工程应用推广的障碍．在此基础上，刘信恩的条件分布仍为一个高斯过程（卡尔曼滤波器），即等【６０进一步研究基于高斯过程响应面模型的贝叶叩（・）ＩＪＢ，盯２，ｄ—ｃＰ（ｍ＋（・），盯２ｃ＋（・，・））（５）
斯参数校准或模型确认问题，提出一种更彻底的简式中：
化方法，不仅让工程分析人员更易理解和实现，且效ｌｉｆｔ‘（ｘ）＝ｈ（ｘ）啊＋ｔ（ｘ）ＴＡ。（ｄ一邱）（６）
果几乎不变，同时，数值更稳定、效率更高，进一步推ｃ‘（，，ｘ’）＝ｃ（ｘ，工７）一ｔ（ｘ）ＴＡ．１ｔ（工’）
（７）
动该方法走向实际工程应用．
式中：ｔ（工）＝ｃ（ｚ，Ｄ）１‘＝［ｃ（工，Ｘ１），ｃ（ｚ，工２），…，
本文简要介绍基于高斯过程响应面模型进行不ｃ（ｘ，Ｘ。）］Ｔ．
确定性分析的基本理论，针对设计点选择问题，提出超参数ＪＢ和盯２的无信息先验分布常取为一种基于输入变量已知概率分布的拉丁超立方设计ｐ（ＪＢ，盯２）ｏｃ盯～，其后验分布根据贝叶斯公式方法，并结合算例对该方法与基于假设均匀分布的ｐ（Ｐ，矿２Ｉｄ）ｏｃＰ（ｄＩＪＢ，矿２）Ｐ（ｐ，盯２）计算，为一个
传统方法展开比较、研究，探讨其对所建立的高斯过正态一逆ｒ分布（自由度为ｎ—ｇ），即
程响应面模型和不确定性分析结果的影响．ｐ
Ｉ，，ｄ一Ⅳ（后，盯２（ＨＴＡ‘１日）。）
（８）盯２
Ｉ
１
高斯过程响应面法基本理论
ｄ一（，Ｉ—ｑ一２）＾２，…－２。
（９）
式中：
确定性的复杂数值模拟模型可视为隐式函数，
卢＝（ＨＴＡ＿日）一百ｒＡ～ｄ
即具体表达式未知的确定性函数ｙ＝田（ｚ）．通常，输入工为Ｐ维矢量，输出Ｙ为标量（或矢量的任意分舀ｚ：虫￡型塑堕窆掣坳（１１）
（１０）
ｎ—ｑ一二
量）．假设叼（・）的先验分布是个平稳高斯过程，即
将式（５），（８）和（９）相乘并对超参数卢和矿２积分，，７（・）ＩＪ６ｌ，矿２一ＧＰ（ｍ（・），盯２ｃ（・，・））（１）可求出函数田（・）的后验分布，这是一个自由度为式中：先验期望函数ｍ（・）和相关函数ｃ（・，・）常ｎ—ｇ的ｔ过程（类似于高斯过程，但需用多元ｔ分布取为
叼（・）Ｉ
ｄ—ＳＰ（ｎ一口，ｍ”（・），生型彦２ｃ”（・，・））
描述），即
（１２）
ｎ—ｇ
式中：
ｍ“（工）＝ｈ（ｘ）１ｐ＋ｆ（ｚ）ＴＡ－１（ｄ一正够）
（１３）ｃ”（Ｘ，Ｘ’）＝ｃ’（Ｘ，Ｘ’）＋（ｈ（ｘ）一百ｒＡ。ｔ（ｘ））’・（日ｒＡ。１日）。１（ｈ（ｘ’）一日’Ａ。１ｔ（ｘ’））
（１４）
通常，后验期望ｍ”（・）用作未知函数叼（・）的预文献［５］对高斯过程响应面法若干技术细节的测（点估计），而后验（协）方差６－２ｃ”（・，・）用于不同处理方法进行研究，并给出先验期望函数选用误差估计（后验方差就是给定输入下的均方误差），低阶多项式、粗糙度参数使用边缘后验众数法估计、也常用某高置信水平（如９５％）下置信区间形式描模型有效性必须经过验算点（新数据）检验等处理述预测值的可能范围．
方法．粗糙度参数∞的边缘后验密度函数为
万方数据
ｈｔｔｐ．／／ｗｗｗ．ｅｈｉｎａｅａｅ．ｅｌｌ
第ｌ期
ｐ（∞Ｉｄ）优Ｐ（∞）Ｉ
刘信恩，等：用于不确定性分析的高斯过程响应面模型设计点选择方法
ＡＩ—Ｔ
１０３
１日‘Ａ－１日Ｉ—Ｔ（存２）一７
（１５）
实际上不可能，因此常用“模拟设计点”技术［２１近似：预先选择不同于原始设计点Ｄ＝｛ｚ，，工：，…，工。｝的若干“模拟设计点”Ｄ＝｛ｉ，，ｉ：，…，ｉ。｝，先抽取出“模拟设计点”上的函数实现ｄｆ＝（叼ｆ（ｉ。），叼ｆ（ｊ：），…，叩ｆ（ｊ。））’，然后将其加入原始设计点集，只要“模拟设计点”足够多且分布合理，其后验方差总能小到足以忽略，此时就可用其后验期望近似任何其他输入Ｊ处的函数实现叼ｉ（工）．该过程不必重建高斯过程响应面模型，可直接通过简单迭代公式计算，即
其最大值所对应的∞值就是边缘后验众数估计值．
∞的无信息先验分布【纠可简单取为Ｐ（∞）ｏｃ
１／（∞ｌ×∞２×…×∞。），经验表明，其先验分布的不同取法一般对结果影响不大．２
基于高斯过程响应面模型的不确定性分析
对于已知确定性函数Ｙ＝叼（工），如果将具有不
确定性的输人工记为随机矢量ｘ（设概率分布已知），则输出Ｙ就成为随机变量Ｙ＝叼（Ｘ）．蒙特卡罗方法（简单随机抽样）首先从ｘ的概率分布中随机
，，ｌｆ”（工）＝ｍ”（工）＋，（工）１Ａ。（ｄｊ一丘）
（１６）
式中：矢量ｔ（工）的第ｉ个元素为ｃ”（ｘ，ｊｉ）；正的第ｉ个元素为ｍ”（ｉ；）；矩阵Ａ的第ｉ行第．『列为ｃ”（ｊｉ，ｊＡ
３
抽取输入样本｛Ｘ１＂，工；，…，石；｝（Ⅳ足够大），然后计
算相应的输出样本｛），？＝叼（ＸＩ’），），ｆ＝叼（卫ｆ），…，，，；＝田（ｘ毒）｝，最后对上述输出样本进行统计分析．
高斯过程响应面法将未知确定性函数田（・）视
为一个随机函数（它的每个实现都是对叼（・）的一个近似），因此输出ｙ的任何总体数字特征（如期望、方差等）都为随机变量，而概率分布和密度函数等为随机函数．可在２个层次上嵌套使用蒙特卡罗方法，即首先从田（・）的后验分布中随机抽取Ｍ个实现｛叼。（・），叼：（・），…，叼村（・）｝（Ｍ足够大），然后对每个实现叼ｊ（・）（歹＝１，２，…，Ｍ）分别使用蒙特卡罗方法进行不确定性分析，最后对不同实现下的分析结果进行统计分析．不仅能获得输出ｙ的期望、方差或概率分布／密度函数等传统的不确定性分析结果，而且还能获得对这些分析结果自身不确定性的估计．
获得随机函数叼（・）的一个精确实现叼，（・）
１．０
高斯过程响应面模型设计点选择
建立高斯过程响应面模型的前提为获得设计点
Ｄ＝｛ｘ。，工：，…，工。｝处的数据ｄ＝［叼（工１），叼（工：），…，叼（ｚ。）］Ｔ．一般使用“空间填充”的试验
设计方法选择设计点（如拉丁超立方设计方法［７１）．该方法将Ｐ维输人中的每一维等概率地分割为ｎ个子区间，每个子区间内依概率随机抽取一个样本作为分点（为简单起见，也可直接使用中点或端点，本文采用该处理方式），然后将不同维上的分点不重复地随机组合形成设计点．该方法设计点分布比较均匀，且投影到任何一维都不重复，代表性强，加密后仍为拉丁超立方设计，很受欢迎．图ｌ为一个二维拉丁超立方设计的例子．
Ｏ．８
０．６犁
０．４
点
密密
Ｏ．２
０
ｏ．２ｏ．４
Ｘｌ
ｏ．６ｏ．８１．０
（ａ）随机取点
（ｂ）取中点（ｃ）取端点
图１一个二维拉丁超立方设计的例子。，ｌ＝５
Ｆｉｇ．１
Ａｎ
ｅｘａｍｐｌｅｏｆ２ＤＬａｔｉｎｈｙｐｅｒｃｕｂｅｄｅｓｉｇｎ。，ｌ＝５
但在实践中，通常是先在标准超立方体（ｊｊ∈［０，１］Ｊ＝１，２，…，Ｐ）中生成拉丁超立方设计点，然
后映射回原始设计空间．具体是：标准超立方体中的
每个设计点ｉｉ＝（ｉｎ＇牙垃，…，ｉｌ；Ｐ）７（ｉ＝１，２，…，ｎ）
万方数据
ｈｔｔｐ：／／ｗｗｗ．ｃｈｉｎａｃａｅ．ｃｎ
１０４
计算机辅助工程
４．１高斯过程响应面模型
２０１１年
的第＿『个分量互ｉ通过其第Ｊ维边缘分布函数‘（・）的逆函数映射回原始设计空间，即并。＝Ｆ／１（ｉ。）．这里所用的概率分布可任意假定，一般使用均匀分布．考虑到开展不确定性分析时输入变量的概率分布已知，本文提出一种新方法，直接基于输入变量的已知概率分布从标准超立方体映射回原始设计空间，即Ｅ（・）使用已知的输入变量边缘分布函数．
为简单起见，本文将基于均匀分布的方法称为均匀取点方法，将基于已知概率分布的方法称为依概率取点方法．均匀取点方法的难点在于如何确定感兴趣的区间，因为不确定性分析常假定输出变量服从正态分布、对数正态分布等无限或半无限区间上的概率分布，基于均匀分布进行映射必须将其截断为有限区间，这是一个两难问题：如果使绝大多数样本点落在不确定性较小的内插区，就必须选择较
选取８个设计点（ｒｔ＝８），在均匀取点时取子区间端点，依概率取点时取子区间中点．高斯过程响应面模型先验期望函数选用ｌ阶多项式，出于演示目的，直接指定粗糙度参数ｔｏ＝０．５（如果使用推荐的边缘后验众数法，估计出的∞值要小许多，函数预测的置信区间太窄，不便于图示，因此本文作了保守处理）．
图２为在２种设计点选择方法下建立的高斯过程响应面模型．在均匀取点时，感兴趣区间取为［一４，４］，即４０＂范围（概率大于９９．９９％），在进行几万次蒙特卡罗抽样时几乎能覆盖所有样本点．可知，此时真实曲线被９５％的置信区间覆盖，但在整个感兴趣区间内，函数叼（・）的不确定性仍比较明显．如果选择３０－或２０－范围（概率分别为９９．７３％和９５．４４％），那么随着设计点间距减小，内插区的不确定性更小而外插区的不确定性明显增大，同时也会有更多样本点落到外插区．与均匀取点相比，依概率取点时９５％的置信区间同样覆盖了真实曲线，且内插区精度极高，置信区间宽度几乎为０，但外插区不确定性明显大得多，同时由于内插区范围比较窄（ｎ＝８时还不足２０＂范围），会有较多的样本点落到外插区．当然，随着设计点增多，其内插区范围会自动拓展从而缓解该问题，但同时高概率区的设计点也可能过密，不仅造成浪费，而且会带来数值上的困难（相关矩阵Ａ接近奇异，难以求逆），因此必要时须舍弃某些设计点．
１０８６４２＾０—２—４—６—８
大的范围，但设计点数量较少又可能使设计点间距
过大，内插区不确定性增大，反而使最终的不确定性分析结果精度不高，反之亦然．与均匀取点方法相比，依概率取点方法不需要人为指定感兴趣的区间，且自动在概率密度较大的区域安排较密的设计点以体现其相对重要性，实现设计点位置的某种“优化”，具有一定优势．
４算例和讨论
假设未知函数Ｙ＝叼（髫）＝１＋省＋ＣＯＳ算（代表实际应用中的复杂数值模拟），输入ｘ—Ｎ（０，１），建立高斯过程响应面模型并用蒙特卡罗方法估计输出ｙ＝ｎ（Ｘ）的概率分布．
１０８６４２＾０
—２—４—６—８
（ａ）均匀取点（ｂ）依概率取点
图２高斯过程响应面模型
Ｆｉｇ．２
Ｇａｕｓｓｉａｎｐｒｏｃｅｓｓｒｅｓｐｏｎｓｅｓｕｒｆａｃｅｍｏｄｅｌｓ
４．２不确定性分析的结果
图３为２种设计点选择方法下输出ｙ的概率分布函数（通过Ｋｅｒｎｅｌ－Ｓｍｏｏｔｈｉｎｇ方法拟合）．输入Ｘ的样本容量Ⅳ和随机函数叼（・）的样本容量肘均取１００００，经重复计算表明，在该样本容量下抽样所
带来的不确定性已很小，可忽略．由图３可知，在２种设计点选择方法下，输出ｌ，的概率分布函数在９５％的置信区间均覆盖其真实值，而且其期望值也与真实值非常接近（该期望值一般作为对概率分布函数的点估计）。相比而言，在均匀取点时概率分布
ｈｔｔｐ：／／ｗｗｗ．ＣｈＪｌ忸Ｃｇｌｅ．ｃ．ｎ万方数据
第１期刘信恩，等：用于不确定性分析的高斯过程响应面模型设计点选择方法
１０５
函数的不确定性更大一些，依概率取点时仅在两端低概率区表现出轻微的不确定性．尝试在均匀取点时将感兴趣的区间缩小为３０＂或２０－范围，结果发现所获得的概率分布函数不确定性减小．其中，当取２０＂范围时，与依概率取点时结果差不多（感兴趣区间取多大合适还与设计点数量、输入维数等有关，不
能轻易下结论）．此外，直接根据高斯过程响应面模型后验期望通过蒙特卡罗方法获得的概率分布函数也与真实值吻合很好，因为高斯过程响应面模型后验期望函数ｍ”（・）常常为未知函数田（・）的一个好的点估计（见图２）．
（ａ）均匀取点（ｂ）依概率取点
图３输出ｙ的概率分布函数
Ｆｉｇ．３
Ｐｍｂａ砌坶ｄｉｓｔｒｉｂｕｔｉｏｎ
ｆｕｎｃｔｉｏｎｓｏｆｏｕｔｐｕｔｖａｒｉａｂｌｅＹ
值得说明的是，Ｋｅｒｎｅｌ・Ｓｍｏｏｔｈｉｎｇ方法拟合出的是光滑的概率分布／密度函数，而本例中函数田（・）在Ｙ＝２．５７处有一个拐点，在该处输出ｙ的概率密度函数奇异，即其概率分布函数斜率为无穷大，拟合结果在该处有明显误差．５
期望、方差或概率分布／密度函数等分析结果自身不确定性的独特优势．
（２）针对高斯过程响应面模型设计点选择这个重要问题，提出一种基于输入变量已知概率分布的拉丁超立方设计方法，与基于假设均匀分布的传统方法相比，新方法不仅可以避开人为指定感兴趣区间的难题，而且能自动实现设计点位置优化．在相同数量设计点情况下，不确定性分析结果精度更高，更具优势．
结论
（１）基于高斯过程响应面模型和蒙特卡罗方法
展开复杂数值模拟不确定性分析，具有可量化输出
参考文献：
［１］ＯＡＫＬＥＹＪＥ，０’ＨＡＧＡＮＡ．Ｐｒｏｂａｂｉｌｉｓｔｉｃ
Ｍｅｔｈｏｄｏｌｏｇｙ。２００４，６６（３）：７５１－７６９．［２］ＯＡＫＬＥＹＪ
７６９．７８４．
Ｅ，Ｏ’ＨＡＧＡＮＡ．Ｂａｙｅｓｉａｎ
ｉｎｆｅｒｅｎｃｅｆｏｒｔｈｅｕｎｃｅｒｔａｉｎｔｙ
ｓｅｎｓｉｔｉｖｉｔｙ
ａｎａｌｙｓｉｓ
ｏｆ
ｃｏｍｐｌｅｘｍｏｄｅｌ：ａ
Ｂａｙｅｓｉａｎ
ａｐｐｒｏａｃｈ［Ｊ］．Ｊ
ＲＳｔａｔ
Ｓｏｅ：Ｓｅｒ
Ｂ：Ｓｔａｒ
ｄｉｓｔｒｉｂｕｔｉｏｎｏｆｃｏｍｐｕｔｅｒｍｏｄｅｌｏｕｔｐｕｔｓ［Ｊ］．Ｂｉｏｍｅｔｒｉｋａ，２００２，８９（４）：
［３］ＫＥＮＮＥＤＹ［４］ＢＡＹＡＲＲＩ
Ｒｅｐ，２００５．
ＭＭ
Ｃ，Ｏ’ＨＡＧＡＮＡ．ＢａｙｅｓｉａｎＪ，ＢＥＲＧＥＲＪ
Ｏ，ＰＡＵＬＯ
ｃａｌｉｂｒａｔｉｏｎｏｆｃｏｍｐｕｔｅｒ
ｆ／ｉｔｍｅｗｏｌ＇ｋ
ｍｏｄｅｌ［Ｊ］．Ｊ
ｖａｌｉｄａｔｉｏｎ
ｏｆ
ＲＳｔａｔＳｏｃ：ＳｅｘＢ：Ｓｔａｔ
Ｍｅｔｈｏｄｏｌｏｇｙ，２００１，６３（３）：４２５－４６４．
ＩｎｓｔＳｔａｒ
Ｒ，ｅｔａ／．Ａ
ｆｏｒ
ｃｏｍｐｕｔｅｒ
ｍｏｄｅｌｓ［Ｒ］．ＮｏｒｔｈＣａｒｏｌｉｎａ：Ｎａｔｌ
ＳｃｉＴｅｃｈ
［５］刘信恩，肖世富，莫军．高斯过程响应面法研究［Ｊ］．应用力学学报，２０１０，２７（１）：１９０—１９５．
ＬＩＵＸｉｎ’ｅｎ，ＸＩＡＯＳｈｉｆｕ，ＭＯＪｕｎ．Ａ
ｎｅｗ
ｒｅｓｐｏｎｓｅ
ｓｕｒｆａｃｅｍｅｔｈｏｄｂａｓｅｄｏｎＧａｕｓｓｉａｎ
ｐｒｏｃｅｓｓ［Ｊ］．ＣｈｉｎＪＡｐｐＩＭｅｃｈ，２０１０，２７（１）：１９０—１９５．
［６］
刘信恩，肖世富，莫军．复杂数值模拟的贝叶斯模型确认框架及其简化［ｃ］∥中国力学学会学术大会’２００９，ｃｃＴＡ№００９瑚３７８９。
郑州，２００９．
Ｍ
Ｄ，ＭＩＴＣＨＥＬＬＴ
［７］ＭＯＲＲＩＳＪ．Ｅｘｐｌｏｒａｔｏｒｙｄｅｓｉｇｎｓ
ｆｏｒｃｏｍｐｕｔａｔｉｏｎａｌｅｘｐｅｒｉｍｅｎｔｓ［Ｊ］．Ｊ
Ｓｔａｒ
Ｐｌａｎｎｉｎｇ＆Ｉｎｆｅｒｅｎｃｅ。１９９５，４３（３）：３８１－４０２．
（编辑陈锋杰）
万方数据
ｈｔｔｐ：／／ｗｗｗ．ｅＭｎａｃａｅ．ｃｎ