底栖生物完整性评价指数结构方程模型的构建方法

技术领域：本发明涉及一种模型的构建方法，尤其是一种底栖生物完整性评价指数结构方程模型的构建方法。

背景技术：生物完整性指数(index of biotic integrity，IBI)是指可定量描述环境状况特别是人为干扰与生物特性之间的关系的一组敏感性生物指数。IBI评价体系是通过对生态系统中某一生物类群所构成的群落(如底栖动物)的物种组成、多样性及功能结构方面进行分析，将其和相应的标准体系比较，然后根据指标评估所观测生态系统的健康状况。鱼类作为生物完整性研究最早应用的对象，在河流监测评价中具有很多优点，目前为止国际众多学者已经对利用鱼类进行生物完整性评价进行了广泛的研究、论证与应用。

1、B-IBI(Benthic Index Of Biotic Integrity)指数的建立:

B-IBI指数的建立方法包括5个重要步骤:

(1)样点底栖动物资料、水体理化数据和栖境质量资料收集;

(2)根据分类学和生态学特征确定候选生物学指数;

(3)根据各生物指数值在参照点或所有样点中的分布，进行生物指数筛选，指数值与分值之间的转换;

(4)建立B-IBI指标体系;

(5)应用B-IBI评价水体的生物完整性程度，与理化指标和栖境质量进行相关性分析。

2、候选生物学指数确定:

构成B-IBI所选用的生物指数必须具有敏感性、可比性、准确性和实用性，拟选用以下20个生物学指数作为候选指数:

反映群落丰富度和个体数量比例的指数——总分类单元数，EPT分类单元数，敏感类群分类单元数，水生昆虫分类单元数，甲壳动物+软体动物分类单元数，摇蚊分类单元数，优势分类单元%，前3位优势分类单元%，毛翅目%，蜉蝣目%，颤蚓%，襀翅目%，摇蚊%，(甲壳动物+软体动物)%，无足类群%；

反映营养级组成的指数——捕食者%，滤食者%；

与生物耐污程度有关的指数——耐污类群%，敏感类群%；

栖境质量参数——粘附者%。

3、现有IBI指标的计算

较多的有传统的1、3、5赋值完整性的评分计算。

1)1、3、5赋值法

是目前应用最多的评分方法。其对各参数的标准化方法为：对各指标测得的实际值在最低到最高范围内三等分，分为3个区域，最好的等级区域记为5，最差的等级区域记为1，中间的为3。每个站点IBI值即为该点各参数标准化记分后的值之和。

2)比值法

比值法中各参数的标准化指数模式分两种情况。

①随干扰增强而值减小的指标，其标准化指数模式为:

$>P_{\mathrm{ij}}=\frac{O_{\mathrm{ij}}}{S_{i95}}\times 100$>

式中，P_ij和O_ij，分别为第j个采样点第i个指标的标准化指数和原始观测值;S_i95为第i个指标的标准化阈值，取全部采样点中第i个指标95th百分位的原始观测值;i=1，2，…，m;j=1，2，…，n。

②随干扰增强而值增大的指标，其标准化指数模式为:

$>P_{\mathrm{ij}}=\frac{{\max O}_{\mathrm{ij}}-O_{\mathrm{ij}}}{{\mathrm{naxO}}_{\mathrm{ij}}-S_{i5}}\times 100$>

式中，max O_ij为第i个指标在n个采样点中的最大值;S_i5为第i个指标的标准化阈值，取全部采样点中第i指标5th百分位的原始观测值;i=1，2，…，m;j=1，2，…，n。

各站点的IBI值是该点各参数标准化指数的均值。即:

$>I_j=\frac{1}{m}\underset{i=1}{\overset{m}{\Sigma }}{P}_{\mathrm{ij}}$>

式中，I_j为第j个采样点的IBI值。

4、生物学指数筛选:

1)指数值分布范围分析,利用参照点数据资料计算各候选生物指数值，分析各指数值的分布范围;

2)判别能力分析,经上述筛选后，符合条件的指数，利用受损点资料计算它们的值,采用箱线图法分析各指数值在参照点和受损点之间的分布情况,主要根据Barbour评价法。

5、相关性分析:

