三维地应力场BP反分析的改进_戴荣

第24卷 第1期
岩石力学与工程学报 V ol.24 No.1
2005年1月 Chinese Journal of Rock Mechanics and Engineering Jan.，2005
收稿日期：2003–05–09；修回日期：2003–07–08
作者简介：戴 荣(1979–)，男，2000年毕业于清华大学水利水电工程系水工结构专业，现为博士研究生，主要从事地下洞室工程的数值模拟和智能分析方法方面的研究工作。

E–mail ：dairong00@ 。

三维地应力场BP 反分析的改进
戴 荣，李仲奎
(清华大学 水利水电工程系，北京 100084)
摘要：从回归分析出发，对地应力场的神经网络反分析进行了改进：采用线弹性有限元计算进行线性回归分析，获得优化参数的大致范围；应用均匀设计来确定计算参数不同水平的组合，进行弹塑性有限元计算获得训练的样本；采用Levenberg-Marquardt 算法来训练BP 神经网络，以提高效率；使用及早停止和正规化方法来避免神经网络的过拟合问题；得到离散化的应力值后，用神经网络拟合出了以坐标为参量的全场应力函数。

关键词：岩石力学；地应力场；BP 神经网络；均匀设计；反分析
中图分类号：TD 311 文献标识码：A 文章编号：1000–6915(2005)01–0083–06
MODIFIED BP BACK ANALYSIS OF 3D IN-SITU STRESSES
DAI Rong ，LI Zhong-kui
(Hydraulic Engineering Department ，Tsinghua University ，Beijing 100084，China )
Abstract ：Based on regression analysis ，an improved method is provided for back propagation (BP) back analysis of in-situ stresses. Linear elastic FEM calculation and linear regression analysis are implemented to determine the general bounds of optimized parameters. Uniform design method is carried out to settle different combinations of factor levels. Training samples are gained by elastoplastic FEM analysis. Levenberg-Marquardt algorithm is applied to achieve better performance during the training process of BP neural networks. Early stopping and regularization are employed to avoid over-fitting problem. When discrete stress values are obtained ，neural networks are also used to construct a global stress function with coordinates as parameters.
Key words ：rock mechanics ；in-situ stresses ；back propagation neural networks ；uniform design ；back analysis
1 引 言
在岩土工程稳定性的数值分析中，初始地应力场是一个非常重要的因素，它是计算分析的初始条件，也是开展进一步分析的基础。

因此，需要对地应力场进行有效的模拟。

地应力场模拟通常有两种方法：一是根据实测的开挖位移进行反分析；二是根据少数实测地应力点的资料，做一定的假定和简化，对初始地应力场进行反分析。

本文将讨论后一
种做法。

根据实测地应力点对初始地应力场进行反分析的常用方法有：边界荷载调整法[1]、应力函数法[2]、位移函数法[3]和有限元多元回归分析[4]等。

边界荷载调整法，是在区域的边界上施加分布形式不同的水平应力和剪切应力，进行多次试算；其缺点是计算结果的收敛性没有保证，计算者的主观随意性比较大。

应力函数法，则采用边界配点的方式，建立一个4次的应力函数来描述域内的2次应力张量场；其缺点是当应力变化剧烈时需要更高阶的应力函
• 84 • 岩石力学与工程学报 2005年
数。

位移函数法与应力函数法类似，它用一个3次的位移函数来描述域内的位移场，进而求得域内的应力场。

有限元多元回归分析，则是建立起初始应力场与待定因素(如自重、构造运动等)之间的多元回归方程，根据计算值与实测应力点的残差平方和达到最小的原则，求得回归方程中各待定因素的系数。

从理论上来说，这种方法是比较合理的。

因此，陆续有不少学者对这一方法进行了研究[5
～8]。

进行多元回归分析时，常用的回归函数为待定因素的线性多项式，这样就忽略了它们之间的互相影响。

如果采用2次以上多项式进行回归，随着待定因素个数的增加和多项式阶数的增长，计算量迅速增加。

而初始应力场和待定因素之间的关系可能比较复杂，用低次多项式函数又是难以描述的。

为了解决这个矛盾，引入了人工神经网络理论。

人工神经网络是由大量的基本信息处理单元(神经元)互相连接而形成的计算体系，它具有以任意精度逼近任何连续函数的能力[9]。

最常用的是BP(back propagation)网络。

