CN112115414B - 一种广布种分布范围的预测方法 - Google Patents

Abstract

本发明涉及物种分布研究领域，具体涉及一种广布种分布范围的预测方法，包括根据二截尾零膨胀泊松模型：对广布种各点分布概率进行预测的步骤。通过采用零膨胀模型对广布种的分布范围进行预测，零膨胀模型中的清楚地刻画了观测数据中的零样本是因为某地确实没有该物种存在，还是因为当次观测没有观测到该物种，从而可以从本质上对上述两种情况进行区分，有效解决了广布种观测数据不平衡、零数据过多的问题，使得对广布种分布范围的预测更加准确。

Description

一种广布种分布范围的预测方法

技术领域

本发明涉及物种分布研究领域，具体涉及一种广布种分布范围的预测方法。

背景技术

广布种，顾名思义，是指能广泛分布于世界各大洋或淡水各区域中的生物种，广布种分布的研究，是生物多样性研究和保护的一个重要方面。广布种的遗传多样性既关系到如何采取科学有效的方法和措施来保护物种，又有助于研究物种的进化历史和适应能力。若能基于样本得到要研究物种在全国范围内的分布范围，就可以对物种进行有效地有针对性的保护，从而极大地减少人力和财力的输出，因此，建立模型对广布种的分布范围进行预测是领域内的一大研究热点。

目前，在对广布种的研究中，应用较广的方法是最大熵模型法，最大熵模型基于热力学第二定律，按照该定律，一个非均衡的生命系统通过与环境的物质和能量交换以保持其存在，也就是说，一个实测存在的系统具有“耗散”的特征，耗散使系统的熵不断增加，直至该生命系统与环境的熵最大，而使熵达到最大的状态，也是系统与环境之间的关系达到平衡的状态。在物种潜在分布的相关研究中，可将物种与其生长环境视为一个系统，通过计算系统具有最大熵时的状态参数确定物种和环境之间的稳定关系，并以此估计物种分布。基于该原则，最大熵模型在已知样本点和对应环境变量的基础上，通过拟合具有熵值最大的概率分布对物种的潜在分布做出估计。最大熵模型在预测物种相对多度、地理尺度上预测物种空间分布、研究宏生态学格局、推断物种间相互作用、解释食物网的度分布等领域的应用十分广泛。

但是，最大熵模型无法区分出观测数据中未观测到目标物种的记录是因为该地真的没有该物种分布，还是因为该地存在该物种，只是因为物种的可移动性等其他原因造成未观测到该物种，因此，采用最大熵模型对广布种的分布范围进行预测的效果并不理想。

发明内容

因此，本发明要解决的技术问题在于克服现有技术中的采用最大熵模型对广布种的分布范围进行预测的效果不理想的缺陷，从而提供一种广布种分布范围的预测方法。

为解决上述技术问题，本发明采用的技术方案为：

一种广布种分布范围的预测方法，包括以下步骤：

获取已知物种分布数量站点的物种数量y；

获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X；

根据二截尾零膨胀泊松模型计算参数γ、β；

所述二截尾零膨胀泊松模型为：

其中， 表示某地某物种真正不存在的概率，/>表示某地某物种存在但没有观测到的概率；

获取待预测物种分布概率站点的环境主成分变量W′、X′；

根据λ′＝exp(β^TX′)计算预测物种分布概率站点的模型参数/>λ′，并根据所述二截尾零膨胀泊松模型对物种分布概率进行预测。

进一步的，所述获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X的步骤包括：

获取已知物种分布数量站点的影响物种分布的环境协变量；

对所述环境协变量进行主成分降维，获得环境主成分变量。

进一步的，所述获取已知物种分布数量站点的影响物种分布的环境协变量的步骤包括：

获取已知物种分布数量站点的原始环境观测数据；

对所述原始环境观测数据进行标准化处理。

进一步的，所述对所述原始环境观测数据进行标准化处理的步骤包括：

对所述原始环境观测数据中的部分缺失项通过插值的方法进行填充；

对所述原始环境观测数据中缺失项较多的数据进行删除；

对所述原始环境观测数据的格式进行统一化整理；

对所述原始环境观测数据中明显异常点值进行处理。

进一步的，所述对所述原始环境观测数据中明显异常点值进行处理选自以下步骤A、B或C：

A、对所述明显异常点值进行取平均值操作；

B、采用合理取值对所述明显异常点值进行替代；

C、对所述明显异常点值进行删除。

进一步的，所述环境协变量包括站点经度、站点纬度、站点天气、站点温度、站点湿度、站点海拔以及站点降水量中的至少一种。

进一步的，所述根据二截尾零膨胀泊松模型计算参数γ、β的步骤包括：