通过Pearson相关性，根据相关系数的大小决定生物指数间所反映信息的重叠程，如果两个指数间的相关系数|r|>0.75,表明两个指数间所反应信息大部份是重叠的，选其中一个,经过上述3步分析进而确定一组构成B-IBI指标体系的生物指数。

6、生物学指数记分标准的建立:

采用指数分值计算的方法建立B-IBI指标体系评价标准。参照点为B-IBI值分布的25%分位数法，如果样点的B-IBI值大于25%分位数值，大于25%分位数值则表示该样点受到的干扰很小，生物完整性是“高”的等级。小于25%分位数值的分布范围，可以4等分，分别代表不同的完整性程度。

7、B-IBI指标体系的评价标准:

根据现状偏离未受干扰河流特有种的程度来划分生物完整性的等级。并把生物完整性的现状划分成5个等级即:高、好、适度、差、恶化。“高”的生物完整性状态是完全或者计划全部与没有受到干扰的参考点情况一致。“好”的生物完整性有着重要的但是轻微偏离没有受到干扰的状态的特征。在“适度”的生物完整性层次，所有的标准都表现出较强的偏离没有受到干扰的状态。“差”的生物完整性则受到很强的偏离，而“恶化”则是极度偏离。采用指数分值计算的方法建立B-IBI指标体系评价标准。参照点B-IBI值分布的25%分位数法，如果样点的B-IBI值大于25%分位数值，大于25%分位数值则表示该样点受到的干扰很小，生物完整性是“高”的等级。小于25%分位数值的分布范围，可以4等分，分别代表不同的完整性程度。一般而言，“高”和“好”的生态系统完整性状态应当成为环境管理者的管理目标，如果生态系统完整性的状态比这个更差的话，就有必要采取措施进行恢复。

8、现有方法存在的问题：

(1)现有方法是建立在常规统计分析的原理之上，虽然其容许统计分析模型中的因变量含有测量误差，但需要假设自变量是没有误差的，而实际上自变量的观察值是存在误差的，因此，建立在传统统计理论上得B-IBI指数算法高估了观测变量的真正变异量。

(2)B-IBI指数涉及到多个因变量，应用传统统计原理计算的回归系数或因素间的路径系数其实是对每个因变量的逐一计算，其貌似对多个因变量同时考虑，实际上在计算对某一个因变量的影响或关系时，都忽略了其他因变量的存在及其影响，无法考虑多个因变量的交互作用。

(3)应用传统统计原理构建的B-IBI指数法，无法了解潜变量之间的相关关系，由于没有考虑到各因子内结构与其他因子的存在的变化情况，所以传统B-IBI算法无法处理一个指标从属多个因子的问题已经高阶影响因子等比较复杂的从属关系问题。

发明内容：针对上述现有技术的缺点，本发明的目的在于提供一种诊断水体沉积物毒性的最优方法，该方法能高效、快速、准确的完成水体沉积物毒性诊断。

为实现上述目的，本发明采用的技术方案是：底栖生物完整性评价指数结构方程模型的构建方法，分成四个步骤进行，具体如下：

(1)B-IBI结构方程模型建构：

①建立观测变量与潜变量之间的关系；

②建立各个潜变量之间的关系；

③确定限制因子负荷及因子相关系数参数的数值。

(2)B-IBI结构方程模型拟合：

求出(1)步骤中建立的模型参数的多套估计值，从中找出能使得模型隐含的协方差矩阵与样本协方差矩阵的“差距”最小。

(3)B-IBI结构方程模型评价：

①检验B-IBI结构方程模型的解是否适当，包括算法的迭代估计是否收敛，各个参数估计值是否在合理范围内；

②检验参数与预设模型的关系是否合理。数据分析可以出现一些预期以外的结果，但各个参数绝不应出现一些互相矛盾，与先验假设有严重冲突的现象；

③验视多个不同类型的整体拟合指数，如非范拟合指数(Non-Normed Fit Index，NNFI)、拟合指数(Comparative Fit Index，CFI)、近似误差指数(Root Mean Square Error of Approximation，RMSEA)和χ²等，以此衡量模型的拟合程度。

(4)B-IBI结构方程模型修正：

