自20世纪80年代引入产品退化系统的概念以来[2]，国内外****对利用部件退化概率计算系统可靠性的问题提出了很多解决方案[3-5]。一些方法着重分析关键节点结构(如串联、并联、k/n系统等)组成的系统形式，如Li等[6]假设系统部件为受权重因素影响的线性退化路径，根据该退化路径公式推导出了基于部件失效概率密度函数和系统可靠性函数的关系数学模型，由此来评估多部件组成系统的可靠性；Faghih-Roohi等[7]结合通用生成函数(Universal Generating Function, UGF)和马尔可夫过程，建立了多状态加权k/n系统可靠性评估的动态模型，并利用遗传算法针对系统各部件在不同状态下的概率随时间变化问题进行优化。另一些方法对组成系统的关键组成部件进行分析，进而以关键部件的退化形式估计系统整体的退化形式，如Li等[8]通过构建退化模型将二元状态可靠性分析推广到连续状态可靠性分析，并利用状态故障树对系统的退化和灾难性失效进行可靠性分析；彭宝华等[9]从元器件性能与系统性能的关系出发，利用元器件的性能退化数据研究系统性能可靠性，给出了系统可靠度计算的马尔可夫链蒙特卡罗(Markov Chain Monte Carlo，MCMC)算法和近似解析方法。虽然上述方法均对部件退化与系统退化的关系进行了深入探讨，但在考虑建立复杂系统结构的基础上，对系统整体进行可靠性分析的成果仍然较少。
本文旨在原有退化系统分析方法基础上建立系统结构模型，提出对整个退化系统进行可靠性分析的有效方法。首先，考虑对组成系统的各部件退化性能进行分析，并建立相应的可靠性模型。对于组成EWIS的高可靠性部件来说，其某些性能指标随着时间而退化，统计和分析这些性能指标退化量，可推测出产品的可靠水平[10]。如汇流条阻抗在使用过程中有变大的趋势，当阻抗退化量达到失效阈值时将发生失效。因此，可以采用退化数据评估其可靠性。其次，考虑根据系统结构建立系统整体的退化可靠性模型。GO法(Goal Oriented methodology)是分析复杂系统可靠性的一种有效方法，且特别适用于实际物流(电流、气流、液流等)构成的多状态复杂系统可靠性分析[11]，因此可采用GO法来描述EWIS的结构关系。最后，考虑结合部件退化可靠性的计算结果分析系统整体退化可靠性问题。UGF法是多状态系统分析领域中的重要方法，结合多状态系统的结构函数，可有效反映系统多状态性能随部件多状态性能变化的关系，且能够降低系统可靠性计算的复杂度[3]。因此结合GO法与UGF的UGF-GO法能较好达到利用部件性能退化分析系统整体退化可靠性的目的。
本文基于上述思路，提出了基于性能退化和UGF-GO法的EWIS退化系统可靠性分析方法，即结合部件性能参数退化可靠性特征，建立EWIS退化状态可靠性分析模型，利用UGF-GO法对系统退化可靠性进行评估，并将该方法应用于某飞机电气线路互联系统，计算结果表明该方法可有效解决系统退化状态可靠性分析问题。
1 退化可靠性建模及参数估计 1.1 基于Wiener过程的汇流条性能退化模型 EWIS中汇流条性能退化一般由微小接触电阻阻值增长累积形成的性能退化表征，可采用如下一元Wiener过程建模[12]：
 | (1) |
式中：β和σ分别为退化速率和扩散参数；Λ(t)为时间t的单调递增函数，ΔΛ(t)=Λ(t+Δt)-Λ(t)为观察时间增量，且满足Λ(0)=0；W(Λ(t))为标准布朗运动；Y(t)具有平稳独立增量，且增量服从正态分布，即Δy(t)~N(βΔΛ(t), σ2ΔΛ(t))。
设Λ(t)=t，易知一元Wiener过程是关于时间的线性函数，可对线性退化过程建模。随机抽取n个产品样本进行m次退化观测，令Y表示样本n，Y=(Y1, Y2, …, Yn)，Yi表示i条退化轨迹，Yi=(yi0, yi1, …, yim)，设Y(t)首次达到失效阈值l时产品失效，则产品寿命T关于首达时间t服从逆Gaussian分布[13]，其概率密度函数和可靠性函数分别为
 | (2) |
 | (3) |
