鱼类自然死亡系数(M)是渔业资源评估中重要的参数。由于补充、捕捞死亡率、生长等参数的影响，自然死亡系数也就成为最难估算的参数之一^[1]。渔业资源学文献中提出的关于自然死亡系数的估算方法，主要有标记-回捕法^[2-5]，渔获量曲线分析方法^[6]以及一些经验公式^[7-9]。然而，任何一种方法都不能得到很好的估算结果^[10]。QUINN和DERISO^[1]曾经阐述了6种估算自然死亡系数的常用方法，但是任何一种方法都存在很大的不确定性^[11]。我们在以前的工作中^[10]研究了用POPE^[12]提出的股分析模型求算M的方法，但是没有考虑模型的过程误差(process error)，并且观测误差(measurement error)也只讨论了正态分布的情况。

广义线性模型(GzLM)是一种基于最大似然法并且拥有一定系统构架的方法，如果模型的误差结构属于指数类型，那么模型的参数就可以通过GzLM估算求得^[13]。JIAO等^[13]讨论了正态，对数正态，泊松和伽马等分布的误差结构对亲体-补充量模型的影响。JIAO和CHEN^[14]将上述4种误差结构应用到了产量模型和连续种群模型当中，并且认为对数正态和伽马分布适用于产量模型，而伽马分布同样适用于连续种群模型。

本研究通过将GzLM应用于Pope的股分析模型探讨了正态，对数正态和伽马3种分布的模型和数据的误差结构对M估算的影响，用皮尔森残差(Pearson residuals)对估算结果进行了同质性检验^[14]。这3种误差结构也同样被应用到利用GzLM估算黄海鳀鱼M的过程中。

1.   方法和数据

1.1   模型

1.1.1   Pope股分析模型

POPE^[12]提出的股分析模型为：

模型当中N_{t, a}是t年初a龄鱼的资源量，N_{t+1, a+1}是下一年初a+1龄鱼的资源量，C_{t, a}是t年a龄鱼的渔获量，所有数量的单位都是尾数。N的数据通过资源调查获得，而C则是统计数据，因此N和C没有相关性。M是自然死亡系数。

1.1.2   广义线性模型(GzLM)

GzLM由3部分组成^[15]。第一部分是随机部分Y，它是一个由n个独立的成分组成的观测向量，其平均值为μ。第二部分是系统部分，根据一些未知参数β₁，β₂，…，β_p对向量μ进行具体的解释。同时给出一个线性的预测值η，并且$\eta=\sum\limits_{j=1}^p X_i \beta_j$，其中X为模型矩阵或者观测数据Y的共变向量(covariates)。第三部分是随机部分和系统部分的链接，通常表示为η=g(μ)，其中g即为链接函数^[15]。如果用i来表示观测值，那么模型的系统部分就可以写成$\eta = \sum\limits_{j = 1}^p {{X_{ij}}} {\beta _j}$，其中X_ij表示第j个共变向量的第i个观测值。

Y的似然函数可以表示成：

其中f是概率分布函数(probability distribution function，PDF)，θ是需要估计的参数向量。对数似然函数可以表示为：

Y代表因变量，X代表自变量，而θ代表参数。

链接函数通常的选择包括恒等式，对数关系式，倒数关系式以及指数关系式等等^[14-16]。GzLM能够灵活地结合不同的链接函数而不必局限于传统的形式。

1.1.3   运用GzLM估算自然死亡系数

本研究采用了恒等式作为链接函数，误差结构为正态，对数正态和伽马分布。当使用对数正态分布时，方程(1)可以改写为：

其中ε服从正态分布N (0, σ²)；lnN_{t+1, a+1}是因变量，N_{t, a}和C_{t, a}是GzLM中的偏移量(offsets)。由数字1组成的与N_{t, a}相同维数的矩阵充当GzLM中的自变量。对于另外2种误差分布，N_{t+1, a+1}是因变量，而N_{t, a}和C_{t, a}是自变量，θ代表方程(1)中的e^-M和$e^{-\frac{M}{2}}$。这样模型中的参数可以通过计算似然函数的最大值求得。

1.1.4   模拟设计