当第i次观测的物种数量观测值y_i来自退化分布，记μ_i＝1；

当第i次观测的物种数量观测值y_i来自非退化分布，记μ_i＝0；

将μ＝(μ₁，μ₂，μ₃，……，μ_n)^T看作缺失数据，记为随机变量Y_m；

记物种观测数据y_i，X_i，W_i(i＝1,2，……,n)为随机变量Y₀；

记完全数据为Y_c＝(Y₀，Y_m)；

则样本的似然函数 其中，I为示性函数，满足下标的条件即为1，否则为0；

对所述似然函数两边同时取对数，得到似然方程；

解似然方程，得到γ、β。

本发明技术方案，具有如下优点：

1.本发明提供的广布种分布范围的预测方法，通过采用零膨胀模型对广布种的分布范围进行预测，零膨胀模型中的清楚地刻画了观测数据中的零样本是因为某地确实没有该物种存在，还是因为当次观测没有观测到该物种，从而可以从本质上对上述两种情况进行区分，有效地解决了广布种观测数据不平衡、零数据过多等问题，进而使得对广布种分布范围的预测更准确，另外，零膨胀模型可以处理不平衡数据集，其通过在模型的假设中对零样本的来源准确辨别，从而在实际上将不平衡数据集划分成了一个全是零样本的数据集和一个有零样本也有非零样本的平衡数据集，这样得到的结果可以很好的避免其他模型直接求解不平衡数据集的影响。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例1中的2017年斑腿树蛙在中国的分布范围预测图；

图2是中国鸟类野外手册中记载的强脚树莺的分布范围图；

图3是本发明可行性评价中对2017年强脚树莺在中国的分布范围预测图。

具体实施方式

提供下述实施例是为了更好地进一步理解本发明，并不局限于所述最佳实施方式，不对本发明的内容和保护范围构成限制，任何人在本发明的启示下或是将本发明与其他现有技术的特征进行组合而得出的任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。

实施例中未注明具体实验步骤或条件者，按照本领域内的文献所描述的常规实验步骤的操作或条件即可进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规试剂产品。

本发明涉及一种广布种分布范围的预测方法，包括以下步骤：

获取已知物种分布数量站点的物种数量y；

获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X；

根据二截尾零膨胀泊松模型计算参数γ、β；

所述二截尾零膨胀泊松模型为：

其中，λ＝exp(β^TX)，/>表示某地某物种真正不存在的概率，表示某地某物种存在但没有观测到的概率；

获取待预测物种分布概率站点的环境主成分变量W′、X′；

根据λ′＝exp(β^TX′)计算预测物种分布概率站点的模型参数/>λ′，并根据所述二截尾零膨胀泊松模型对物种分布概率进行预测。

需要说明的是，“T”表示的是向量转置。

优选的，所述获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X的步骤包括：

获取已知物种分布数量站点的影响物种分布的环境协变量；

对所述环境协变量进行主成分降维，获得环境主成分变量。

具体的，由于一条数据中包含的协变量过多，所以通过主成分分析的方法对协变量进行降维，可以达到筛选变量的目的。在选取主成分时，可以利用累积贡献率来进行主成分数量的选择，当主成分的累积贡献率大于0.9且主成分的数量少于等于3个时，此时选择使累积贡献率大于0.9的主成分变量，否则选取前三个主成分变量，由此得到影响物种分布的主成分变量。

优选的，所述获取已知物种分布数量站点的影响物种分布的环境协变量的步骤包括：

获取已知物种分布数量站点的原始环境观测数据；

对所述原始环境观测数据进行标准化处理。

优选的，所述对所述原始环境观测数据进行标准化处理的步骤包括：

对所述原始环境观测数据中的部分缺失项通过插值的方法进行填充；

对所述原始环境观测数据中缺失项较多的数据进行删除；

对所述原始环境观测数据的格式进行统一化整理；

对所述原始环境观测数据中明显异常点值进行处理。

优选的，所述对所述原始环境观测数据中明显异常点值进行处理选自以下步骤A、B或C：

A、对所述明显异常点值进行取平均值操作；

B、采用合理取值对所述明显异常点值进行替代；

C、对所述明显异常点值进行删除。

优选的，所述环境协变量包括站点经度、站点纬度、站点天气、站点温度、站点湿度、站点海拔以及站点降水量中的至少一种。

优选的，所述根据二截尾零膨胀泊松模型计算参数γ、β的步骤包括：

当第i次观测的物种数量观测值y_i来自退化分布，记μ_i＝1；

当第i次观测的物种数量观测值y_i来自非退化分布，记μ_i＝0；