式中：Φ()为正态分布函数。
同类产品个体的退化过程不可能完全一致，即使同型号产品，由于各产品自身性能及工作环境的不同，其退化程度也不完全相同。为描述这种差异性，设式(3)中退化速率β~N(μβ，σβ2)为随机变量，相应的失效函数和可靠性函数为
 | (4) |
 | (5) |
式中：μβ、σβ2和σ2为待估参数。
易知在确定退化参数和失效阈值的情况下，可得到任意时刻含个体差异的退化产品可靠性。
1.2 基于MCMC算法的模型参数估计 根据Bayes统计推断，设参数μβ、σβ2和σ2统计独立且服从给定的共轭先验分布，μβ~N(θμ, θσ2)，σβ2~Ga(θ1, γ1)，σ2~Ga(θ2, γ2)，其中，伽马分布可以表示较大范围的密度函数特征，包括减函数、常数及增函数，能够较全面地描述分布特征的灵活性[14-15]。Θ=(μβ, σβ2, σ2)的联合先验分布为
 | (6) |
只考虑分布形式的情况下，式(6)可简化为
 | (7) |
令Ω=(θμ, θσ, θ1, θ2, γ1, γ2)为先验分布的初始设定参数。
根据Wiener过程独立增量服从正态分布的性质，得到退化数据增量Δy(t)关于参数μβ、σβ2和σ2的最大似然函数为
 | (8) |
式中：Δyij是第i条退化轨迹第j次观测数据增量; Δtij为相应时间增量，根据Bayes公式推导Θ的后验概率分布形式为
 | (9) |
式中：分母看成是一个标量因子，代入式(7)和式(8)，简化后得到关于参数μβ、σβ2和σ2的后验概率分布为[16]
 | (10) |
式(10)为关于参数μβ、σβ2和σ2的先验概率分布和产品退化过程观测数据增量的后验分布，不易直接计算参数μβ、σβ2和σ2值，采用MCMC算法估计[17]。该算法将2个或多个随机变量的条件概率分布分解为模型中每个参数的条件概率分布进行采样，从中生成一个参数样本序列，MCMC算法迭代算式为
 | (11) |
式中：k为当前随机参数值; k+1为迭代后更新参数值。各参数值迭代计算，具体算法步骤如下：
步骤1?设置初值μβ(0)、σβ2(0)和σ2(0)，迭代计数k=0。
步骤2?从N(μβ(k), 1)中随机获得采样点μβ*，计算
 |
步骤3?从U(0, 1)中随机抽取u1，若u1≤α1，则令μβ(k+1)=μβ*，否则μβ(k+1)=μβ(k)。
步骤4?从N(σβ2(k), 1)中随机获得采样点σβ2*，计算
 |
步骤5?从U(0, 1)中随机抽取u2，若u2 ≤ α2，则令σβ2(k+1)=σβ2*，否则σβ2(k+1)=σβ2(k)。
步骤6?从N(σ2(k), 1)中随机获得采样点σ2*，计算
 |
步骤7?从U(0, 1)中随机抽取u3，若u3≤α3，则令σ2(k)=σ2*，否则σ2(k+1)=σ2(k)，k=k+1。
步骤8 ?k=5 000迭代终止，计算各参数采样序列的期望，作为参数μβ、σβ2和σ2估计值。
1.3 汇流条可靠性模型及参数估计 用随机Wiener过程仿真某飞机EWIS关键部件性能退化。设主设备汇流条失效阈值为5 mΩ，同型号10组部件样本经10次观测，其退化数据如图 1所示。
 |
| 图 1 主设备汇流条性能退化仿真数据 Fig. 1 Performance degradation simulation data of main equipment busbar |
| 图选项 |
利用仿真数据按式(3)建立主汇流条退化可靠性模型，得失效阈值为5 mΩ时，可靠性计算公式为
 | (12) |