讨论了3种分布的数据误差结构，包括正态、对数正态和伽马分布。随机变量根据PRESS等^[17]提供的方法产生。误差结构服从以上3种分布的各年龄的资源量和渔获量模拟数据分别由下列的方程(5)和(6)产生。

其中F代表捕捞死亡系数。

参数的估计值与真实值之间的差异用相对估计误差(relative estimate error，REE)来衡量

其中M_i^*表示第i次模拟得到的自然死亡系数的估计值，M是自然死亡系数的真实值，n代表模拟的次数。REE的数值越小表示M的估计值越精确^[13]。

1.1.5   皮尔森残差和同质性

皮尔森残差

被广泛用来分析模型误差结构的同质性，上式中y和u分别表示变量的观测值和估计值，V(u)是估计值的标准偏差^[14-15]。由于文中所用的样本比较小(8年×8龄=64，见下文)，用肉眼即可判断皮尔森残差相对变量的估计值所作散点图的同质性。当残差值呈现出关于x轴对称的左三角或右三角形图案时，就认为模型使用了不适当的误差结构^[14]。

1.2   数据

1.2.1   模拟数据

由方程(5)和(6)产生包括8龄和8年的N和C的模拟数据矩阵。假定M固定不变等于0.4，F等于0.6。表 1给出了一个世代8个年龄的模拟数据。对每一种误差结构我们假定了5个水平的白色噪音(即变异系数，CV=1%，5%，10%，20%和30%)，并且渔获量的噪音水平等于资源量的一半(通常渔获量数据比资源量数据更加准确)。对不同水平的CV每组模拟数据重复计算1 000次以得到稳定结果。

表  1  捕捞死亡系数(F)等于0.6 a^-1自然死亡系数(M)等于0.4 a^-1条件下利用方程(5)和(6)得到的模拟数据[没有白色噪音(CV)]

Table  1.  Example of the simulated data using Eq. (5) and (6) [without white noise(CV)] when fishing mortality (F) was 0.6 a^-1, natural mortality coefficient(M) was 0.4 a^-1

	年龄/a age
	1	2	3	4	5	6	7	8
资源量 abundance	109 663.3	40 342.9	14 841.3	5 459.8	2 008.6	738.9	271.8	100.0
渔获量 catch	41 592.3	15 300.9	5 628.9	2 070.8	761.8	280.2	103.1	37.9

下载: 导出CSV 

| 显示表格

为了研究不同生长状况和M对估算结果的影响，构造了2个鱼类种群：一个M相对较小(0.2)但寿命较长(11龄)，另一个则相反，M较大(0.8)寿命较短(4龄)^[18]。在进行这2个鱼类种群的模拟分析时，F假定等于0.3。

1.2.2   黄海鳀鱼(Engraulis japonicus)渔业真实数据

鳀鱼是一种生长快，寿命短的鱼种。本研究所用的是16年(1987~2002)4个年龄(1~4)的黄海鳀鱼的数据^[19]，其中包括资源量和渔获量(ZHAO等^[19]的表 1和表 2)。但因未知1997和1998年的资源量以及1987和1988年的渔获量数据，故本研究仅使用了8年的数据(1989~1996)。

表  2  不同模拟数据的白色噪音(CV, CV_N=2CV_C)水平和不同广义线性模型(GzLM)及数据的误差结构条件下自然死亡系数(M)的相对估计误差(REE, %)

Table  2.  Summary of the relative estimate error (REE, %) of the estimated natural mortality coefficient (M) for different white noises (CV) of simulated data (CV_N=2CV_C) and error structures in generalized linear model (GzLM) and simulated data

误差结构 error structure			资源量白色噪音 white noise of abundance (CVN)
GzLM	数据 data		1%	5%	10%	20%	30%
正态 normal	正态 normal		2.551	5.626	9.684	21.239	36.816
	对数正态 lognormal		17.424	70.366	123.789	168.306	255.133
	伽马 gamma		3.887	7.262	16.513	33.373	50.460
对数正态 lognormal	正态 normal		3.945	5.807	9.947	29.482	41.941
	对数正态 lognormal		7.909	31.448	41.828	53.649	60.307
	伽马 gamma		2.482	2.882	5.471	9.125	11.035