将μ＝(μ₁，μ₂，μ₃，……，μ_n)^T看作缺失数据，记为随机变量Y_m；

记物种观测数据y_i，X_i，W_i(i＝1,2，……,n)为随机变量Y₀；

记完全数据为Y_c＝(Y₀，Y_m)；

则样本的似然函数 其中，I为示性函数，满足下标的条件即为1，否则为0；

对所述似然函数两边同时取对数，得到似然方程；

解似然方程，得到γ、β。

优选的，通过EM算法对样本的似然函数进行处理，得到γ、β。EM算法包括两步：

(1)E步，即求期望，根据已知物种分布数量站点的观测数据计算出的概率最大化Q函数：

其中，表示EM算法过程中第t步的参数估计值，且/>

(2)M步，即求最大值：

由于参数γ和β恰好分离在函数Q₁、Q₂中，所以执行M步，只要分别极大化Q₁、Q₂即可。通过计算，得到下面两个迭代公式：

其中，

重复上述EM算法的步骤直到收敛，得到γ、β。

实施例1

本实施例涉及2017年斑腿树蛙的分布范围预测方法，监测数据从生态环境部南京环境科学研究所获得，具体包括以下步骤：

S1、从生态环境部南京环境科学研究所的数据库中获取2017年度的若干个已知物种分布数量站点所观测记录的物种数量y；

S2、从生态环境部南京环境科学研究所的数据库中获取2017年度的已知物种分布数量站点的原始环境观测数据，包括站点经度、站点纬度、站点天气、站点温度、站点湿度、站点水温、站点pH、站点起始海拔、站点终点海拔、观测场海拔高度、气压传感器海拔高度、观测有效时段、站点气温、站点海拔、站点风速、站点风向以及站点降水量。

S3、对S2步骤中的原始环境观测数据进行标准化处理，即将斑腿树蛙的原始观测数据与国家气象数据网站上全国的气象站点的气象观测数据按照距离最近的原则进行拼接，得到拼接好的数据，并对数据的格式进行统一化整理，得到影响物种分布的环境协变量。其中，拼接是指两栖物种观测数据中缺少某一观测站点的气象数据，此时则以距离观测站点最近的气象站的气象数据作为观测站点的气象数据。

S4、对S3步骤的环境协变量进行主成分分析降维，得到环境主成分变量W、X。

S5、计算参数γ、β：

当第i次观测的物种数量观测值y_i来自退化分布，记μ_i＝1；

当第i次观测的物种数量观测值y_i来自非退化分布，记μ_i＝0；

将μ＝(μ₁，μ₂，μ₃，……，μ_n)^T看作缺失数据，记为随机变量Y_m；

记物种观测数据y_i，X_i，W_i(i＝1,2，……,n)为随机变量Y₀；

记完全数据为Y_c＝(Y₀，Y_m)；

则样本的似然函数 其中，I为示性函数，满足下标的条件即为1，否则为0；

对所述似然函数两边同时取对数，得到似然方程；

根据已知物种分布数量站点的观测数据计算出的概率最大化Q函数：

其中，表示EM算法过程中第t步的参数估计值，且/>

(2)M步，即求最大值：

由于参数γ和β恰好分离在函数Q₁、Q₂中，所以执行M步，只要分别极大化Q₁、Q₂即可。通过计算，得到下面两个迭代公式：

其中，

重复上述EM算法的步骤直到收敛，最终，得到

γ＝(1.18×10^-6，-1.38×10^-3，1.27×10^-2)

β＝(2.00×10^-4，5.77×10^-5，-9.07×10^-3)

S6、从生态环境部南京环境科学研究所的数据库中获取2017年度的待预测物种分布概率站点的原始环境观测数据，包括站点经度、站点纬度、站点天气、站点温度、站点湿度、站点水温、站点pH、站点起始海拔、站点终点海拔、观测场海拔高度、气压传感器海拔高度、观测有效时段、站点气温、站点海拔、站点风速、站点风向以及站点降水量。

S7、对S6步骤中的原始环境观测数据进行标准化处理，即将斑腿树蛙的原始观测数据与国家气象数据网站上全国的气象站点的气象观测数据按照距离最近的原则进行拼接，得到拼接好的数据，并对数据的格式进行统一化整理，得到影响物种分布的环境协变量。其中，拼接是指两栖物种观测数据中缺少某一观测站点的气象数据，此时则以距离观测站点最近的气象站的气象数据作为观测站点的气象数据。

S8、对S7步骤的环境协变量进行主成分分析降维，得到环境主成分变量W′、X′。

S9、根据λ′＝exp(β^TX′)计算预测物种分布概率站点的模型参数/>λ′；