①依据结构方程理论和有关假设，提出一个或数个合理的先验模型；

②检查潜变量与指标之间的关系，建立测量模型；

③若模型中含有多个因子，可以循序渐进地，每次检查含两个因子的模型，确立测量模型部分的合理后，最后再将所有因子合并成预设的先验模型，作一个总体检查。

④对每个模型，检查标准误、t值、标准化残差、修正指数、参数期望改变值、χ²及各种拟合指数，并据此修改模型并重复第③和④步骤。

⑤最后的模型将依据某一个样本数据修改而成，用另一个独立的样本交互确定。

B-IBI结构方程模型算法主要被分解成模型建构、模型拟合、模型评价、模型修正四个部分，通过合理、严谨的算法设计以解决之前现有IBI指数算法存在的问题。

附图说明：

图1为依据底栖生态系统理论构建的B-IBI结构方程机理模型。

具体实施方式：

如图1所示：假设构成B-IBI所选用的生物指数如下：

反映群落丰富度和个体数量比例的指数——总分类单元数，EPT分类单元数，敏感类群分类单元数，水生昆虫分类单元数，甲壳动物+软体动物分类单元数，摇蚊分类单元数，优势分类单元%，前3位优势分类单元%，毛翅目%，蜉蝣目%，颤蚓%，襀翅目%，摇蚊%，(甲壳动物+软体动物)%，无足类群%；

反映营养级组成的指数——捕食者%，滤食者%；

与生物耐污程度有关的指数——耐污类群%，敏感类群%；

栖境质量参数——粘附者%。

其中：

总分类单元数，EPT分类单元数，敏感类群分类单元数，水生昆虫分类单元数，甲壳动物+软体动物分类单元数，摇蚊分类单元数，优势分类单元%，前3位优势分类单元%，毛翅目%，蜉蝣目%，颤蚓%，襀翅目%，摇蚊%，(甲壳动物+软体动物)%，无足类群%，捕食者%，滤食者%，耐污类群%，敏感类群%，粘附者%为内因观察变量Y₁、Y₂、......、Y₇；

相应的内因观察变量的误差项为ε₁、ε₂、......、ε₂₀；

反映群落丰富度和个体数量比例的指数、反映营养级组成的指数、与生物耐污程度有关的指数、栖境质量参数为一阶内因潜在变量η₁、η₂、η₃、η₄；

相应的因潜在变量的误差项为ζ₁、ζ₂、ζ₃、ζ₄；

内因观察变量与一阶内因潜在变量之间的负荷因子为λ₁₁、λ₁₂、λ₁₃、λ₁₄、λ₁₅、λ₁₆、λ₁₇、λ₁₈、λ₁₉、λ₁₁₀、λ₁₁₁、λ₁₁₂、λ₁₁₃、λ₁₁₄、λ₁₁₅、λ₂₁₆、λ₂₁₇、λ₃₁₈、λ₃₁₉、λ₄₂₀；

IBI指数为二阶内因潜在变量ξ1；

二阶内因潜在变量与一阶内因潜在变量之间的负荷因子为γ₁₁、γ₂₁、γ₃₁、γ₄₁。

依据机理模型构建IBI指数的测量方程矩阵如下：