式中：未知参数μβ、σβ2和σ2先验概率分布可以通过式(6)和式(7)得到，其中参数Ω=(θμ, θσ, θ1, θ2, γ1, γ2)通过经验确定，本文设置Ω=(2, 1, 2, 2, 3, 20)，关于未知参数μβ、σβ2和σ2的后验概率分布形式可由式(9)确定。MCMC算法迭代结果与参数初始值无关，初值可随机给定，设置Θ=(0, 0, 0)，利用MCMC算法对参数进行估计，经5 000次采样，得到各参数样本序列，文中将各参数样本序列期望作为参数估计值，参数的期望、方差及置信区间等估计量见表 1。
表 1 主设备汇流条可靠性模型参数估计 Table 1 Parameter estimate for reliability model of main equipment busbar
| 参数 | 期望 | 方差 | 置信区间(置信水平为0.95) |
| μβ | 0.447 0 | 0.017 7 | [0.340 0，0.545 6] |
| σβ2 | 0.005 0 | 0.009 6 | [0.000 3，0.007 6] |
| σ2 | 0.222 3 | 0.371 6 | [0.129 0，0.301 0] |
表选项
2 系统可靠性估计 2.1 GO法系统建模 GO法是从系统原理图分析系统成功运行的概率，用操作符和信号流表示系统各单元及逻辑关系建立GO图，进行可靠性分析[18-19]。图 2为根据电子设备系统建立的GO图。
 |
| 图 2 电子设备系统GO图 Fig. 2 GO figure of electronic equipment system |
| 图选项 |
2.2 UGF-GO法基本运算方法 UGF法是按系统的结构形式，进行递归计算求出系统可靠性的常用方法，已知系统部件状态概率以及系统结构情况下，建立u函数来描述部件状态组合[20-21]，t时刻部件X的u函数为
 | (13) |
式中：xi, j为i部件具有的j种状态；p为X处于各状态的概率；z是x的z变换。在式中引入描述GO图系统结构的 ?g算子，则系统的U函数定义为
 | (14) |
UGF-GO法在计算思想上与GO法定量分析时常用的状态组合方法类似，但在形式上通过UGF改进的GO法更为直观，便于计算。本文以某飞机EWIS中电子设备系统二状态系统为例说明UGF-GO法的基本运算过程。系统有1个输入操作符，4个功能操作符，1个或门操作符，其成功和故障状态概率数据如表 2所示。
表 2 电子设备系统操作符数据 Table 2 Operator data of electronic equipment system
| 操作符编号 | 操作符类型 | 状态概率 | |
| pi, 0 | pi, 1 | ||
| 1 | 5 | 0.01 | 0.99 |
| 2 | 1 | 0.05 | 0.95 |
| 3 | 1 | 0.10 | 0.90 |
| 4 | 1 | 0.05 | 0.95 |
| 5 | 1 | 0.10 | 0.90 |
表选项
据UGF-GO法，各操作符u函数可表示为
 | (15) |
将各操作符的u函数代入GO图进行可靠性系统分析，根据GO图的系统结构得到系统 ?g算子，则系统U函数可表示为
 | (16) |
系统GO图形式对称，u3(z)和u5(z)具有相同的u函数。系统GO图输出信号3、5未计算共有信号1，因此在计算输出信号6时进行共有信号修正。设共有信号故障状态为f，状态概率为v，设含j个共有信号的子系统u函数uj(z)之间相互独立，则子系统输出信号修正运算符ξ定义为
 | (17) |
根据u函数计算规则进行运算。
 |
输出信号6状态概率即为系统的状态概率，根据运算过程得到系统状态值1，即成功状态的可靠性概率为0.969 185。该结果与文献[10]中概率公式计算方法结果一致，验证了UGF-GO法用于GO图定量计算的有效性，部分操作符输出状态概率如表 3所示，P为对应操作符的输出概率。
表 3 电子设备系统可靠性结果对比 Table 3 Reliability result comparison of electronic equipment system
| 输出信号 | 概率公式算法 | UGF-GO法 | |||
| pi, 0 | pi, 1 | pi, 0 | pi, 1 | ||
| Po, 3 | 0.846 450 | 0.153 550 | 0.855 000 | 0.145 000 | |
| Po, 5 | 0.846 450 | 0.153 550 | 0.855 000 | 0.145 000 | |
| Po, 6(未修正) | 0.976 422 | 0.023 577 | 0.978 975 | 0.021 025 | |
| Po, 6(修正后) | 0.969 185 | 0.030 814 | 0.969 185 | 0.030 815 | |
表选项
从上例可见，在计算形式上，利用UGF-GO方法本身的计算规则可以在不列出全部状态组合表的情况下，直接计算GO图各信号流的状态概率，因此本文在分析复杂系统退化可靠性时也采用UGF-GO法。进一步将UGF-GO法应用于含部件退化问题的退化系统可靠性评估中。退化模型算子δ定义为给定阈值l下部件各状态概率的和，根据文献[13]中的描述，该概率状态u函数可表示为在普通u函数中引入δ算子，即
 | (18) |