伽马 gamma	正态 normal		2.947	5.302	9.787	21.218	38.474
	对数正态 lognormal		18.867	72.724	159.738	264.500	289.795
	伽马 gamma		3.757	5.912	12.716	27.832	44.471

下载: 导出CSV 

| 显示表格

2.   结果

2.1   不同误差结构和不同白色噪音水平下的模拟分析

图 1为3种模型误差和数据误差条件下皮尔森残差对资源量估计值所作的散点图。当GzLM的误差结构为对数正态分布时，3种分布的数据误差结构均得到了同质的皮尔森残差(图 1-b)。而当GzLM的误差结构为正态和伽马分布时皮尔森残差值是异质的(图 1-a和图 1-c)。同一种GzLM误差结构条件下不同数据误差的计算结果表明，与正态和伽马分布相比由对数正态分布的数据误差所作的皮尔森残差散点图的异质性最为明显(图 1)。对3种GzLM误差分布而言，随着数据噪音水平的增加，对应的皮尔森散点图的异质性也随之增加，然而与另外2种分布相比，对数正态分布的GzLM误差结构对数据的CV具有较小的敏感度(图 1)。表 2显示，当数据的噪音水平增加时M估计值的相对估计误差(REE)值也随之增大；当数据的噪音水平低于大约10%时，数据误差为正态和伽马分布的2种情况都能够得到可靠的M估计值。不同数据误差结构和GzLM误差结构的组合所得到的结果存在差异，正态分布的数据误差适合于正态和伽马分布的GzLM误差结构，而伽马分布的数据误差则更适合于对数正态分布的GzLM误差结构。

图  1  3种模型误差条件下的皮尔森残差(Pearson residuals)

a.正态分布; b.对数正态分布; c.伽马分布
1-3列分别代表数据的正态，对数正态和伽马分布；1-5行分别代表资源量的白色噪音(CVN)水平1%，5%，10%，20%和30%(CV_N=2CV_C)

Figure  1.  Pearson residuals diagnostic plots for three error structures in the models

a. normal distribution; b. lognormal; c. gamma in GzLM
The 1st to 3rd columns represent normal, lognormal and gamma distributions in the simulated data; the 1st to 5th rows represent white noise levels of abundance data (CVN): 1%, 5%, 10%, 20% and 30% (CV_N=2CV_C).

下载: 全尺寸图片 幻灯片

表 3总结了2个不同寿命鱼类种群M估计结果的REE值。由于对数正态分布的GzLM误差结构得到了最同质性的皮尔森残差，并且在该条件下伽马分布的数据误差得到了最好的M估计值(表 2)，所以在对这2个种群的模拟计算过程中我们采用了这2种误差结构的组合。2个种群M估计值的REE都随着数据CV的增加而变大。所得的短寿命种群M估计值的REE值小于长寿命种群(表 3)，但是两者之间的差异不显著(P=0.161)。

表  3  2个模拟种群在不同模拟数据的白色噪音条件下(CV, CV_N=2CV_C)自然死亡系数(M)估计值的相对估计误差(REE, %)

Table  3.  The relative error estimation (REE, %) of the estimated natural mortality coefficient (M) of two simulated populations for different white noises (CV) simulated data (CV_N=2CV_C)

种群 population	1%	5%	10%	20%	30%
长寿命小M long-lived with small M	2.448	5.304	9.570	14.470	17.195
短寿命大M short-lived with large M	1.426	2.276	3.962	6.783	8.020

下载: 导出CSV 

| 显示表格

2.2   黄海鳀鱼自然死亡系数的估算