在地应力场的反分析问题上，也有学者进行这方面的探索
[10，11]。

其基本思
想是利用神经网络来表达初始应力场和待定因素之间的非线性映射关系。

其基本步骤是：首先以地应力测点的计算值为输入，待定因素为输出，建立神经网络模型；然后通过多次正演计算获得神经网络的训练样本；根据样本对网络进行训练；训练完成后，以地应力测点的实测值作为输入，得到待定因素的输出；最后以此输出参数进行一次正演计算，即可得到初始应力场。

在建立神经网络模型的过程中，需要解决一些问题。

譬如：如何建立神经网络模型，选择哪些待定因素，其取值范围如何，怎么设计训练样本，怎样训练神经网络，如何构造全场应力函数等。

本文在进行三维初始地应力场的BP 神经网络反分析时，对这些方面进行了一些探索和改进，下面介绍一下处理方法。

2 三维地应力场BP 反分析与拟合
在介绍三维地应力场BP 反分析与拟合之前，有必要简单地介绍一下多元线性回归的概念[4]。

2.1 多元线性回归的概念
一般认为，三维初始应力场是如下参变量的函数：
)(z y x U U U G E z y x f ，，，，，，，，，γµσ= (1)
式中：σ为初始应力值；x ，y ，z 为应力点空间位置的坐标；E ，µ，γ分别为岩体的弹性模量、泊松比和容重；G 为自重因素；U x ，U y ，U z 为地质构造位移因素。

其中，x ，y ，z 可通过地质勘测获得；E ，µ，
γ可通过材料试验获得，在对初始应力场进行反分
析时，认为是定值；通常认为，地应力场由自重应力和构造应力2部分组成，因此将G ，U x ，U y ，U z 视为待定因素。

经验表明，垂向地应力基本上由自重决定，而构造作用主要是水平方向，因此在反分析中需考虑在区域边界上的水平构造运动U x 和U z ，必要时也可将垂向构造运动U y 纳入考虑范围。

通常认为构造运动在边界上为均匀分布。

根据待定因素的选择，应力场函数(式(1))可表示为
)(z y x U U U G f ，，，=σ (2)
在弹性材料模型的情况下，在计算域上分别施加重力G 和任意构造位移U x ，U y ，U z 进行正分析，得到基本初始应力场G σ，x σ，y σ，z σ：
⎪⎪⎭
⎪
⎪⎬⎫=
===)000()000()000()000(z z y y x x G U f U f U f G f ，，，，，，，，，，，，σσσσ (3)
然后根据叠加原理，运用最小二乘法，进行回归分析，可以得到
==′)(z z y y x x G U b U b U b G b f ，，，σ
z z y y x x G G b b b b σσσσ+++ (4)
式中：σ′为初始应力的拟合值；G σ为自重应力；
x σ，y σ，z σ为基本初始应力；G b ，x b ，y b ，z b 为
相应的拟合系数。

一般说来，自重应力的计算比较精确，因此其回归系数G b 通常取1.0。

x b ，y b ，z b 的数值等于构造位移U x ，U y ，U z 应缩放的比例。

拟合应力场等于基本应力场的线性叠加。

2.2 神经网络模型的建立
在线弹性回归中，x σ，y σ，z σ之间相互独立，所以可以应用弹性力学的叠加原理，得到共同作用下的初始应力场(式(4))。

而实际上，各个构造作用是同时施加在计算域上的，在采用弹塑性材料模型的情况下，它们的作用互相耦合，也即
≠=′)(z z y y x x G U b U b U b G b f ，，，σ
z z y y x x G G b b b b σσσσ+++ (5)
第24卷 第1期 戴 荣等. 三维地应力场BP 反分析的改进 • 85 •
此时，f 是一个非线性函数，且不能应用叠加原理，因此式(4)不再成立了。

可以看到，在重力保持不变的情况下，拟合的初始应力场σ′取决于U x ，U y ，U z 的大小，即
)(z y x U U U Φf ，，=′σ (6)
式中：Φ为一个复杂的非线性函数。

如果在此引入神经网络模型，利用它能够以任意精度逼近任何连续函数的能力来表达该函数关系，则可以回避叠加原理是否适用的问题。

同样，利用神经网络可以得到式(6)的反函数关系，也即
)(][1T σ′=−ΦU U U z y x (7) 式(7)表述的即是初始应力场的反分析问题。

神经网络的反分析，就是以此为依据进行的。

2.3 参数及其取值范围的确定
利用神经网络进行反分析时，首先要进行多次正分析，以获得训练样本。

从式(6)可见，计算中有3个参数U x ，U y ，U z 可以进行调整。

为了大致确定其取值范围，可先按式(3)，(4)进行一次线弹性回归分析。

经过回归分析可得到
x b ，y b ，z b ，即构造位移应取值为b x U x ，b y U y ， b z U z 。

可以认为，回归分析的结果是问题的一个近似解，构造位移的真实值应在b x U x ，b y U y ，b z U z 周围变动，可以根据经验取一个合适的范围。

2.4 训练样本的设计与获得
在U x ，U y ，U z 的取值范围内，可对它们进行调整。

但参数应在问题空间中分布均匀，才能充分反映问题的特性。

也就是说，要在参数的取值范围内挑选出有代表性的点，这是一个实验设计的问题。

常用的实验设计方法有正交设计和均匀设 计[12]。

正交设计具有“均匀分散，整齐可比”的特点，是科学实验中经常采用的方法。

若一项实验中有s 项因素，每个因素有q 个水平，采用正交实验安排设计的话，至少需做q 2次实验。

当q 较大时，实验次数明显太多。

为了减少实验次数，均匀设计去掉了对“整齐可比”的要求，只考虑实验点在实验范围内的均匀散布。

这样，只需做q 次实验便可满足要求。

正交设计和均匀设计的实验方案均制成了表格，使用时根据需要查表即可，非常方便。

在地应力的反分析中，笔者选择了均匀设计的方案，主要是从计算量上考虑，希望计算的次数越
少越好。

而与正交设计相比，均匀设计所需的计算次数要少得多，但与此同时，获得的样本数也更少，对神经网络的训练，可能会有不利影响。

从这个意义上讲，正交设计得到的多个样本，是有价值的。

具体分析时，应在计算量和样本数之间权衡，选择合适样本的设计方案。

对U x ，U y ，U z 的每一组取值，都进行一次有限元的正向计算，可以获得相应测点的地应力值。

计算结果即作为神经网络的训练样本。

2.5 神经网络的训练[9]
在通过数值计算获得神经网络的训练样本之后，要对神经网络进行训练，以得到测点的应力值和构造位移值之间的一个非线性映射关系。

传统的BP 算法，存在收敛慢、易陷于局部极小等问题，因此人们提出了各种改进的办法。

一般的办法是动量法和自适应调整学习率法。

BP 算法的本质是用最速下降法来极小化误差函数，因此，引入共轭梯度法、拟牛顿法、Levenberg–Marquardt 算法等效果更好的最优化算法，可以提高学习收敛的速度。

笔者采用了Levenberg–Marquardt 算法，网络的权值按如下公式进行调整：
e J I J J T 1T 1)(−++−=µωωk k (8)
式中：1+k ω和k ω为网络的权值，e 为误差向量，J
为误差向量e 对权值ω的Jacobian 矩阵，µ为1个标量。

当µ很小时，变为牛顿法，µ很大则成为最速下降法(此时步长为1−µ)。

因为牛顿法收敛更快，因此每次成功迭代(误差函数减小)之后，就减小µ值；反之则增大µ值。

这样可保证误差函数总是减小的。

另外，神经网络也存在过拟合问题，主要是因为使用了规模过大的网络结构的缘故。

这影响了它的泛化推广能力。

通常可以采用2种办法解决：及早停止(early stopping)和正规化(regularization)方法。

及早停止方法将训练数据集分为3个部分：训练样本、验证样本和测试样本。

使用训练样本进行学习，同时监测验证样本的误差。

验证样本的误差在网络训练的初期是逐渐减小的，而当网络出现过拟合时，误差开始上升。

此时就应停止训练。

测试样本用于比较不同的网络模型，它的误差开始增加的时间应接近于验证样本。

正规化方法对神经网络的误差函数进行了修改，在其中加入了网络权值的影响，即
∑∑==−+−=N
i n j j i i n x f y N
E 1
12
2
1)1())((1ωαα
(9)
• 86 • 岩石力学与工程学报 2005年
式中：E 为误差函数，α为[0，1]之间的系数，N 为样本数，y i 为目标值，x i 为神经网络的输入，)(i x f 为神经网络的输出，n 为神经网络的权值个数，j ω为权值。

式(9)中右边的第1项是误差均方和，即为一般的误差函数；右边的第2项是权值的均方和。

这样会强迫网络取较小的权值，从而使拟合的函数更加平滑。

一般说来，用正规化方法得到的函数，比用及早停止方法得到的函数平滑性要好。

本文对这两种方法的效果进行了比较。

神经网络的训练完成之后，得到了测点的应力值和构造位移之间的映射关系。

此时以地应力测点的实测值作为输入，输出即为构造位移U x ，U y ，U z 的真实值。

根据它进行一次有限元的正演计算，
即可得到初始应力场。

2.6 全场应力函数的拟合
通过反分析，获得了求解区域中整体的应力分布情况。

但是现在得到的应力场是用离散的数值表示的，而且依赖于单元划分。

如果要查询区域中某点的应力状况，必须找到该点处的单元，然后进行插值。

如果能用一个以坐标为参量的函数来描述整个应力场，对于使用者而言，将会非常方便。

即需构造函数：
)(z y x f ，，=σ (10)
式中：σ为地应力值；x ，y ，z 为应力点的坐标。

式(10)是一个比较复杂的向量函数，可以采用正交多项式的方法构造[6
，7]
，但是对神经网络而言，
其对数据的出色拟合能力，正可在这个问题的求解上发挥作用。

只要提取一部分有代表性的节点的计算结果作为训练样本，以应力点的坐标x ，y ，z 为输入，地应力值σ为输出，进行训练即可得到全场应力函数。

为了评价神经网络的拟合程度，可以定义地应
力场的拟合值σ′与实测值σ之间的相关系数r ，得
2
2
)
()())((i i i i i i i i r σσσσσσσσ′−′Σ−Σ′−′−Σ=
(11)
式中：i σ′为地应力的拟合值，i σ′为其均值，i σ为地应力的实测值，i σ为其均值。

上面介绍了神经网络进行地应力反分析和拟合的基本方法，下面看一个具体的例子。

3 应用实例
某水电站地下厂房位于基岩坡体内。

为了了解
山体中的地应力分布情况，在平硐PD2及PD5中分别进行了3个方向钻孔水压致裂三维应力测量，共获得了4个地应力测点的应力值。

最大主应力方向为N55°～68°W/SE ∠33°～35°。

为了便于处理，将应力值转换到计算坐标系中，见表1。

表1 测点的实测应力值
Table 1 In-situ measured stress values MPa
测点
σx σy σz τxy τyz τxz 1 －6.88 －5.32－5.94 1.44 －2.52 1.44 2
－10.97
－8.00－9.09 2.35 －4.48 2.32 3 －7.38 －7.61－7.37 0.12 －2.32 2.28 4 －8.99
－10.01
－9.87 0.51 －2.71
2.45
为了进行有限元正向分析，将整个计算区域离散为9 845个单元，14 573个节点，如图1所示。

采用ANSYS 进行有限元计算分析。

图1 有限元网格图 Fig.1 FEM mesh
用神经网络的方法进行地应力反分析时，如前所述，首先进行一次线性回归分析，以大致确定构造位移U x ，U y ，U z 的取值范围。

回归分析的结果：U x = 0.255 6，U y = 0.399 6，U z = 0.207 0，r = 0.91。

回归分析得到的地应力拟合值如表2所示。

表2 测点的回归地应力拟合值
Table 2 In –situ stress values by regression MPa
测点σx σy σz τxy τyz τxz 1 －7.57－7.08－6.89 0.83 －0.95－0.522 －7.81－7.30－7.27 0.87 －0.81－0.433 －9.94－8.46－9.34 0.99 －0.37－0.264
－8.80
－7.88
－8.13 0.48 －0.42
－0.29
根据经验，将U x ，U y ，U z 的取值范围确定为：
第24卷第1期戴荣等. 三维地应力场BP反分析的改进• 87 • U x∈[0.1，0.5]，U y∈[0.2，0.8]，U z∈[0.1，0.5]。

通过调整U x，U y，U z的数值，进行正向计算，
可以获得神经网络的训练样本。

训练样本的数值调
整采用均匀设计的思想。

选择U9(9d)表进行设计，
每个因素(构造位移)取9个水平。

各因素水平的组
合如表3所示。

表3 各因素水平组合
Table 3 Combinations of factor levels
水平U x /m 水平U y/m水平U z/m
1 0.10 6 0.575 3 0.20
2 0.15 2 0.275 8 0.45
3 0.20 9 0.800 6 0.35
4 0.2
5 3 0.350 1 0.10
5 0.30 5 0.500 5 0.30
6 0.35
7 0.650 9 0.50
7 0.40 1 0.200 4 0.25
8 0.45 8 0.725 2 0.15
9 0.50 4 0.425 7 0.40
首先按照表3中各因素水平的组合，按弹塑性
材料进行计算。

通常情况下，区域中的应力状况会
与使用弹性材料计算有所不同，但变化幅值不应太大。

将这9个计算结果作为样本，用BP神经网络
进行训练。

以这4个测点共24个地应力分量的计算
值作为输入，3个构造位移U x，U y，U z作为输出。

网络结构选择3层网络，隐含层有8个单元，网络
的结构为24–8–3，采用Levenberg–Marquardt算
法进行训练。

训练的结果与样本的相关系数均达到
了1.0。

可见训练结果是相当不错的。

当以实测应力值作为训练好的神经网络的输入时，网络的输出结果为U x = 0.232 7，U y = 0.365 7，
U z = 0.264 0，这就是所对应的3个构造位移值。

以
这3个参数作为计算区域的边界条件，再进行一次
正向计算。

最终反分析得到的各实测点的地应力拟
合值如表4所示，相关系数r = 0.89。

虽然实测地应力数值的规律性不强，但采用了
弹塑性材料模型之后，神经网络反分析得到的地应
力拟合值与实测值相比，依然能达到0.89左右的相
关度。

由于计算量的限制，样本数比较少，分析的
结果受到一定程度的影响。

如果每个因素多取几个
水平，多进行几次计算以增加样本的数目，那么计
算结果的相关度还可以提高。

表4 测点的神经网络地应力拟合值
Table 4 In-situ stress values by neural networks MPa 测点σxσyσzτxyτyzτxz
1
－6.45－6.15－7.74 0.40 －0.660.05
2
－8.15－7.76－8.79 0.75 －0.24－1.19
3 －10.09－8.38－10.6
4 1.8
5 0.22－0.33
4
－8.25－7.60－9.26 0.37 －0.97－1.04从计算结果中选择其中853个节点(约5.85%)的应力值，将节点的3个坐标作为神经网络的输入，而节点的6个应力分量作为神经网络的输出，可以建立坐标和应力分量的映射关系。

注意：为了保证选取的样本点的代表性，应尽量保证节点在整个空间中的分布比较均匀。

选择最常用的BP神经网络进行拟合。

网络结构选择3层网络，因为数据量较大，选择隐含层有14个单元，网络的结构为3–14–6。

为了防止神经网络的过拟合问题，提高其泛化能力，分别采用及早停止和正规化方法对网络进行训练。

采用Levenberg-Marquardt算法进行训练，以提高效率。

拟合结果如表5所示。

表5 不同训练方法的相关系数
Table 5 Correlation coefficients of different training
methods
相关系数
训练方法
σxσyσzτxyτyzτxz
及早停止0.919 0.970 0.921 0.815 0.7590.731正规化0.939 0.973 0.940 0.867 0.8060.755由表5可见，2种训练方法得到的结果是接近的，正规化的方法结果稍好一点。

在正应力的拟合上，两者都能达到0.9以上的相关度，剪应力的拟合也能达到0.7以上的相关度。

通常剪应力的变化规律比正应力更加难以描述，拟合结果恰好说明了这一点。

从拟合结果来看，用神经网络来拟合全局应力函数还是有一定精度的。

4 结论
本文从初始地应力场的回归拟合出发，应用神经网络有关理论，对三维地应力场进行了反分析。

• 88 • 岩石力学与工程学报 2005年
对BP反分析中的问题，提出了一些改进：
(1) 利用线弹性有限元进行线性回归分析，获
得设计参数的大致范围；
(2) 应用均匀设计来确定计算参数不同水平的
组合，并用弹塑性有限元数值计算获得训练的样
本；
(3) 采用Levenberg-Marquardt算法来训练神经
网络，以提高效率；
(4) 使用及早停止和正规化方法来避免神经网
络的过拟合问题；
(5) 得到离散化的应力值之后，用神经网络拟
合出了以坐标为参量的全场应力函数。

实例计算结果表明，用BP神经网络对地应力
场进行反分析，不必进行繁琐的参数敏感性分析，
具有操作简便、程式化强的特点，并且反分析结果
也具有较好的精度。

参考文献(References)：
[1] 丰定祥，谷先荣，杨家岭等. 关于地下工程有限元分析中初始地
应力场的假定[J]. 地下工程，1982，(2)：20–27.(Feng Dingxiang，
Gu Xianrong，Yang Jialing，et al. Assumptions on initial stress field
in FEM analysis of underground engineering[J]. Underground
Engineering，1982，(2)：20–27.(in Chinese))
[2] 张有天，胡惠昌. 地应力场的趋势分析[J]. 水利学报，1984，(4)：
31–38.(Zhang Youtian，Hu Huichang. Trend analysis of residual
stress distribution in rock mass[J]. Journal of Hydraulic Engineering，
1984，(4)：31–38. (in Chinese))
[3] 杨柯，张立翔，李仲奎. 地下洞室群有限元分析的地应力场计
算方法[J]. 岩石力学与工程学报，2002，21(11)：1 639–
1 644.(Yang Ke，Zhang Lixiang，Li Zhongkui. New method for
calculating geostresses in FEM analysis of underground houses[J].
Chinese Journal of Rock Mechanics and Engineering，2002，21(11)：
1 639–1 644.(in Chinese))
[4] 郭怀志，马启智，薛玺成等. 岩体初始应力场的分析方法[J]. 岩
土工程学报，1983，5(3)：64–75.(Guo Huizhi，Ma Qizhi，Xue Xichen，et al. The analytical method of the initial stress field for rock mass[J]. Chinese Journal of Geotechnical Engineering，1983，5(3)：64–75.(in Chinese))
[5] 朱焕春，赵海滨. 河谷地应力场的数值模拟[J]. 水利学报，1996，
(5)：29–36.(Zhu Huanchun，Zhao Haibin. Numerical simulation for
rock stress field in a valley[J]. Journal of Hydraulic Engineering，1996，
(5)：29–36.(in Chinese))
[6] 李青麟. 初始应力的回归与三维拟合[J]. 岩土工程学报，1998，
20(5)：68–71.(Li Qinglin. Regression analysis and 3D fitting of initial stress[J]. Chinese Journal of Geotechnical Engineering，1998，20(5)：68～71.(in Chinese))
[7] 肖明，刘志明. 锦屏二级水电站三维地应力场反演回归分析[J].
人民长江，2000，31(9)：42–44.(Xiao Ming，Liu Zhiming. Feedback regression analysis of 3D geostress field in Jinping II hydropower station[J]. Yangtze River，2000，31(9)：42–44.(in Chinese)) [8] 邵国建. 岩体初始地应力场的反演回归分析[J]. 水利水电科技进
展，2002，20(5)：36–38.(Shao Guojian. Feedback regression analysis of initial stress field for rock masses[J]. Advances in Science and Technology of Water Resources，2002，20(5)：36–38.(in Chinese)) [9] 阎平凡，张长水. 人工神经网络与模拟进化计算[M]. 北京：清华
大学出版社，2000.(Yan Pingfan，Zhang Changshui. Artificial Neural Networks and Simulated Evolutionary Computation[M]. Beijing：
Tsinghua University Press，2000.(in Chinese))
[10] 蒋中明，徐卫亚，邵建富. 基于人工神经网络的初始地应力场三
维反分析[J]. 河海大学学报，2002，30(3)：52–56.(Jiang Zhongming，Xu Weiya，Shao Jianfu. ANN based 3D back analysis of initial stress in rock masses[J]. Journal of Hohai University，2002，
30(3)：52–56.(in Chinese))
[11] 蒋中明，徐卫亚. 三维初始地应力场反分析的径向基函数法[J].
岩土力学，2002，23(6)：737–741(Jiang Zhongming，Xu Weiya. 3D back analysis of initial ground stress based on ANN–RBF[J]. Rock and Soil Mechanics，2002，23(6)：737–741.(in Chinese))
[12] 方开泰. 均匀设计与均匀设计表[M]. 北京：科学出版社，
1994.(Fang Kaitai. Uniform Design and Design Tables[M]. Beijing：
Science Press，1994.(in Chinese))。