则系统部件X在失效阈值l下退化可靠性计算模型的u函数可表示为
 | (19) |
式中：Ri(t)为系统i部件随时间变化的可靠性。
3 EWIS退化可靠性分析 某型号飞机电气线路互联系统配置方案如图 3所示，方案中汇流条配置包括主设备汇流条、电子设备汇流条、非必需设备汇流条、蓄电池设备汇流条和连接汇流条等。主设备汇流条连接飞机飞行任务关键用电设备线路，如发动机、起落架等。非必需设备汇流条连接一般用电设备线路，如着陆灯、备用灯等。电子设备汇流条连接两组无线电通信设备线路。蓄电池设备汇流条连接外部供电电源线路与内部蓄电池供电线路[22]。
 |
| 图 3 某飞机电气线路互联系统设计方案 Fig. 3 Design scheme of an aircraft's EWIS |
| 图选项 |
3.1 汇流条退化可靠性 利用Wiener退化过程仿真汇流条接触电阻退化过程，建立EWIS系统各汇流条可靠性模型，分别采用二步法[21]和MCMC算法对EWIS各汇流条可靠性模型参数进行对比估计，得μβ、σβ2和σ2如表 4所示。
表 4 二步法与MCMC算法参数估计值对比 Table 4 Comparison of parameter estimate between two-step algorithm and MCMC algorithm
| 汇流条 | 二步法 | MCMC算法 | |||||
| μβ | σβ2 | σ2 | μβ | σβ2 | σ2 | ||
| 主设备 | 0.441 4 | 0.011 9 | 0.278 6 | 0.447 0 | 0.005 0 | 0.222 3 | |
| 电子设备 | 0.431 1 | 0.014 1 | 0.319 7 | 0.421 9 | 0.005 1 | 0.221 7 | |
| 非必需设备 | 0.473 9 | 0.005 5 | 0.207 0 | 0.464 1 | 0.004 5 | 0.201 0 | |
| 蓄电池设备 | 0.499 8 | 0.005 1 | 0.224 8 | 0.479 3 | 0.005 0 | 0.218 9 | |
表选项
表 4可见，2种方法得到的估计值接近，但MCMC算法估计值与二步法相比波动较小，这是由于在MCMC算法预定义先验概率类型，使估计值范围缩小。仿真得到系统各汇流条可靠性曲线如图 4所示。
 |
| 图 4 四种汇流条退化可靠性曲线 Fig. 4 Degraded reliability curve of four kinds of busbar |
| 图选项 |
由图 4可见，含随机退化速率变量的Wiener过程能够更好地描述同型号产品不同工作环境下的差异性，因此，考虑部件退化状态的系统可靠性更符合实际情况。
3.2 EWIS可靠性建模分析 分析图 2中某飞机EWIS结构，建立描述该系统结构GO图，如图 5所示。GO图中EWIS各操作符及状态概率见表 5。
 |
| 图 5 某飞机电气线路互联系统GO图 Fig. 5 GO figure of an aircraft's EWIS |
| 图选项 |
表 5 某飞机电气线路互联系统操作符数据 Table 5 Operator data of an aircraft's EWIS
| 编号 | 类型 | 名称 | 状态概率 | ||
| pi, 0 | pi, 1 | pi, 2 | |||
| 1 | 5 | 电源 | 0.999 984 | - | 0.000 016 |
| 2~6 | 1 | 导通控制 | 0.999 995 | - | 0.000 005 |
| 7 | 1 | 主设备汇流条 | - | R1(t) | 1-R1(t) |
| 8, 9 | 1 | 电子设备汇流条 | - | R2(t) | 1-R2(t) |
| 10 | 2 | 或门 | - | - | - |
| 11 | 1 | 非必需设备汇流条 | - | R3(t) | 1-R3(t) |
| 12 | 1 | 内部供电开关 | 0.999 995 | - | 0.000 005 |
| 13 | 1 | 外部供电开关 | 0.999 989 | - | 0.000 011 |
| 14 | 2 | 或门 | - | - | - |
| 15 | 1 | 蓄电池设备汇流条 | - | R4(t) | 1-R4(t) |
| 16 | 10 | 与门 | - | - | - |
表选项
表 5中系统部件状态值pi, j，(j=0, 1, 2)分别表示部件处于完好、退化、故障3种状态。电源信号通断及输入、输出控制只考虑有信号通(完好)和无信号断(故障)2种状态，即0、2状态。根据UGF-GO法系统在失效阈值为5 mΩ时可靠性的u函数为
 | (20) |