本研究应用股分析(CA)模型分别求算了3种GzLM误差结构条件下不同年龄黄海鳀鱼的M。表 4列出了不同误差结构下黄海鳀鱼各龄的M以及2003年ZHAO等^[19]所估算的结果。从中可以看到研究所得到的1龄鱼的M小于ZHAO等^[19]的估计，而两者2龄鱼M的结果则比较接近。当使用全部观测数据来计算3龄鱼的M时本研究得到了较大的估计值，其原因可能是由于4龄鱼资源量的观测数据中存在的异常值(本文资源量异常值定义为大于低年龄组或者小于低年龄组1%的资源量观测值)。用删除2个异常值(1993和1995年4龄鱼的资源量)后的观测数据计算3龄鱼的M，此时得到了较为合理的结果。在所假设的3种GzLM误差结构当中对数正态分布得到了最好的估计结果。在这里应该指出，ZHAO等^[19]在计算3龄鱼M的时候也删除了1993年4龄鱼的资源量数据(ZHAO个人通讯)。比较3种GzLM误差结构下删除异常值前后3龄鱼M的计算结果，发现删除异常值后在对数正态分布的误差结构条件下计算的结果与ZHAO等^[19]所给出的结果最为接近。这表明对数正态分布的误差结构较另外2种更适合于POPE ^[12]的股分析模型。另外，也提醒我们即使计算过程中使用了适当的误差结构仍要注意异常值的存在。

表  4  不同广义线性模型误差结构条件下黄海鳀鱼自然死亡系数(M, a^-1)估计值

Table  4.  Estimated natural mortality coefficient (M, a^-1) of anchovy in the Yellow Sea for different error distribution in generalized linear model (GzLM)

误差结构 error structure	年龄 age			删除异常值后3龄的M M of age 3 without outliers
	1	2	3
正态 normal	0.26	0.50	1.68	1.19
对数正态 lognormal	0.15	0.56	2.04	0.94
伽玛 Gamma	0.22	0.52	2.39	1.23
ZHAO等[19]	0.09	0.45	0.92
注：最后一列为删除异常值(1993和1995年4龄鱼的资源量)后3龄鱼的M Note: The last column is the estimated M of age 3 without outliers (abundances of age 4 for 1993 and 1995).

下载: 导出CSV 

| 显示表格

3.   讨论

不同分布的GzLM和数据误差结构得到了精确程度不同的M的估计值。在应用Pope股分析模型求算M时，对数正态分布的GzLM误差结构得到了较好的结果(图 1)，而正态和伽马分布的数据误差结构则好于对数正态(表 2)。这说明渔业种群评估的效果既受到渔业数据的影响也要受到评估模型的影响，也就是说数据和模型共同决定评估结果的好坏^[14]。因此，在每一次评估过程当中应该尽可能多地检验假设的误差分布类型，并分析找出适当的数据和模型的误差结构。

文中构造的具有不同寿命的2个鱼类种群的模拟结果显示，短寿命鱼类(M较大)得到了较好的M的估计值(表 3)，这可能主要是由于M自身数值的影响，因为在模拟过程中我们对2个种群采用相等的F。此外，短寿命种群拥有较少的观测数据，这样每年的所累积CV就会较小，这样也能够提高估计结果的精确度。增加长寿命种群所遭受的F到0.9(使2个种群的总死亡率都达到1.1)，其M估计结果的REE值也随之变大。然而，F增大前与增大后所得到的M估计值之间的差异并不显著(P=0.658)，但是这个差异性系数(P)值要大于F不变时2个拥有不同M的种群所求得的M估计值的P值(P=0.161)。这些模拟结果表明，M自身数值的大小对其估计结果质量产生的影响要大于F，因此当使用Pope股分析模型求算M时，应该更加注意参数M自身的性质。

黄海鳀鱼1龄鱼M的计算结果比我们设想的值要小，原因可能是冬季在黄海进行声学调查期间1龄鱼还没有完全补充，导致对鱼群观测不完全^[19]。2龄鱼则得到了比较合理的M值，并且3种GzLM误差结构条件所计算的结果比较接近。3龄鱼M的计算结果不够精确，这可能是调查过程中4龄鱼的数量少而且波动大的原因所导致。删除4龄鱼资源量的异常值之后便得到了较为合理的结果，尤其是对数正态分布的误差结构得到了M最好的估计值。这也又一次证实了文中的模拟分析结果(图 1、表 2)，适当的误差分布能够改善计算的结果，另外不可靠的数据或者异常值将会降低计算结果的精确度。

本研究讨论了使用Pope股分析模型求算鱼类M时不同分布的GzLM和数据误差结构对求算结果产生的影响。使用GzLM和皮尔森残差能够较好地分析和确定适当的模型和数据的误差结构，进而改善参数的求算质量。