$>\left(\begin{array}{c}{Y}_1\\ Y_2\\ Y_3\\ Y_4\\ Y_5\\ Y_6\\ Y_7\\ Y_8\\ Y_9\\ Y_{10}\\ Y_{11}\\ Y_{12}\\ Y_{13}\\ Y_{14}\\ Y_{15}\\ Y_{16}\\ Y_{17}\\ Y_{18}\\ Y_{19}\\ Y_{20}\end{array}\right)=\left(\begin{array}{cccc}{\lambda }_{11}& 0& 0& 0\\ {\lambda }_{12}& 0& 0& 0\\ {\lambda }_{13}& 0& 0& 0\\ {\lambda }_{14}& 0& 0& 0\\ {\lambda }_{15}& 0& 0& 0\\ {\lambda }_{16}& 0& 0& 0\\ {\lambda }_{17}& 0& 0& 0\\ {\lambda }_{18}& 0& 0& 0\\ {\lambda }_{19}& 0& 0& 0\\ {\lambda }_{100}& 0& 0& 0\\ {\lambda }_{111}& 0& 0& 0\\ {\lambda }_{112}& 0& 0& 0\\ {\lambda }_{113}& 0& 0& 0\\ {\lambda }_{114}& 0& 0& 0\\ {\lambda }_{115}& 0& 0& 0\\ 0& {\lambda }_{216}& 0& 0\\ 0& {\lambda }_{217}& 0& 0\\ 0& 0& {\lambda }_{318}& 0\\ 0& 0& {\lambda }_{319}& 0\\ 0& 0& 0& {\lambda }_{420}\end{array}\right)\left(\begin{array}{c}\\ \end{array},\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\eta }_4\end{array}\right),+,\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ {ϵ}_4\\ {ϵ}_5\\ {ϵ}_6\\ {ϵ}_7\\ {ϵ}_8\\ {ϵ}_9\\ {ϵ}_{10}\\ {ϵ}_{11}\\ {ϵ}_{12}\\ {ϵ}_{13}\\ {ϵ}_{14}\\ {ϵ}_{15}\\ {ϵ}_{16}\\ {ϵ}_{17}\\ {ϵ}_{18}\\ {ϵ}_{19}\\ {ϵ}_{20}\end{array}\right)\right)---\left(1\right)$>

依据机理模型构建IBI指数的结构方程矩阵如下：

$>\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\eta }_4\end{array}\right)=\left(\begin{array}{c}{\gamma }_{11}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{14}\end{array}\right)\left[{\xi }_1\right]+\left(\begin{array}{c}{\zeta }_1\\ {\zeta }_2\\ {\zeta }_3\\ {\zeta }_4\end{array}\right)---\left(2\right)$>

设

$>y=\left(\begin{array}{c}{Y}_1\\ Y_2\\ Y_3\\ Y_4\\ Y_5\\ Y_6\\ Y_7\\ Y_8\\ Y_9\\ Y_{10}\\ Y_{11}\\ Y_{12}\\ Y_{13}\\ Y_{14}\\ Y_{15}\\ Y_{16}\\ Y_{17}\\ Y_{18}\\ Y_{19}\\ Y_{20}\end{array}\right),{\Lambda }_y=\left(\begin{array}{cccc}{\lambda }_{11}& 0& 0& 0\\ {\lambda }_{12}& 0& 0& 0\\ {\lambda }_{13}& 0& 0& 0\\ {\lambda }_{14}& 0& 0& 0\\ {\lambda }_{15}& 0& 0& 0\\ {\lambda }_{16}& 0& 0& 0\\ {\lambda }_{17}& 0& 0& 0\\ {\lambda }_{18}& 0& 0& 0\\ {\lambda }_{19}& 0& 0& 0\\ {\lambda }_{100}& 0& 0& 0\\ {\lambda }_{111}& 0& 0& 0\\ {\lambda }_{112}& 0& 0& 0\\ {\lambda }_{113}& 0& 0& 0\\ {\lambda }_{114}& 0& 0& 0\\ {\lambda }_{115}& 0& 0& 0\\ 0& {\lambda }_{216}& 0& 0\\ 0& {\lambda }_{217}& 0& 0\\ 0& 0& {\lambda }_{318}& 0\\ 0& 0& {\lambda }_{319}& 0\\ 0& 0& 0& {\lambda }_{420}\end{array}\right)\left(\begin{array}{c},\eta =\\ \end{array},\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\eta }_4\end{array}\right),,,ϵ,=,\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ {ϵ}_4\\ {ϵ}_5\\ {ϵ}_6\\ {ϵ}_7\\ {ϵ}_8\\ {ϵ}_9\\ {ϵ}_{10}\\ {ϵ}_{11}\\ {ϵ}_{12}\\ {ϵ}_{13}\\ {ϵ}_{14}\\ {ϵ}_{15}\\ {ϵ}_{16}\\ {ϵ}_{17}\\ {ϵ}_{18}\\ {ϵ}_{19}\\ {ϵ}_{20}\end{array}\right)\right)$>

则公式(1)可写成：

y＝Λ_yη+ε    (3)

设