将表 5中各部件状态概率代入式(20)中计算系统随时间变化可靠性为
 | (21) |
输出信号10、14分别含有共有信号1、6，在计算过程中已经修正。系统输出信号16计算时未考虑共有信号1的影响，需在信号流16输出端进行修正，修正后系统输出状态u函数为
 | (22) |
根据UGF-GO法，式(22)中U函数在1状态的可靠性即为系统退化状态可靠性。为比较不同失效阈值下可靠性曲线的不同变化，选取10、15和20 mΩ为系统失效阈值，得到不同阈值下系统退化可靠性曲线如图 6所示。
 |
| 图 6 不同阈值下的系统可靠性曲线 Fig. 6 System reliability curves under different thresholds |
| 图选项 |
由图 6可知，选取不同失效阈值时，系统可靠性评估结果差异较大。系统可靠性随运行时间可靠性不同程度降低。失效阈值选取越大，则利用当前仿真数据进行的系统可靠性评估结果越可靠，这也较符合实际情况。不同失效阈值下，系统部分可靠性随时间-阈值变化数据如表 6所示。
表 6 不同阈值系统可靠性结果对比 Table 6 Comparison of system reliability results under different thresholds
| 时间/ (103h) | 失效阈值 | |||
| 5 mΩ | 10 mΩ | 15 mΩ | 20 mΩ | |
| 0.5 | 0.999 869 23 | 0.999 869 23 | 0.999 911 11 | 0.999 952 99 |
| 1.0 | 0.983 131 70 | 0.983 201 29 | 0.991 532 16 | 0.999 933 62 |
| 1.5 | 0.807 300 25 | 0.814 910 50 | 0.900 868 99 | 0.995 894 56 |
| 2.0 | 0.410 376 10 | 0.455 760 56 | 0.650 608 32 | 0.928 757 82 |
| 2.5 | 0.111 053 74 | 0.161 576 30 | 0.328 864 11 | 0.669 353 15 |
| 3.0 | 0.017 385 46 | 0.037 127 48 | 0.107 027 91 | 0.308 530 88 |
| 3.5 | 0.001 890 51 | 0.005 873 75 | 0.022 638 82 | 0.087 255 26 |
表选项
4 结论 1) 某飞机EWIS退化严重，且受环境影响较大，采用随机Wiener过程建立部件退化可靠性模型，由MCMC算法估计模型参数，可得到随部件性能退化，复杂系统可靠性变化情况。
2) 在结合GO法及UGF法优势的基础上，提出UGF-GO法。分别利用UGF-GO法与概率公式法计算电子设备电气线路系统GO图，最终得到的系统可靠性结果一致，验证了UGF-GO法计算GO图的正确性和简便性，说明UGF-GO法适用于复杂系统可靠性计算。
参考文献
| [1] | WANG W, WANG P. Effect of environments on EWIS failure rate estimate[J]. Procedia Engineering, 2011, 17: 473-479. DOI:10.1016/j.proeng.2011.10.057 |
| [2] | MEEKER W Q, HAMADA M. Statistical tools for the rapid development and evaluation of high-reliability products[J]. IEEE Transactions on Reliability, 1995, 44(2): 187-198. DOI:10.1109/24.387370 |
| [3] | LEVITEN G. The universal generating function in reliability analysis and optimization[M]. Berlin: Springer, 2005. |
| [4] | LI W J, PHAM H. Reliability modeling of multi-state degraded systems with multi-competing failures and random shocks[J]. IEEE Transactions on Reliability, 2005, 54(2): 297-303. DOI:10.1109/TR.2005.847278 |
| [5] | WANG G, DUAN F, ZHOU Y. Reliability evaluation of multi-state series systems with performance sharing[J]. Reliability Engineering & System Safety, 2018, 173: 58-63. |
| [6] | LI J, COIT D W, ELSAYED E A.Reliability modeling of a series system with correlated or dependent component degradation processes[C]//International Conference on Quality, Reliability, Risk, Maintenance, and Safety Engineering.Piscataway, NJ: IEEE Press, 2011: 388-393. |
| [7] | FAGHIH-ROOHI S F, MIN X, KIEN M G, et al. Dynamic availability assessment and optimal component design of multi-state weighted k-out-of-n systems[J]. Reliability Engineering & System Safety, 2014, 123: 57-62. |
| [8] | LI X, JIANG T M, MA J, et al.State tree analysis of FOG based on drift brownian motion[C]//International Conference on Reliability, Maintainability and Safety.Piscataway, NJ: IEEE Press, 2009: 1322-1326. |
| [9] | 彭宝华, 周经伦, 刘学敏. 元器件性能退化信息的系统性能可靠性评估[J]. 火力与指挥控制, 2011, 36(10): 148-151. PENG B H, ZHOU J L, LIU X M. System performance reliability assessment based on degradation data from components[J]. Fire Control & Command Control, 2011, 36(10): 148-151. DOI:10.3969/j.issn.1002-0640.2011.10.039 (in Chinese) |
| [10] | NIKULIN M S, LIMNIOS N, BALARISHNAN N, et al. Advances in degradation modeling[M]. Boston: Birkh?user Boston, 2010. |
| [11] | 沈祖培, 黄祥瑞. GO法原理及应用:一种系统可靠性分析方法[M]. 北京: 清华大学出版社, 2004. SHEN Z P, HUANG X R. Principle and application of GO methodology:A system reliability analysis methodology[M]. Beijing: Tsinghua University Press, 2004. (in Chinese) |
| [12] | YE Z S, XIE M. Stochastic modelling and analysis of degradation for highly reliable products[J]. Applied Stochastic Models in Business & Industry, 2015, 31(1): 16-32. |
| [13] | SHERIF Y.S, SMITH M L. First-passage time distribution of brownian motion as a reliability model[J]. IEEE Transactions on Reliability, 2009, R-29(5): 425-426. |
| [14] | HAO H B, LI C P. Reliability assessment and residual life prediction method based on Wiener process and current degradation quantity[J]. Engineering Letters, 2016, 24(2): 172-177. |
| [15] | LI C P, HAO H B. Degradation data analysis using Wiener process and MCMC approach[J]. Engineering Letters, 2017, 25(3): 234-238. |
| [16] | FAN T H, CHEN C H. A Bayesian predictive analysis of step-stress accelerated tests in gamma degradation-based processes[J]. Quality & Reliability Engineering International, 2017, 33(7): 1417-1424. |
| [17] | JAIN V. Computational statistics handbook with MATLAB[J]. Journal of the Royal Statistial Society:Series A(Statistics in Society), 2009, 172(4): 942-943. |
| [18] | 沈祖培, 高佳. GO法原理和改进的定量分析方法[J]. 清华大学学报, 1999, 39(6): 15-19. SHEN Z P, GAO J. GO methodology and improved quantification of system reliability[J]. Journal of Tsinghua University, 1999, 39(6): 15-19. DOI:10.3321/j.issn:1000-0054.1999.06.005 (in Chinese) |
| [19] | REN Y, FAN D, WANG Z, et al. System dynamic behavior modelling based on extended GO methodology[J]. IEEE Access, 2018, 6: 22513-22523. DOI:10.1109/ACCESS.2018.2816165 |
| [20] | LENVITIN G. A universal generating function approach for the analysis of multi-state systems with dependent elements[J]. Reliability Engineering & System Safety, 2004, 84(3): 285-292. |
| [21] | JOSEPHLU C J, MEEKER W Q. Using degradation measures to estimate a time-to-failure distribution[J]. Technometrics, 1993, 35(2): 161-174. DOI:10.1080/00401706.1993.10485038 |
| [22] | 朱新宇, 王有隆, 胡焱. 民航飞机电气仪表及通信系统[M]. 成都: 西南交通大学出版社, 2006. ZHU X Y, WANG Y L, HU Y. Electrical instrument and communications systems in civil aircraft[M]. Chengdu: Southwest Jiaotong University Press, 2006. (in Chinese) |