S10、根据二截尾零膨胀泊松模型

计算各待预测物种分布概率站点的物种分布概率。

S10、根据各站点对应的经纬度和概率值进行散点图的绘制，将概率阈值定为0.5，若概率大于0.5，则预测此地斑腿泛树蛙存在，否则为不存在，得到斑腿泛树蛙在中国范围内的分布图，如图1所示。

需要说明的是，上述步骤S2-S9均是采用计算机直接计算的，即通过导入原始数据库来得到最后的各站点物种分布概率；本实施例的原始数据来源均来自于生态环境部南京环境科学研究所。

可行性评价

根据中国鸟类野外手册，强脚树莺的分布范围为喜马拉雅山脉至中国南方、东南亚以及大巽他群岛，如图2所示。目前，只有国内的数据被收集。因此若本方法的模型能将预测的范围与该分布图国内部分的实际分布区域相吻合，则说明采用该方法对广布种的分布进行预测是可行有效的。

按照实施例1的步骤，根据2017年全国气象站点的气象观测数据预测全国范围内强脚树莺的分布情况，得到对应站点(给定经纬度)的强脚树莺的分布概率，对2017年强脚树莺在中国范围内的分布范围的预测结果如图3所示。

根据图2与图3的对比，强脚树莺分布分布范围的预测范围与实际分布范围吻合度很高，说明本发明提供的预测方法对物种分布范围的预测是可行的。

另外，需要说明的一点是，《中国鸟类野外手册》为2000年出版，距今已有近20年的时间，在此期间，强脚树莺的分布范围已经向北扩展，这点可以从2014年-2017年观测到的数据得到验证。而零膨胀模型预测的结果也与之相同，因此，更进一步说明零膨胀模型在广布种的分布预测是可行的。

显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (7)

1.一种广布种分布范围的预测方法，其特征在于，包括以下步骤：

获取已知物种分布数量站点的物种数量y；

获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X；

根据二截尾零膨胀泊松模型计算参数γ、β；

所述二截尾零膨胀泊松模型为：

其中，λ＝exp(β^TX)，/>表示某地某物种真正不存在的概率，表示某地某物种存在但没有观测到的概率；

获取待预测物种分布概率站点的环境主成分变量W′、X′；

根据λ′＝exp(β^TX′)计算待预测物种分布概率站点的模型参数/>λ′，并根据所述二截尾零膨胀泊松模型对物种分布概率进行预测。

2.根据权利要求1所述的预测方法，其特征在于，所述获取已知物种分布数量站点的影响物种分布的环境主成分变量W、X的步骤包括：

获取已知物种分布数量站点的影响物种分布的环境协变量；

对所述环境协变量进行主成分降维，获得环境主成分变量。

3.根据权利要求2所述的预测方法，其特征在于，所述获取已知物种分布数量站点的影响物种分布的环境协变量的步骤包括：

获取已知物种分布数量站点的原始环境观测数据；

对所述原始环境观测数据进行标准化处理。

4.根据权利要求3所述的预测方法，其特征在于，所述对所述原始环境观测数据进行标准化处理的步骤包括：

对所述原始环境观测数据中的部分缺失项通过插值的方法进行填充；

对所述原始环境观测数据中缺失项较多的数据进行删除；

对所述原始环境观测数据的格式进行统一化整理；

对所述原始环境观测数据中明显异常点值进行处理。

5.根据权利要求4所述的预测方法，其特征在于，所述对所述原始环境观测数据中明显异常点值进行处理选自以下步骤A、B或C：

A、对所述明显异常点值进行取平均值操作；

B、采用合理取值对所述明显异常点值进行替代；

C、对所述明显异常点值进行删除。

6.根据权利要求2-5中任一项所述的预测方法，其特征在于，所述环境协变量包括站点经度、站点纬度、站点天气、站点温度、站点湿度、站点海拔以及站点降水量中的至少一种。

7.根据权利要求1-6中任一项所述的预测方法，其特征在于，所述根据二截尾零膨胀泊松模型计算参数γ、β的步骤包括：

当第i次观测的物种数量观测值y_i来自退化分布，记μ_i＝1；

当第i次观测的物种数量观测值y_i来自非退化分布，记μ_i＝0；

将μ＝(μ₁，μ₂，μ₃，……，μ_n)^T看作缺失数据，记为随机变量Y_m；

记物种观测数据y_i，X_i，W_i(i＝1,2，……,n)为随机变量Y₀；

记完全数据为Y_c＝(Y₀，Y_m)；

则样本的似然函数 其中，I为示性函数，满足下标的条件即为1，否则为0；

对所述似然函数两边同时取对数，得到似然方程；

解似然方程，得到γ、β。