$>\eta =\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\eta }_4\end{array}\right),\Gamma =\left(\begin{array}{c}{\gamma }_{11}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{14}\end{array}\right),\xi =\left[{\xi }_1\right],\zeta =\left(\begin{array}{c}{\zeta }_1\\ {\zeta }_2\\ {\zeta }_3\\ {\zeta }_4\end{array}\right)$>

则公式(2)可写成：

η＝Γξ+ζ    (4)

其中：

y为内因观察变量向量；

Λ_y为内因观测变量与一阶内因潜在变量之间的因子负荷矩阵；

η为一阶内因潜在变量向量；

ε为一阶内因潜在变量的误差向量；

ξ为外源潜在变量的向量；

Γ为外源潜在变量对内因观察变量影响向量；

ζ为因观察变量的误差向量。

设内因观察变量y之间的协方差矩阵为S(因为协方差矩阵为对称阵，所以只写出其下三角矩阵表示)，则

θ表示模型中全体未知参数组成的向量，表示θ的估计，由所研究的模型推出的总体协方差矩阵记为Σ(θ)，根据样本估计出参数后得到的协方差矩阵记为则有:

则S与Σ(θ)之间的差距函数如下:

F(S,Σ(θ))    (5)

F是观察样本共变数矩阵S与理论假设模式矩阵Σ(θ)间距离的数值，求解其要满足以下条件：

①F值最小化；

②F(S,Σ(θ))≥0；

③F(S,Σ(θ))＝0,若且为若Σ(θ)＝S；

④在S与Σ(θ)中，F(S,Σ(θ))是一个连续函数。

满足以上4个条件而获得的θ的一致性估计值即为IBI模型的指数参数值。

其求解方法选择广义最小二乘法(Generalized LeastSquares,GLS),则公式(5)可写成：

$>F_{\mathrm{GLS}}=\frac{1}{2}\mathrm{tr}\left\{{\left[(S-\Sigma \left(\theta \right))W^{-1}\right]}^2\right\}---\left(6\right)$>

其中

W^-1为残差矩阵的加权矩阵，为正定矩阵。

当W^-1=S^-1时，则：

$>\left(\begin{array}{c}{F}_{\mathrm{GLS}}=\frac{1}{2}\mathrm{tr}\left\{{\left[(S=\Sigma \left(\theta \right))S^{-1}\right]}^2\right\}\\ =\frac{1}{2}\mathrm{tr}\left\{{\left[(I-\Sigma \left(\theta \right))S^{-1}\right]}^2\right\}\end{array}\right)---\left(7\right)$>

算法通过以下参数确定模型的优劣：

(1)卡方检验

χ²＝(N-1)min{F_GLS}    (8)

其中：

N为样本的容量；

min{F_GLS}为F函数的最小值，其值为

对于GLS算法，min{F_GLS}这个距离的(N-1)倍可用于检验假设:H₀:Σ＝Σ(θ)。

(2)近似误差均方根(Root Mean Square Error of Approximation,RMSEA)

$>\mathrm{RMSEA}={\{\max \frac{{\chi }^2-\mathrm{df}}{N-1},0]/\mathrm{df}\}}^2---\left(9\right)$>

其中：

Df为卡方的自由度。

(3)非范拟合指数(Non-Normed Fit Index,NNFI)

$>\mathrm{NNFI}=\frac{{\chi }_N^2/{\mathrm{df}}_N-{\chi }_T^2/{\mathrm{df}}_T}{{\chi }_N^2/{\mathrm{df}}_N-1}---\left(10\right)$>

(4)比较拟合指数(Comparative Fit Index,CFI)

$>\mathrm{CFI}=1-\frac{\max [({\chi }_T^2-{\mathrm{df}}_T),0]}{\max [({\chi }_T^2-{\mathrm{df}}_T),({\chi }_N^2-{\mathrm{df}}_N),0]}---\left(11\right)$>

模型验证原则:

χ²值必须未达到显著性水平，即p值必须>0.10；

RMSEA值小于或等于0.05表示“模型良好”，0.05-0.08表示“不错的模型”，0.08-0.10表示“中度的模型”，大于0.08表示“不良适配模型”；

NNFI和CFI在0.9以上(越大越好)；

满足这些条件的IBI模式可以认为是一个好模型。