时间：2023-11-04 13:09:38
[关键词] 风险 波动性 非参数GARCH模型GARCH模型
自从马科维茨的资产组合理论问世后,资产组合的收益率波动成为各类投资者和金融经济学家们长期关注的一个焦点问题。Robert Engle提出了著名的用来刻画资产波动的ARCH模型。Bollerslev在Engle的研究基础上进一步发展了GARCH(p,q)模型。实证结果验证, GARCH模型比较成熟,具有很多优点,也易于估计。然而, GARCH模型却依赖于模型的具体设定形式和误差分布形式,因此存在着模型误设的问题,一些GARCH模型并不能完美的估计波动性。而非参数GARCH模型在这方面可以弥补参数GARCH模型的这些劣势。
一、GARCH模型的表现形式
对数收益率序列,我们假定其均值方程是一个ARMA模型,设是均值修正的对数收益率,若:
(1)
其中是一个独立同分布的随机变量序列,均值为0,方差为1,参数,(这里对 ),则称服从GARCH(m,s)模型。对的限制条件保证的无条件方差是有限的,同时它的条件方差是随时间变化的。如前面一样,通常假设是标准正态分布或标准化的学生-t分布。若s=0,则上述方程就变为一个纯ARCH(m)模型。
二、非参数GRACH (p ,q) 模型思想及估计算法
设为一平稳时间序列, 表示由直到时刻 的过去信息产生的域, 的非参数GARCH(p ,q)模型可描述为:
(2)
(3)
其中为独立同分布的随机误差序列,均值为零,方差为,并且存在有限四阶矩。与相互独立,并且的严格正值函数:就是的条件方差,就是通常所说的波动率。
三、非参数GARCH模型的具体应用步骤
将3 改写成具有可加性误差形式的模型:
(4)
其中为一个鞅差序列,其均值;时,协方差,因而,以及,这表明,我们可以通过作 对滞后变量的非参数回归来获取函数 的非参数估计,然而波动率现在是一个无法观测的潜在变量。假设为来自过程GARCH模型的样本,Buhlmann 和McNeil (2002) 给出如下在比较弱条件下具有一致收敛性的估计算法:
(1)对一般的标准参数GARCH(p,q) 模型运用极大似然法计算波动率的第一次估计{},令m=1。
(2) 以为权数,用对和做作个变量的非参数加权回归,获得 的估计。
(3) 计算,边界值用代替。
(4) 增加m,如果,回到第(2)步,其中M为事先设定的最大迭代次数。
(5) 考虑最后K步波动率的平均,然后作对和的非参数回归,获得函数的最终估计以及条件以方差的最终估计形式。
四、非参数GARCH模型和GARCH模型对比
非参数GARCH模型具有广泛的应用范围,同时,非参数GARCH模型放宽了参数GARCH模型的很多限制,例如模型的形式和误差分布,为波动性问题的研究提供了非常有力的工具。使用非参数GARCH模型来拟合数据可以减少模型由于设定误差带来的偏差,在较大的范围内拟合模型,寻找更适合的线性形式。GARCH 模型簇已经成为度量金融市场波动性的强有力工具。许多专家和学者做了许多关于各种GARCH模型的理论及应用、波动率预测应用。可以说在研究股市波动方面,非参数GARCH模型已被证明是非常成功的。在其他方面,非参数GARCH模型的优势还有待发掘。尤其是我国基金市场,其波动性和股市的波动性有非常紧密的联系,具有许多相似的性质,因此我们利用非参数GARCH模型研究我国基金等投资市场的波动性具有很好前景。
参考文献:
[1]Engle, R. F.. Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation[J]. Econometrica, 1982, 50, 987-1008
[2]Bollerslev,T.. Generalized autoregression conditional heteroskedasticity[J]. Journal of Econometrics, 1986, 31, 307-327
[3]Xiangdong Long Aman Ullah.Nonparametric and Semiparametric Multivariate GARCH Model[J]
[4]鲁万波:基于非参数GARCH 模型的中国股市波动性预测[J].数理统计与管理,2006,4(25),455-461
[5]许冰,任军峰:基于非参数GARCH一种波动率估计方法[J].统计与决策,2006,10,6-10
关键词：瞬时测频；多项式拟合法；脉内参数估计
雷达信号的脉内特征是雷达信号相位特征的重要体现，是电子侦察中对雷达信号分选和识别的重要参数。因此，要可靠的分选和识别雷达信号就必须对雷达信号进行脉内特征分析。
在各种脉内分析方法中，分析信号的瞬时频率变化规律是进行脉内调制分析的有效手段。由于信号的频率是信号相位的变化率，所以对信号相位的估计自然可以用来对信号瞬时频率进行估计，从而进行相关脉内参数的估计，因此该方法的前提是能够得到准确的瞬时频率。
文中以线性调频信号以及相位编码信号为例，从侦察接收机的信号相位入手，利用瞬时相位法实现对瞬时频率的快速、准确的推算，然后利用多项式拟合法[1]来估计信号的相位参数。在Matlab中进行了仿真验证，同时给出了在某电子侦察系统中工程实现的过程。
4 结束语
本文讨论了利用信号相位建模进行瞬时频率估计算法，将频率的估计转化为求解信号模型参数的问题。在满足白噪声的情况下，由于在参数的估计过程中使用了最小二乘方法，因此能够有效的消除噪声的干扰，获得对信号瞬时频率的准确估计。基于瞬时频率进行脉内信号的调制特征的识别，还能够实现对常规脉冲、线性调频信号、非线性调频信号以及相位编码信号等调制类型信号的识别，提取相应信号的调制参数。该方法在一定信噪比条件下有较高的正确识别率，对脉宽变化不敏感，更加符合雷达侦察数字接收机上高速实现，并且在某雷达侦察数字接收机系统中已经取得成功验证，但是该类算法的问题在于，必须精确的建立相位模型，否则估计性能会受到很大的影响。
参考文献
[1] 林琳. 多项式相位信号的参数估计.实验科学与技术. 2006年6月第三期.
[2] 张春杰，郜丽鹏，司锡才. 瞬时相位法线性调频信号瞬时频率提取技术研究，弹箭与制导学报 TN911 2006.
[3] 齐国清，贾欣乐. 插值FFT估计正弦信号频率的精度分析. 电子学报 Vol.32 No.4 Apr. 2004.
[4] 胡来招. 数字瞬时测频—相位推算法测频[ J] . 电子对抗，2005.
课堂满意度 结构方程模型 参数估计
一、引言
课堂满意度是高校教学管理以及教育理论中重要的一个内容。目前很多高校都有学评教系统，课程结束时都要由学生对该课程作课堂满意度评价。在这些系统中，评价的指标数量比较少，各指标的权重分配有时不很合理，参评的学生人数也比较少，更多时候，由于种种原因学生会随意填写评价结果，这极大影响了评价结果的准确性和客观性。如果要从定量角度客观衡量学生课堂满意度就要建立定量模型。结构方程模型不失为一种有效的评价方法，可以将课堂满意度指数和无法直接观察的潜在变量以及可观察变量联系起来，根据样本数据情况对模型检验并修正后计算变量之间相关程度。
二、研究方法及变量选择
我们将设计主要针对在校大学生的调查问卷获得样本数据。根据实际情况，选择浙江省杭州市几所大学教学楼、图书馆、生活区等地为发放地点，随机向大学生发放问卷表，收回并筛选整理，据此分别为可观察变量赋值，从而建立课堂满意度结构方程模型。结构方程模型参数估计方法很多，在AMOS软件中，默认是极大似然估计法（Maximum likelihood，ML）。考虑到我们这次问卷设置的变量间关系不很明确，且变量不一定能完全满足正态分布，我们决定用偏最小二乘法（PLS）估计模型参数。
此次调查，共收回问卷356份，通过整理筛选，确定可用的338份。在问卷中，题目和答案的设置借鉴李克特量表（Likert scale）5级评价法，即最低1分，然后2分、3分，4分，最高5分，表示从非常不满意到非常满意[1][2]。
四、模型检验
利用AMOS软件，可以求得模型拟合优度评价指标值，其中GFI（非常规拟合指数）=0.903，CFI（比较拟合指数）=0.932，RMSEA（近似误差均方根）=0.023，GFI、CFI指标值越接近1，RMSEA小于0.1，表明模型拟合程度越好，从计算结果看，我们构建的模型是可以接受的。
对于各指标的效度，计算卡方值为441.52（sig=0.000），说明模型的变量具有充分代表性。
信度的检验，需要计算各潜变量的Cronbach's Alpha系数，而且该系数值要大于0.7才能通过信度检验。我们计算的5个潜变量α值都超过了0.8，说明都达到了信度标准。
五、模型参数估计
在AMOS中通过最大似然估计法，可以显示输出一个路径图，同时显示各项参数值即路径系数。根据前面研究设计，我们采用偏最小二乘回归法（PLS）估计模型的参数，各路径系数通过以下矩阵方程的形式显示[4]。
1.内生潜变量与观察变量的测量方程
对课堂满意度影响较大的因素，是Y2和Y3，学生对该课程的兴趣及接受程度决定了对课堂是否满意。比较而言，课程对生活实践是否有帮助及是否能从老师身上学到课外的东西显得不特别重要，这种情况应该和现在大学里很多课程注重期末考试且过分强调期末分数有关。
关键词：生长曲线，参数估计，伴随同化
0 引言
生长曲线(Logistic curve)也称S曲线，它是描述单一种群空间约束的生长过程曲线。其特点是开始生长较为缓慢，以后随着某些条件的变化，在某一段时间内增长速度较快，当达到某一界限之后，生长速度又趋于缓慢，以至最后停止增长，生长曲线的特征决定了其在生命科学领域中的广泛应用。论文大全。目前，生长曲线在其他领域中也得到广泛应用。例如向前忠将生长曲线模型用于高速公路诱增交通量预测，王吉权等将生长曲线用于电力负荷预测中。生长曲线的一般形式为
(1)
这里是某物种数量，、、是三个参数，应用时通常需要识别。参数识别是生长曲线
模型应用的前提，目前已有一些研究结果。如果令，，，则是如下方程的解
(2)
通常已知，于是只需要识别参数和，方程式（2）即为著名的Logistic模型。这里利用伴随同化方法对生长曲线的参数进行识别，同时将该方法用于文献[1]和美国1790-1950年人口数据。
1伴随同化参数识别方法
令为的观测，定义代价函数
(3)
这里为权重，为观测算子，为研究区间。代价函数是度量观测与模型解之间的距离函数，它反映在区间上与的拟合程度。于是模型参数识别问题就转换为以(2)为约束，以(3)目标函数的约束的极小值问题
(4)
构造拉格朗日函数
(5)
这里为的伴随变量。依据取极值的条件，容易得到满足
(6)
方程(6)称为方程(2)的伴随方程，需要逆向求解。依据(5)可计算代价函数关于模型参数的梯度
(7)
为方便，记
,(8)
于是可对模型参数进行校正
(9)
从而达到识别模型参数的目的。通常采用差分方法数值求解(2)和(6)，这里采用精度较高的4阶Rounge-Kutta方法，但要注意(6)要逆向求解。归纳起来利用伴随同化方法识别生长曲线参数的步骤如下：
1) 正向积分方程 (2)；
2) 逆向积分方程 (6)；
3) 计算梯度和代价函数；
4) 调整参数，为步长；
5) 如果则迭代终止(为事先给定的迭代终止参数)，否则转(1)。
2 数值实验
2.1 基于文献[1-2]数据的数值实验
由文献[1-2]可知，某种大豆的叶面指数y(t)与生育日数t的关系如表1的第一行和第2行。第3行为本文的结果，第4行为文献[2]的结果。通过表1可看出，本文方法可以较好地识别出参数值，本文得到生长曲线
(10)
表1 数值实验结果
L 1+ae-bx
中参数L，a，b的估计往往是基于一定的经济意义或生物统计的背景，但从纯粹的数学方程角度来看，当L，a，b的取值范围不加任何限制时，是否也能估计出L，a，b呢？文章证明了当（L，a，b）
∈R3
时，y=
L 1+ae-bx
中参数L，a，b不能用最小二乘法估计，并给出了其参数在一定的限制条件下，可以用最小二乘法估计.
[关键词] 逻辑曲线 参数 最小二乘估计
[中图分类号] G633.6 [文献标识码] A [文章编号] 1674 6058（2016）17 0049
一、引言
基于残差平方和最小的思想，最小二乘法已在社会生产和实践中得到了广泛的应用.在解决非线性回归模型的参数估计问题时，我们往往将其转化为线性回归模型，进而用最小二乘法估计出其参数，但有的模型并不能通过直观的观察和简单的计算就能转化为线性回归模型，也就谈不上用最小二乘法估计其参数了.这就引发了我们对非线性回归模型的参数是否能用最小二乘法估计的沉思，如常见的逻辑曲线y=
L 1+ae-bx
（（L，a，b∈R3）），其参数L，a，b能否用最小二乘法的思路估计出或至少能被确定属于一定的范围呢？这就是本文所要研究的重点.
二、猜想
任意取n个不同的样本（xi，yi），i=1，2，…，n.只要样本中不含用（0，y（0））；那么逻辑曲线y= L 1+ae-bx （（L，a，b）∈R3）的参数L，a，b就不能用最小二乘法估计出且不能估出（L，a，b）属于R3中某一真子空间，即无法确定（L，a，b）∈H，HR3.
三、证明
假设用最小二乘法可以估计出（L，a，b），尽管（L，a，b）并不一定唯一，即（L，a，b）也可属于R3中一真子空间.
第一步，观察y= L 1+ae-bx ，
（xi，yi）可取n个不同的样本，
可推断出L的估计，L≠0.
第二步，（Ⅰ）y= L 1+ae-bx （*）两边同时对x求导，得
dy dx =by（1-
y L
）
，其中，y（0）=
L 1+a
（**），
参数（L，a，b）的估计若满足（*）式，则也满足（**）.
（Ⅱ）考查微分方程y= L 1+ae-bx ，y（0）= L 1+a （**）.
令f（x，y）=by（
1- y L ）
，取R3平面上一区域D：|x| ≤h1|y|≤h2.
不难发现：f（x，y）在D上连续，且
|f（x1，y1）-f（x2，y2）|= 1 L |bL
（y1-y2）
-b（y21-y22）|≤
b|y1-y2|+2
bh2 L
|y1-y2|=（b+2 bh2 L ）|y1-y2|.
故f（x，y）在D上满足Lipschitz条件.
由常系数微分方程的解存在唯一性定理知，（**） 式满足y（0）= L 1+a 的解在|x|≤h上存在且唯一.
这里的h=min{h1，h2}，M=max|f（x，y）|，（x，y∈D）.
（Ⅲ）现在解（**）式：
1）观察f（x，y），得：y=0或y=L为（**）的解.
2）在f（x，y）中，令z= 1 y ，则
dz dx =
dz dy × dy dx =
- 1 y2 ・ b L y（L-y）
=-bz+
b L .
（1）
再令z=c（x）e-bx，则
dz dx =
dc（x） dx e-bx-
be-bxc（x）
= dc（x） dx e-bx
-bz. （2）
比较（1）（2）得： dc（x） dx = b L ebx.
利用变量分离法，求得：c（x）= 1 L ebx
+C，C为任一常数.
y= 1 z =
1 c（x）
ebx=
L ebx+CL ebx=
L 1+CLe-bx .
（3）
将y（0）=
L 1+a
代入（3），得：C= a L .
y= L 1+ae-bx .
综上，（**）式的解为y=0或L或 L 1+ae-bx .
但由（Ⅰ）知，样本（xi，yi）对（**）式中参数L，a，b的估计问题同样适用.故实际操作中，由于n个样本（xi，yi）是不同的，所以y=0，L要排除.
综合（Ⅱ）（Ⅲ）知，在整个实数域 R 上，（**）式的解都是存在参数（L，a，b）作最小二乘估计.
第三步，多元线性回归方程的一般形式是y=β0+β1x1+…+βpxp，不失一般性，总可以设β0=0.因为可以形式上引入自变量x0=1，所以在理论研究时，可以假设多元线性回归模型的形式为y=β1x1+…+βpxp.
由假设知，y= L 1+ae-bx 中参数a，b可作最小二乘估计，故y= L 1+ae-bx 必可划为线性回归形式：g（x，y）=p（L）u（x，y）+q（a）u（x，y）+h（b）w（x，y）.（4）
由于我们是先估计出p（L），q（a），h（b）后，利用L=p（L），
q（a）=q（a），h（b）=h（b），进而估计出，L，a，b的，所以（4）式与y= L 1+ae-bx 实质上对L，a，b的估计并没有任何区别.
故不妨设y= L 1+ae-bx 可划为y=Lx1+ax2+bx3（5）.
（5）式中，残差θi=
（yi-Lxi1-axi2-bxi3）=（L，a，b）xx′（Lab）-2y′x（Lab）+y′y.
其中，x=x11x22…xnn，y=y1y2…yn，
2xx′（Lab）-2x′y=0.
由假设知，L，a，b可用最小二乘估计出：若（L，a，b）唯一，则r（xx′）=3；若（L，a，b）不唯一，则（L，a，b）属于R3中一真子空间，此时r（xx′）=1或2.
综上，由假设得到r（xx′）>0.
第四步，由第二步知，对y= L 1+ae-bx 中参数L，a，b的估计等价于对 dy dx
=by（1- y L ），y（0）= L 1+a 中参数L，a，b的估计.
分析后者，我们发现通过样本（xi，yi）至多只能直接估出b，L，而a与样本（xi，yi）并没有直接关系，它只与y（0）有直接联系，a= L y（0） -1（由第一步知，L≠0，y（0）= L 1+a ≠0）.
但实际的样本数据中，由于（0，y（0））并未给出，故y（0）的取值范围为R＼{0}.
现任意取定L，b，一方面残差θ退化为只含a的二次函数，即r，s，t∈R1θ=θ（a）=ra2+sa+t，
则a的估计值必须满足θ′（a）=0.即2ra+s=0. （6）
另一方面，a=
L y（0）
-1，由L≠0，y（0）的取值范围为R＼{0}，得到a的取值范围为R＼{-1}，故R＼{-1}中任一点都可作为a，显然也要满足（6）式.
r=s=0，
θ（a）=t为常数，故残差θ与a无关.
同理，任意取定a，b，也可得到残差θ与L无关.
故残差θ至多与b有关.
如果θ与b有关，则θ=θ（b），即y= L 1+ae-bx 可转化为线性回归形式，g（x，y）=h（b）w（x，y）；但前式含有3个参数L，a，b，后式只含有1个参数b，显然，这种转化不成立.
故残差θ=t为常数.
反馈到第三步，可知此时的xx′=03×3，即r（xx′）=0，与假设得到的r（xx′）>0矛盾.
原假设不成立，故猜想成立.
证毕.
四、总结与归纳
尽管L，a，b在非限制条件下，y= L 1+ae-bx 中的参数L，a，b不能用最小二乘法估计，但在实际经济模型或生物统计模型中，L，a，b往往具有特殊的经济含义和现实意义.我们可以用其他方法先估计出L或b等，然后再利用最小二乘法估计剩余的两个参数.比如，著名的Logistic人口发展模型y= L 1+ae-bx 中的参数L代表自然资源和环境条件所能容纳的最大人口数量，我们可以先根据生态学的知识预测出L，然后将y= L 1+ae-bx 转化为线性回归形式，即
-bx+lna=ln（ L y -1）
，再利用不同年份的人口统计数据和最小二乘法估计出a，b.
[ 参 考 文 献 ]
[1] 王高雄.周之铭等.常微分方程[M].北京：高等教育出版社，2007.
[2] 北大数学系几何与代数教研室前代数小组.高等代数[M].北京：高等教育出版社，2003.
[3] 陈家鼎，孙山泽等.数量统计学讲义[M].北京：高等教育出版社，2011.
【摘要】
提出了药物动力学(PK)参数估计的非线性加权最小一乘法(ML1)。目标函数意义为使加权绝对残差和或相对平均绝对误差(AARE)趋于最小。采用的遗传算法(GA)对目标函数性质几无要求，易得全局优化参数估计值，且稳健性好。应用GA作PK模型参数估计，效果较佳。
【关键词】 药物动力学 参数估计 稳健性 权重 遗传算法
在目前的文献［1～4］中，药物动力学(PK)模型拟合时的目标函数Q一般为：minQ=ΣWiε2i (i=1，2，3，…，n)(1)
式(1)中，Q为残差或加权残差平方和；Wi为权重(一般为1，1／Ci或1／C2i等)；εi为残差即药物浓度测定值与相应PK模型(参数为θj，j=1，2，3，…，p，p为独立参数个数)计算值或预测值之差，εi独立同分布，εi的期望即均值为零，各点的εi方差相等或近似相等均为σ2(此时不需加权或可认为需加权，只是Wi均为1)或不相等分别为σ2i(此时需加权处理)；n为样本含量即药(C)时(t)数据对个数。根据非线性最小二乘法(LS)及最优化原理，使Q趋于最小值，则得到欲求PK参数θ的LS估计值。
一般地，Ct曲线高低浓度相差较大，可相差数倍、数十倍乃至数百倍、在绝大多数情况下，各εi方差σ2i不相等或差异较大，故PK模型拟合时必需作加权处理，否则，低浓度数据的作用将会被忽视，常会出现低浓度处的拟合结果失真，拟合曲线与观察点有较大差距。作加权处理后，不但可在一定程度上解决上述问题，而且增强了LS估计的稳健性［5］。这是PK模型参数辨识的一个显著特点。在估算PK参数时，如何选择Wi，在实践和理论上都是一个重要问题，但目前尚无统一的办法加以解决。只能按数据的实际情况，以及在不同权重下拟合后，考察并比较各自的一些统计学指标，择优选取。
作者认为，以式(1)为目标函数，由于其固有的缺陷，无法统一Wi。因此，我们所用的目标函数(见式2)，意义明确，权重统一，而且用基于Matlab的遗传算法(GA)［6，7］作参数估计，对目标函数的性质几无要求，稳健性更强，易得全局优化的参数估计值，从而解决了式(2)的计算问题。
1 原理与方法
定义目标函数J为：
J(θj，t)=Σ(Wi×|εi|)=Σ(1／Ci×|εi|)min (i=1，2，3，…，n，n为观察点数)(2) 式(2)中，Wi及εi含义：
① Wi为权重，因不能得到“最优”权重，故只能寻求“较优”权重。原则是：对大的|εi|，应予小的权重；对小的|εi|，应予大的权重。若Wi=1，则式(2)为最小绝对残差和即最小一乘(L1)，该法未对|εi|加权，实际上把各|εi|的重要性同等看待，不尽合理；若Wi=|εi|，则式(2)为最小残差平方和即普通最小二乘(LS)，实际上是对大的|εi|，予大的权重，对小的|εi|，予小的权重，这显然更不合理。结合药物动力学及体内药物测定的特点，一般地，药物浓度高时，|εi|大；药物浓度低时，|εi|小。假定数据的误差|εi|与数据的大小Ci基本成比例变化，故取Wi=1／Ci较为合理。
② εi为观测值与似合值之差即残差：
εi=Ci-F(θj，t)
(j=1，2，3，…，p，p为独立参数个数)(3)
式(3)中，Ci-t为药物实验观测值时间数据，F为一非线性函数即具有p个独立参数θj的已知PK模型的似合值或预测值，n为样本含量即药物(C)时间(t)实验观察点数。因Σεi0，故每个εi取其绝对值。
从权重的角度理解式(2)，|εi|为误差(残差)εi的绝对值，对于大的误差，应给予小的权重；对于小的误差，应给予大的权重。这样|εi|的权重为1／Ci是恰当的。故式(2)可视为使加权绝对残差和趋于最小(加权最小一乘WL1)。F(θj，t)为一非线性函数，因而式(2)为非线性WL1。
从误差或预测精度角度理解式(2)，J为拟合值与观测值的相对误 差之和。J越小，则总相对误差越小。由于对ct数据而言，实验点数n已知，故平均相对误差亦越小，表明拟合效果较佳。因此，式(2)亦可理解为使用相对平均绝对误差(AARE)［8，9］趋于最小。
根据最优化原理，求得使J(θj，t)趋于最小值时的θj的估计值，即为欲求PK参数估计值。但对式(2)编程计算是困难的，应用基于Matlab的遗传算法(GA)可较易解决所面临的问题。
转贴于
GA在解决优化问题时，具有许多优点如全局搜索与寻优、稳健性好、不受优化函数连续可微的约束、不依赖于梯度信息、不要求优化函数导数必须存在、对目标函数的性质几无要求、搜索具有探索性或启发性等。因此，这里采用GA搜索目标函数式(2)的参数θj的最优化估计值。方法如下：
在P4计算机及Matlab 7.0平台上，首先根据式(2)或(3)，编制适应度函数即目标函数，然后利用遗传算法与直接搜索工具箱有关函数，设定适当运算参数值及PK参数的大致初始范围。一般数分钟内即可得到目标函数近似最小值及θj估计值。
2 结果
2.1 实例
大鼠按30mg/kg静脉注射盐酸川芎嗪，体内PK过程呈开放2房室模型［10］。不同时间的血浆药物浓度平均值及拟合结果见表1。表1 血药浓度(C)时间(t)实验数据及用遗传算法按不同目标函数(OBJ)拟合后结果（略）
由表1可见，以式(2)为目标函数，拟合后相对平均绝对误差最小，表明本法预测精度高。
2.2 数值模拟试验
用常见静脉注射2房室PK模型C=A*exp(-α*t)+B*exp(-β*t)，设定参数真值后拟合结果见表2。表2 模拟Ct数据及用遗传算法按不同目标函数(略)
由表2可见，在只有一个异常点C(1.5)为+20%误差时，本法拟合精度最高，参数估计值与设定真值最为接近，表明本法有较强的稳健性和可靠性。其他数值模拟试验亦证明了这一点。
3 讨论
由上述分析及结果可知，本研究提出的PK参数估计的非线性WL1或最小AARE法，具有3个特点：首先，解决了PK参数估计时的权重选择问题，即把权重Wi统一为1／Ci；其次，本法在估计PK参数时LS或加权LS更具稳健性；最后，目标函数或评价函数尚具有更确切涵义即使残差εi的相对平均绝对误差或加权绝对残差和趋于最小。对式(2)作参数估计时，由于涉及绝对值运算，而GA对目标函数性质无要求且全局优化功能强，因此，本研究均用基于Matlab的GA作计算，效果极佳。而且单纯型算法作式(2)最优化计算，对参数初值要求高，常出现过早收敛，只能得到局部优化参数估计值。
【参考文献】
1 宋振玉.药物代谢研究—意义，方法，应用.第1版.北京：人民卫生出版社，1990，159～165.
2 凌树森.治疗药物监测新理论与新方法.第1版.北京：中国医药科技出版社，2002，69～72.
3 曾衍霖.药物代谢动力学中的二个问题—原始数据的权重与线性数学模型中房室数的确定.药学学报，1980，15(9)：571.
4 刘昌孝，孙瑞元.药物评价实验设计与统计学基础.第1版.北京：军事医学科学出版社，1999，119～130.
5 《现代数学手册》编纂委员会.现代数学手册(随机数学卷).第1版.武汉：华中科技大学出版社，2000，619～619.
6 雷英杰，等.MATLAB遗传算法工具箱及应用.第1版.西安：西安电子科技大学出版社，2005，146～207.
7 薛定宇，陈阳泉.高等应用数学问题的MATLAB求解.第1版.北京：清华大学出版社，2004，357～369.
8 李华，胡奇英.预测与决策.第1版.西安：西安电子科技大学出版社，2005，150～153.
9 王晓燕，等.医学预测方法的精度分析.数理医药学杂志，2001，14(5)：390.
关键词：参数估计 适线法 模糊加权优化 相对隶属度
1 问题的提出
P—Ⅲ型分布做为我国水文分析计算中规定的线型，长期以来得到广泛应用，我国水文工作者对其参数估计问题作了大量研究工作，先后提出了适线法、权函数法、数值积分权函数法、Monte—Carlo法、概率权重矩法、下限法、动点动线适线法、密度函数法等等.金光炎对适线法作了大量工作［1］.适线法实际是一种优选法，即优选一组参数使之相应的理论曲线与经验点据最好吻合.但由于目估适线没有一个明确定量的拟合优度标准，使结果因人而异，往往达不到优选目的，找不到最佳结果.为了有一个客观定量标准，并利用计算机适线，近年来产生了一种新的方法——模型搜索适线法，其原理同目估适线法一样，寻找一组参数，使这组参数所决定的理论曲线与经验点据拟合最优.由于不同适线准则与不同的经验频率计算公式得到的结果有时差异很大，为了确定合理的适线准则和相应的经验频率计算公式，丛树铮等［2］对此问题作了大量的水文统计试验研究.
其实，在采用优化适线法时，尚有一个经验点与所配曲线之间离差的加权问题，由于各经验点的精度不同、在频率曲线上的位置不同、配线的目的不同，其在优化适线时的重要程度理应不同.现有的优化适线法均不能考虑这些因素，只能将各经验点与所配频率曲线之间的离差视为等权，因此计算结果受个别异常点影响十分敏感，致使所配曲线受少数异常点控制而时常偏离大多数经验点的分布趋势，这正是现行优化适线法计算结果欠稳定、平差能力低的症结所在.如何合理的确定各经验点与所配曲线之间离差在优化适线中的权重，使参数估计具有更好的稳健性与良好的平差能力，是一个值得探讨的问题.本文根据模糊数学的基本原理［3］，首次提出了经验点对所配曲线的相对隶属度概念及该相对隶属度的分析计算方法，进而提出了以该隶属度为权重的P—Ⅲ型分布参数模糊加权优化适线法，以期克服现行优化适线法之不足.
2 模糊加权优化适线法
若随机变量X服从P—Ⅲ型分布，则其概率度函数和分布函数分别为：
(1)(2)
其中α，β，γ为3个未知参数，它们与常用的3个特征参数Ex,Cv,Cs有如下关系
(3)
设想自总体X中抽取K组容量为M的样本系列，则可用适线法选配出K条频率曲线，若有某条频率曲线恰好完全通过Exm(即X(pm)=Exm,m=1,2,…，M)，则称此条频率曲线为随机变量X的理想最优频率曲线.因为只有通过它外延得出的设计值才是无偏的，其它频率曲线都是有偏的.
在理想最优频率曲线概念的基础上，我们进一步给出经验点据(xm,pm)对理想最优频率曲线的相对隶属度定义及其计算方法.
若以随机独立方式从随机变量X的总体中抽取一组容量为M的样本，将其按大小为序排列，则序位为m的变量xm(相应的经验频率为pm)称为次序统计量.它也是一个随机变量，若xm恰好落在理想最优频率曲线上，则xm对该曲线的相对隶属度为1；若xm距离理想最优频率曲线愈远则相对隶属度愈小，若xm距离曲线无穷远，则相对隶属度为0.
次序统计量xm是原随机变量X衍生出来的，其概率分布特性可以通过条件分布密度函数f(xm｜pm)来描述，按陈守煜教授提出的概率分析与模糊集分析相结合的模糊水文学基本思想［4］，本文假定条件分布密度函数f(xm｜pm)近似服从均值为Exm的正态分布，即f(xm｜pm)=N(Exm,σm)，并按下式定义xm对理想最优频率曲线的相对隶属度.
定义：若随机变量X的理想最优频率曲线为f(x,Ex,Cv,Cs)，则某经验点据(xm,pm)对f(x,Ex,Cv,Cs)的相对隶属度为
(4)
f(x,Ex,Cv,Cs)为理想最优频率曲线密度函数，其它符号同前.
在f(xm｜pm)服从正态分布的条件下，可按文献［2，5］所述的水文统计试验方法推求不同适线准则下式(4)中的σm，金光炎对此问题作了深入细致的研究，并绘制了绝对值准则和最小二乘法适线的B值诺模图［1，5］.
本文采用了金光炎的以上成果，并利用下式计算σm
式中S为X系列的标准差，n为样本容量，B为B值模图中的取值，它是Cs和pm的函数.
因为样本均值对期望值Ex的估计具有无偏性且方差最小，因而不需通过适线法来确定，于是便可得到本文提出的模糊加权优化适线法的数学模型。
(5)
式中c为优化准则参数.
3 各种参数估计方法的对比检验
3.1 计算结果 任何参数估计方法，即使在理论上正确，也必须通过大量适当的数据充分检验合格以后，才有应用的价值，本文采用刘光文教授［6］提出的与水文安全因数有关的绝对精度标准，即x1%和x0.1%(作为代表)的真误差不应超过8%—10%，用本文提出的模糊加权优化适线法和现行的几种参数估计方法，对大量的理想样本资料同时进行了对比计算，以下仅将与文献［6］中20组理想样本对
比的计算结果列于表1.
表1 假想样本对各种P—Ⅲ型参数估计方法的对比检验 样本序号* 样本容量(n) 项 目 假想理论值
(TV) 数值积分*
双权函数法
(NIDWF) 模糊加权最小
平方准则
(FWLSC) 模糊加权绝对值
最小准则
(FWLAVC) 1 19 cv
cs
x1%
x0.1%
1
1/2
1
2.51128
3.26556 0.99344(-0.7%)
0.48154(-3.7%)
0.96948(-3.1%)
2.42987(-3.2%)
3.14018(-3.8% 0.97942(2.1%)
0.5046(+1.3%)
0.8573(14.3%)
2.4184(-1.2%)
3.18867(-2.4% 0.97942(2.1%)
0.5089(+1.8%)
0.7474(25.3%)
2.44141(2.4%)
3.11913(4.5%) 2 29 cv
cs
x1%
x0.1% 1
1/2
1
2.51128
3.26556 0.99561(-0.4%)
0.48796(-2.4%)
0.98480(-1.5%)
2.45920(-2.1%)
3.18634(-2.4%) 0.98454(1.5%)
0.5097(+1.9%)
0.9041(-9.6%)
2.51049(0.0%)
3.24205(1.0%) 0.98454(1.5%)
0.4953(-0.9%)
0.6418(35.8%)
2.38095(5.2%)
2.99158(8.4%) 3 19 cv
cs
x1%
x0.1% 1
1/2
3/2
2.66518
3.61676 0.99072(-0.9%)
0.47775(-4.5%)
1.47686(-1.5%)
2.56057(-3.9%)
3.45267(-4.5%) 0.97031(3.0%)
0.5119(+2.4%)
1.2878(14.1%)
2.63823(1.0%)
3.52551(2.5%) 0.97031(-3.0%)
0.5084(+1.7%)
1.1152(-25.7%)
2.57268(-3.5%)
3.38809(-6.3%) 4 29
cv
cs
x1%
x0.1% 1
1/2
3/2
2.66518
3.61676 0.909381(-0.6%)
0.48511(-3.0%)
1.48752(-0.8%)
2.59586(-2.6%)
3.50864(-3.0%) 0.97775(-2.2%)
0.5095(+1.9%)
1.3440(-10.4%)
2.65240(-0.5%)
3.55990(-1.6%) 0.97775(-2.2%)
0.5115(+2.3%)
1.2061(-19.6%)
2.61539(-1.9%)
3.47005(-4.1%) 5 19
cv
cs
x1%
x0.1% 1
3/4
3/2
3.49776
4.92515 0.98608(-1.4%)
0.72573(-3.2%)
1.47852(-1.4%)
3.56030(-3.9%)
4.71006(-4.4%) 0.95547(-4.5%)
0.7839(+4.5%)
1.1547(-23.0%)
3.44339(-1.6%)
4.74307(-3.7%) 0.95547(-4.5%)
0.7689(+2.5%)
1.0853(-29.4%)
3.3642(-3.8%)
4.57124(-7.2%) 6 29
cv
cs
x1%
x0.1% 1
3/4
3/2
3.49776
4.92515 0.99072(-0.9%)
0.73903(-1.5%)
1.48737(-0.8%)
3.42364(-2.1%)
4.80974(-2.3%) 0.96663(-3.3%)
0.7724(+3.0%)
1.2649(-15.7%)
3.46916(-0.8%)
4.70729(-2.6%) 0.96663(-3.3%)
0.7580(+1.1%)
1.2038(-19.7%)
3.39449(-3.0%)
4.66096(-5.4%) 7 19
cv
cs
x1%
x0.1% 1
1/2
2
2.83259
3.95388 0.98850(-1.2%)
0.47137(-5.7%)
2.00927(+0.5%)
2.67059(-4.7%)
3.74683(-5.2%) 0.96261(-3.7%)
0.5066(+1.3%)
1.8071(-9.6%)
2.77295(-1.6%)
3.86085(-2.4%) 0.96262(-3.7%)
0.5058(+1.2%)
1.5219(-23.9%)
2.68934(-4.0%)
3.66066(-7.4%) 8 29
cv
cs
x1%
x0.1% 1
1/2
2
2.80259
3.95388 0.99233(-0.8%)
0.48052(-3.9%)
2.00328(+0.2%)
2.71220(-3.2%)
3.81141(-3.6%) 0.97203(-2.8%)
0.5049(+1.0%)
1.8811(-5.9%)
2.79284(-0.3%)
3.90880(-1.1%) 0.97203(-2.8%)
0.5073(+1.5%)
1.6324(-18.4%)
2.73194(-2.5%)
3.75257(-5.1%) 9 19
cv
cs
x1%
x0.1% 1
2/3
2
3.40345
4.93850 0.98466(-1.5%)
0.64245(-3.6%)
2.00993(+0.5%)
3.26847(-4.0%)
4.72963(-4.2%) 0.94014(-5.0%)
0.6951(+4.2%)
1.5256(-2.37%)
3.32663(-2.3%)
4.83837(-2.0%) 0.95014(-5.0%)
0.6954(+4.3%)
1.5155(-24.2%)
3.31947(-2.5%)
4.65110(-5.8%) 续表1
样本序号* 样本容量
(n) 项 目 假想理论值
(TV) 数值积分*
双权函数法
(NIDWF) 模糊加权最小
平方准则
(FWLSC) 模糊加权绝对值
最小准则
(FWLAVC) 10 29 cv
cs
x1%
x0.1% 1
2/3
2
3.40345
4.93850 0.98978(-1.0%)
0.65060(-2.4%)
2.00250(+0.1%)
3.31215(-2.7%)
4.79620(-2.9%) 0.96271(-3.7%)
0.6940(+4.0%)
1.7077(-14.6%)
3.39963(-0.1%)
4.83837(-2.0%) 0.96271(-3.7%)
0.6930(+3.9%)
1.6239(-18.8%)
3.36387(-1.2%)
4.75381(-3.7%) 11 19 cv
cs
x1%
x0.1% 1
1
2
4.60517
6.90776 0.97699(-2.3%)
0.99364(-0.6%)
2.01400(+0.7%)
4.48373(-2.6%)
6.72992(-2.6%) 0.92521(-7.5%)
1.0652(+6.5%)
1.9021(-4.9%)
4.77998(+3.8%)
7.14573(+3.4%) 0.92521(-7.5%)
0.8599(-14.0%)
1.4378(-28.1%)
3.82762(-16.9%)
5.41929(-21.5%) 12 29
cv
cs
x1%
x0.1% 1
1
2
4.60517
6.90076 0.98466(-1.5%)
0.99481(-0.5%)
2.00226(+0.1%)
4.51724(-1.9%)
6.77452(-1.9%) 0.94406(-5.6%)
1.0472(+4.7%)
1.9021(+4.9%)
4.75414(+3.2%)
7.11889(+3.1%) 0.94406(-5.6%)
0.9684(-3.2%)
1.7022(-14.9%)
4.34908(-5.6%)
6.35447(-8.0%) 13 19
cv
cs
x1%
x0.1% 1
5/2
5/8
3.40337
5.09259 0.98316(-1.7%)
0.59299(-5.1%)
2.60188(+4.1%)
3.25111(-4.5%)
4.87426(-4.8%) 0.94549(-5.4%)
0.6308(+0.9%)
2.2220(-11.1%)
3.34111(-0.2%)
4.90459(-3.7%) 0.94549(-5.4%)
0.5990(-4.2%)
1.9462(-22.2%)
3.14045(-7.7%)
4.49276(-11.8%) 14 29
cv
cs
x1%
x0.1% 1
5/1/2
5/8
3.40337
5.09259 0.98873(-1.1%)
0.60223(-3.6%)
2.52765(+1.1%)
3.28902(-3.4%)
4.90818(-3.6%) 0.95920(-4.1%)
0.6354(-1.7%)
2.3285(-6.9%)
3.40037(0.0%)
5.03408(-1.1%) 0.95920(-4.1%)
0.6209(-0.7%)
2.0570(-17.7%)
3.26343(-4.1%)
4.72475(-7.2%) 15 19
cv
cs
x1%
x0.1% 1
3/4
3
4.03853
6.36426 0.97690(-2.3%)
0.72046(-3.9%)
3.14988(+5.0%)
3.86713(-4.2%)
6.13317(-3.6%) 0.92712(-7.3%)
0.7508(+0.1%)
2.5727(-14.2%)
3.90632(-3.3%)
5.97705(-6.1%) 0.92712(-7.3%)
0.6686(-10.9%)
2.4218(-19.3%)
3.54320(-12.3%)
5.30654(-16.6%) 16 29
cv
cs
x1%
x0.1% 1
3/4
3
4.03853
6.36426 0.98448(-1.6%)
0.73826(-2.9%)
3.07152(+2.4%)
3.90830(-3.3%)
6.17219(-2.0%) 0.94530(-5.5%)
0.7713(+2.8%)
2.7660(-7.8%)
4.06518(+0.7%)
6.31981(-0.7%) 0.94530(-5.5%)
0.7119(-5.1%)
2.4904(-17.0%)
3.74470(-7.3%)
5.66888(-10.9%) 17 19
cv
cs
x1%
x0.1% 1
1
3
5.05138
8.15235 0.96920(-3.1%)
0.99675(-0.3%)
3.16031(+5.3%)
4.97404(-1.5%)
8.16870(+0.2%) 0.90283(-9.7%)
1.1035(+10.4%)
2.2418(-25.3%)
5.09906(+0.9%)
7.84849(-3.7%) 0.90283(-9.7%)
1.0692(+6.9%)
2.3048(-23.2%)
5.00685(-0.9%)
7.72542(-5.2%) 18 29
cv
cs
x1%
x0.1% 1
1
3
5.05138
8.15235 0.97931(-2.1%)
0.99508(-0.5%)
3.08385(+2.8%)
4.95773(-1.9%)
8.04434(-1.3%) 0.92707(-7.3%)
1.0778(+7.8%)
2.4421(-18.6%)
5.13616(+1.6%)
8.01065(-1.7%) 0.92707(-7.3%)
1.0472(+4.7%)
2.4619(-17.9%)
5.02789(-0.5%)
7.83777(-3.9%) 19 19
cv
cs
x1%
x0.1% 1
3/2
3
7.07706
11.72853 0.95380(-4.6%)
1.55446(+3.6%)
3.18574(+6.2%)
7.06141(-0.2%)
11.87708(+1.3%) 0.85424(-14.6%)
1.6422(+9.5%)
2.7897(-7.0%)
7.49812(+5.9%)
12.30353(+4.9%) 0.85424(-14.6%)
1.3063(-12.9%)
2.43470(-12.9%)
5.96459(-15.7%)
9.41022(-19.8%) 20 29
cv
cs
x1%
x0.1% 1
3/2
3
7.07706
11.72853 0.96896(-3.1%)
1.52306(+1.5%)
3.08189(+2.7%)
6.99291(-1.2%)
11.66507(-0.5%) 0.89060(-10.9%)
1.5864(+5.7%)
2.8152(-6.2%)
7.35641(+3.9%)
12.07083(+2.9%) 0.89060(-10.9%)
1.3903(-7.3%)
2.6342(-12.2%)
6.46617(-8.6%)
10.3966(-11.4%) * 计算结果引自文献［6］
又以淮河蚌埠站平均流量系列为例，进行了计算，其结果见表2，图1.
表2 各种方法计算结果对照
估计方法 Ex Cv Cs 矩 法 827.5 0.65 1.20 极大似然法 827.5 0.72 1.88 适 线 法 827.5 0.75 1.88 概率权矩法* 827.5 0.68 1.69 模糊加权适线法 827.5 0.71 1.52 *计算结果引自文献［8］
3.2 分析与讨论 本文仅对模糊加权优化适线法进行讨论，关于其他方法的讨论请参考文献［6］.
表1的计算结果表明：模糊加权最小二乘准则，对20组理想样本资料全都能达到精度标准.x1%的真误差绝对值平均为1.6%；而且正负误差基本持平，其中误差小于5%的有19组，最大误差为-6.1%；x0.1%的真误差绝对值平均为2.4%，其中正误差1.0%.负误差1.4%，误差小于5%的有19组，最大一组误差+5.9%，没有随Cv,Cs变化而出现明显的系统偏大或偏小现象，对不同参数的资料适应性强，计算结果稳定.
图1 淮河蚌埠站年平均流量频率曲线
模糊加权绝对值和最小准则，设计值估值全是负误差.x1%的真误差平均为5.2%，误差小于5%的有12组，5%—10%的5组，大于10%的3组，最大误差16.9%；x0.1%的误差均值为8.7%，小于5%的有4组，5%—10%的10组，大于10%的6组，最大误差为21.5%，不能达到绝对精度标准，计算结果系统偏小，对工程设计偏于不安全.计算结果还表明：模糊加权最小二乘准则适线法所估计的参数Cv有系统偏大，Cs有系统偏小的现象，这主要是由于用矩法估计有系统偏小的原因所致.但最终的设计估计值x1%,x0.1%没有出现明显的系统偏差，计算精度很高，这正是由于本方法有在，Cv,Cs估计值偏小的情况下，使Cv适当偏大，在一定程度上补偿了,Cs估值偏小的影响，使最终设计估值有很高的稳定性和计算精度，及较强的资料适应能力的优点.而模糊加权绝对值和最小准则适线法在估计偏小的情况下，不能很好地协调Cv,Cs估值，补偿各参数的估计值偏差，使最终设计估值系统偏小明显，难以达到绝对精度要求.
淮河蚌埠站实测资料的计算也表明，模糊加权最小二乘准则适线法所获得的理论频率曲线与实测资料拟合良好.
4 结 语
P—Ⅲ型曲线作为我国水文计算规范所规定的频率分析线型，在我国水文频率计算中得到广泛应用.但对于P—Ⅲ型分布参数的估计问题长期以来没有得到很好的解决，不同学派众说纷云，难达一致.本文提出了理想最优频率曲线和经验点据对理想最优频率曲线的相对隶属度概念及该相对隶属度的分析计算方法.然后从充分利用样本信息，增强参数估计方法的稳健性和平差能力等方面考虑，提出了以经验点对理想最优频率曲线相对隶属度为权重的模糊加权优化适线法.
模糊加权优化适线法可以充分利用样本信息(包括经验点的大小，排位和可能存在的误差等信息)，区分不同样本点在配线分析中的作用，较好的反映了“参数估计应以精度精高的实测点据为主，不要过于崇信可能存在较大误差的历史洪水，但又要适当考虑历史洪水点据所带来的信息”等配线分析过程中的模糊性思想，并把这些思想模型化、规范化，克服了一般优化适线方法在配线过程中等同看待精度不同的各经验点据，最终导致计算结果受个别异常点影响甚大，计算结果不稳定的缺点.
利用与水文安全因素有关的绝对精度检验标准，以理想样本和实测样本对不同模糊加权优化准则(模糊加平方和最小、模糊加权绝对值和最小)进行计算的结果表明：模糊加权绝对值和最小优化适线法，Cv,Cs,x1%,x0.1%均系统偏小，最大误差为-16.9%—-21.5%，难以达到容许误差标准；模糊加权平方和最小优化适线法，虽,Cs系统偏小，Cv系统偏大，但x1%,x0.1%计算精度很高，计算结果稳定，分析参数有系统偏差的原因，主要因用矩法估计系统偏小所致，Cv偏大是为补偿因，Cs估值偏小，使最终设计估计值稳定.
模糊加权平方和最小优化适线法计算结果稳定，对不同参数资料适应性强且参数,Cv,Cs间有很好的补偿效果，设计估值精度能达到容许误差标准，能很好拟合实测资料，是一种较好的P—Ⅲ型分布参数估计方法.
参 考 文 献 1 金光炎. 论水文频率计算中的适线法.水文，1992，(4).
2 从树铮，等. 水文频率计算中参数估计的统计试验研究.水利学报，1980，(3).
3 陈守煜. 系统模糊决策理论与应用.大连：大连理工大学出版社，1994.
4 陈守煜. 模糊水文学与水资源系统模糊优化原理.大连理工大学出版社，1990.
5 金光炎，等. Γ分布保证率修正值参数B的确定.水文，1991，(6).
6 刘光文. 皮尔逊Ⅲ型分布参数估计.水文，1990，(4，5).
7 刘治中. 数值积分权函数法推求P—Ⅲ型分布参数.水文，1987，(7).
8 宋德敦，等. 概率权重矩法及其在P—Ⅲ型分布中的应用.水利学报，1988，(3).
9 邱 林，等. 水文水资源系统模糊集分析理论及应用.郑州：黄河水利出版社，1995.
Weighted optimum curve－fitting method for estimating the
parametersof Pearson type－Ⅲ distribution
Abstract The weighted optimum curve－fitting method for estimating the parameters of Pearson type－Ⅲ distribution is presented in this paper. This method can increase the stability of calculation and reduce the influnce of errors in hydrologic data on estimation of the parameters as well as utilizing information in real hydrologic data.The method is tested by safety standard of hydraulic engineering with a great deal of ideal and real hydrologic data. The results show that the method of least squares criterion is stable and its precision satisfies the safety standard for hydraulic engineering.
关键词：煤层瓦斯；吸附常数；分布参数；极大似然估计；极大似然原理 文献标识码：A
中图分类号：O647 文章编号：1009-2374（2015）36-0145-02 DOI：10.13535/ki.11-4406/n.2015.36.072
在科学技术水平不断发展的时代背景下，极大似然估计算法对于明确产品的可靠性来说具有积极的作用，本研究据此而展开。首先从某矿区进行煤块采摘取样，并随机分成五组数据，用容器测量法将这五组样本放在同样的容器下，并释放不同的压力，进行瓦斯吸附实验，利用极大似然估计方法对常数的平均数和差值进行研究，并将最接近的数值作为煤层吸附的参考值，为了确保实验结果的准确性，本次实验从煤层进行了全方面的随机取样。
1 煤层瓦斯吸附常数分布参数的主要特点
首先，煤层瓦斯是成煤反应的过程，高压高温是其发展的前提条件，经过物理和化学的双层作用，在这些客观环境因素的影响下形成的一种具有腐蚀性的有机物，煤层在形成过程中并不是很均匀，最后导致瓦斯在煤层释放出来时是通过散布和渗出等途径进行采取，如果在采取的过程中没有良好的理论知识和实际操作作为基础很容易发生火灾等意外。目前，实验室进行相关的吸附常数测量是了解煤层瓦斯吸附数情况的主要手段和途径，在实验的过程中采取的煤层可能根据分布点的不同，导致吸附常数也不相同，严格来说，这是因为吸附常数主要包括体积和压力常数，这两点的数值是不固定的，所以在研究中应该严格按照相关的瓦斯定义和规律进行计算，由于我国对相关内容的研究并不深入，导致研究受到了一定的阻碍，容易受到环境影响，如温度、瓦斯中的水分或者是电磁环境等方面，由此可见，现阶段的煤层瓦斯吸附常数分布参数的极大似然估计在实验过程中操作严谨。
2 煤层的瓦斯极大似然估计以及吸附性质参量
压力和体积参数是煤层瓦斯吸附的主要内容，将体积和压力分别用简写的英文字母VL和PL表示，并将我们在这次研究中用大写字母AB进行标示，采取正向分布的方法描述，即AI（i=1，2，3，…，n）和BI（i=1，2，3，…，n），根据概率，从而计算出样本的平均值和差值，在本次试验中，共记录了五组数据，即n=5，根据极大似然估计方法，得到吸附参数。具体的实验参考步骤和方法如下：
第一，分组，将采取好的煤层瓦斯进行一定数量的分组，在本次试验中共分为五组。
第二，测试，将分好数量的每一组煤块，放置在不同的压力下进行吸附测试，从而得到整体的煤块吸附量值即压力P和体积V的数值。
第三，计算，通过AI，BI（I=1，2，3，…5）样本参数，根据测量方法，进行极大似然估计，从而得出相近的数值。
第四，估计，最后将三步骤中的煤块瓦斯的平均值和平均方差的近似值作为参考数值。
3 煤块瓦斯的吸附值测量途径
实验室测量瓦斯吸附数值主要有三种方法，通常被我们所使用的是重量和容器法，这次煤样瓦斯吸附值测试主要采用后一种方法，该实验总共分为八个步骤进行，最后得到的结果即是压力和体积的极大似然估计值。具体操作可以分为以下八点：
第一，取样，严格遵守相关规定，采取没有受到外部环境干扰污染的煤块，并用完好无损的双层密封袋进行紧实的包装，按照GB/T211、GB/T212、GB/T217要求制作煤块，最后用破碎机器弄碎，用标准直径进行严格的挑选。通常情况下，煤块样本一般选择煤层随机厚度中块度大于20mm的样本。
第二，充气，将已经选好的样本放入吸附容器中，在压力不变的环境下，将瓦斯冲入容器中，并将容器密封好。此过程中对压力值进行恒定测量，确保压力始终介于临界值范围之内。
第三，保压，为了保证容器中的煤块样本能够充分地跟瓦斯进行融合，通常要保压在3个小时左右。
第四，测算，V1（容器中的瓦斯体积）-V2（未使用的瓦斯体积）=V3（煤块吸附的瓦斯体积）。
第五，改变压力，重复第二步和第四步。
第六，回归压力体积常数值，将吸附量设为S的话，那么与瓦斯压力P之间的计算公式为PVL-SPL-SP=0，从而计算出公式中的体积和压力参数值。
第七，样本改变，并回到第二步和第六步。
第八，将煤层的压力和体积参数进行极大似然
估计。
在进行以上步骤之后，将本组样本的Langmuir体积常数与压力常数设定为随机变量值进行观测并记录，假设该随机变量值为Y和Z，则当取值Zi，Yi（i=常数：1～6）。之后用观测值对随机变量InZ和InY进行极大似然估计，得出InZ与InY的平均值与平均方差。
上述结果能够有效估计出煤块的Langmuir体积与压力常数分布参数下的极大似然值，然而一个不容忽视的问题是，整个煤层的吸附常数不可能仅通过几组样本来进行估计，因此与所得的值的准确性亦有待进一步研究。而要想使得估计值准确，则所需之样本数量必然要远远大于样本组。但鉴于样本采集存在一定的困难，同时实验周期亦难以在短时间内完成，成本也较高，因此本研究仅将煤块的吸附常数所估计出的近似值判定为是对煤层吸附常数的估计值。
4 煤块的瓦斯性质测试结论
将样本采取100～200g的数量进行分组实验，除掉所得结果偏差较大的数值，仅保留正常结果样本，将实验的压力值和体积值进行显示，得到如下数据：
在瓦斯压力分别为0.5MPa、1.3MPa、2.1MPa、3.1MPa和4.3MPa下时，第一组分组的瓦斯吸附量分别为6.32cm3/g、11.9cm3/g、14.8cm3/g、16.1cm3/g和17.2cm3/g；第二组瓦斯吸附量分别为7.03cm3/g、12.3cm3/g、15.8cm3/g、16.9cm3/g和17.1cm3/g；第三组瓦斯吸附量分别为7.1cm3/g、12.8cm3/g、16.1cm3/g、17.2cm3/g和18.9cm3/g；第四组瓦斯吸附量分别为6.12cm3/g、10.8cm3/g、14.4cm3/g、16.1cm3/g和16.9cm3/g；第五组瓦斯吸附量分别为7.1cm3/g、10.4cm3/g、14.2cm3/g、16.4cm3/g和17.1cm3/g。
第一组至第五组Langmuir体积的随机变量InZ值分别为3.115、3.110、3.230、3.071和3.131，得出均值为3.21，标准差值为6.12×10-2；第一组至第五组Langmuir压力的随机变量InY值分别为0.12、0.211、0.225、0.423和0.432，均值为0.312，标准差为1.13×10-1。
在不同样本组的Langmuir体积与Langmuir压力常数值进行估计之后所得出的随机变量InZ和InY的平均值与标准差的极大似然估计值亦得出相应的结果，
5 讨论
第一，将直径不同的煤块在相同压力下进行瓦斯吸附研究，样本总数为50个，将不成功的实验除出去以后，成功的实验为32个，那么实验数据的整体利用率为64%，将这五组煤炭进行具体的结果分析和研究，并从中调取吸附数值，并作为煤块观测的估计数值。但是不同的采集地点会出现不同的常数值，为了保证煤块的瓦斯常数离散性能够更好地反映出煤层瓦斯吸附常数的离散性进而保证煤层的瓦斯吸附常数更加准确，应当以多处煤层作为样本采集地进行多处样本采集。
第二，样本的离散数值差距不大，也就是说吸附数值的离散性不高，在一定程度上来说，煤块瓦斯吸附数值的离散性与煤层瓦斯吸附数值的离散性不划等号，两者没有什么必然联系，发生这种情况的原因大部分是因为煤块的采集地点不够全面，导致数值测量不准确，在实际研究中为了能够保证数据的有效性，应该加大采取每块瓦斯的面积。
第三，由于此类实验的进行需要花费大量的精力和成本，此次研究只对范围内的每块瓦斯吸附特性进行了研究。
第四，考虑到煤块瓦斯的生成环境和影响因素，忽略了其他温度和湿度的影响，因此技术手段有待进一步提升。
第五，通过先进科学技术和手段的运用，大大减少了实验时间，能够对同一份样品在不同压力下的平衡吸附量进行测定，提高了实验效率。
6 结语
现阶段，极大似然估计方法是我国研究煤层瓦斯吸附常数分布参数的一种重要方法，由于煤层瓦斯的成长环境的关系，具体的实验通常都是在实验室中进行，为了确保实验结果的准确性应该进行全方面的煤块取样。本文以容器法作为测量瓦斯吸附数值的主要手段，并将抽取的煤块随机的分成五组，准确地记录了研究的数值，并客观地评价了煤块瓦斯吸附常数的离散度，最后对抽取样本中的体积系数和压力系数进行了简要的探析，希望对以后的煤层瓦斯吸附常数分布参数的极大似然估计研究提供参考性意见。
参考文献
[1] 张力，邢平伟.煤体瓦斯吸附和解吸特性的研究[J].江苏煤炭，2000，（4）.
[2] 张晓静，温百根，赵耀江.煤层瓦斯吸附能力实验研究[J].太原理工大学学报，2008，（z2）.
[3] 辛海会，徐立华，杜欣.煤层瓦斯含量的实验室测定[J].工业安全与环保，2009，（12）.
[4] 曹林.高压吸附下的瓦斯放散初速度研究[J].煤矿安全，2004，（9）.
[5] 曹林，仇海生.碎屑状煤芯瓦斯解吸规律研究[J].中国矿业，2007，（12）.
[6] 张占存.水分对煤的瓦斯吸附特性影响的实验研究
[D].煤炭科学研究总院，2004.