森林转基因林木抽样策略的研究

赵毅强, 蔡新宇, 秦朝秋, 高秉铮, 徐君怡*

1.中国农业大学生物学院，北京 100193；2.大连海关技术中心，辽宁 大连 116001

林木转基因技术的主要工作之一是对其性状进行改良。目前对林木基因性状的改良主要集中在抗虫、抗病原微生物、抗除草剂、抗非生物胁迫、降低木质素合成以及增加纤维素合成等方面[1]。转基因技术是林业前所未有的机遇，但是其潜在的负面影响同样不容小觑。虽然转基因林木的主要用途是在木材或生态景观方面，不像转基因农作物一样涉及人类食品健康与安全等相关问题，但是由于林木生产周期长、经营管理粗放、风媒传粉等特点，可能存在转基因林木对于林场以及周边生态环境的环境释放以及潜在的生态风险性[2]。一方面，如果转基因林木的转基因成分发生转移并与天然林木的遗传物质进行交流，则会使天然林木受到基因污染。另一方面，对转基因林木抗虫、抗病等优势性状的自然选择会增加相关抗性基因型在群体中的频率，从而导致遗传多态性减少，由此进一步造成生态系统生物多样性的破坏。此外，抗虫抗病的转基因林木可能会使害虫、病菌进行协同进化变成“超级病虫”[3]，进一步打破生态平衡。因此，在这些问题被解决之前，建立一种在森林中转基因林木监测的抽样策略十分必要。

由于森林面积辽阔，森林监测中通常采用抽样的办法。传统的抽样方法，如五点取样法、平行线取样法、对角线取样法、“Z”字形取样法和棋盘式取样法等，大多根据实验人员的经验提出，更适合规则区域的田间抽样[4]。这类方法无法满足复杂条件下的抽样均匀性，比如抽样区域的不规则形状、抽样点的密度差异等。另一方面，在实际工作中抽样点的确定需要配合森林的测绘信息。森林面积辽阔和抽样区复杂的地质特征导致测绘难度和人工成本增加[5]。无人机航测是传统测量手段的补充，具有机动灵活、作业成本低、适用范围广等优点[6-7]。将无人机引入林业监测可以大大节省人力成本，并对林区野外作业提供有力的支撑。

本研究期望开发一套森林转基因林木抽样方案，从而获取真实抽样区地形、制定抽样策略到获得抽样点。在使用无人机获得森林区域影像的基础上，我们尝试把聚类和空间分析的思想引入到森林转基因林木抽样方案中来，期望开发适用于不规则抽样区域内随机抽样的通用方法。

1 实验方法

1.1 抽样样本量计算

样本量指研究包含的样本例数。由于大多数研究都是通过样本统计量推断总体参数，样本量过小，会导致由于抽样误差的影响不能真实地反映总体规律；样本量过大，会增加不必要的资源开销。合理的样本量对于做出有效、可靠的科学结论至关重要。

假设转基因的出现率在0.5%的水平属于正常可控，达到2%的水平则提示存在扩散风险，到达20%提示风险失控，实验人员可以通过预抽样或者其他途径获得当前可能的转基因出现率作为备择检验率，通过二项式分布的正态近似来估计正式检验所需的样本量，从而判断当前的风险等级。在本研究中，我们通过SAS中的power过程进行所需样本量的计算。样本量取决于以下4个条件：①假设检验的一类错误率α；②假设检验的二类错误率β；③零假设的率；④备择假设的率。1-β即为通常所说的检验效能，定义为备择假设为真的情况下拒绝零假设的概率。使用下面的SAS程序给出相关参数不同取值的情况下，所需样本量的梯度表。

proc power;

onesamplefreq test=z method=normal

nullproportion = 0.005 0.02 0.2

proportion = 0.005 to 0.4 by 0.005

sides = u

ntotal =.

alpha = 0.01 0.05

power = 0.8 0.9;

ods output output=outtable(where=(NTotal>0)keep=NullProportion Alpha Proportion Power NTotal);

run;

proc export data= outtable

outfile= "C:sample_size.xls"

dbms=xls replace;

run;

上面的代码中，nullproportion为设置的零假设，proportion为设置的备择假设，sides=u表示备择假设高于零假设的单尾检验，ntotal=表示要估计样本量，alpha为一类错误率，power为检验效能，即1-β。实验人员也可通过SAS的power过程获得现有抽样量情况下的检验效能，从而判断是否需要增加抽样量。

1.2 无人机林区平面图获取和预处理

本研究使用无人机拍照获取森林区域的影像。由于本研究并不需要进行高精度定位，根据经验10米级定位精度即可满足要求，因此采用机载GPS数据进行定位。用于定位和拍照无人机为DJ Mavic 2专业版，挂载高分辨率相机，使用飞控软件Pix4Dcapture编程自动飞行和拍照。将获取的照片导入Agisoft PhotoScan Professional V1.43软件，使用Align Photos功能对齐照片，并人工检查对齐效果。软件生成点云数据，并输出正射影像用于后期分析。

获得影像后对影像的测绘参数进行设置。以图片左上角设置为笛卡尔坐标的原点，以图片的宽度方向设为笛卡尔坐标的X轴，以图片高度方向设置为笛卡尔坐标的Y轴，以一个像素为基本单位构建笛卡尔坐标系。如果数据本身来自测绘，则无须进行这一步处理。我们以坐标点的形式标注出森林区域多边形的顶点，顺时针或逆时针方向依次排列构成多边形外边框，并记录于文件中。同样的方法，标注出森林区域内部无法抽样的孔洞多边形的顶点，顺时针或逆时针方向依次排列构成孔洞多边形的外边框，并记录于文件中。

注：公路左右两侧的森林区域独立区分。黄色线条为标记出的森林区域边缘，红色线条为标记出的森林内部无法抽样的区域(仅右侧图中有)。

1.3 抽样点算法

传统的抽样方法大多设计在规则多边形，比如正方形、矩形等抽样区域中。为了提高抽样方法在实际抽样过程中的可用性，抽样方法应适用于各种不规则的抽样区域，包括不规则多边形，以及带有不可抽样孔洞的不规则多边形。为此，我们开发了不同的算法来实现这个目标。对于抽样点，我们规定其要满足如下条件：①不要在多边形的边上；②抽样点在多边形内部足够均匀；③避开多边形内部无法抽样的区域。

K-means是一种广泛使用的聚类算法，其算法的基本思想是随机选取K个空间坐标点作为初始中心，并将空间中的每个离散点按度量距离分配到对应的簇中。重新计算每个簇的质心，作为新的中心并再次将空间中的散点按度量距离分配到对应的簇中。重复上述过程，直到每个簇类中心不再发生变化而达到收敛。本研究中，度量距离使用欧式距离，抽样点的数目通过SAS的power过程查表获得，初始中心的位置在多边形内部随机产生。为了适配不规则多边形，在读入多边形顶点坐标后，遍历多边形的每一个顶点，找到多边形的包围盒。为了简化后续的抽样点迭代算法，我们使用一组小正方形平铺包围盒。如果正方形的中心点在多边形内部，则认为这个多边形包含这个正方形。我们使用下面的算法来判断正方形的中心点A是否在多边形内部：①随机任一方向C，做一条AC射线；②遍历整个多边形顶点，判断是否有顶点在AC射线上；③若有多边形顶点在AC上则重复1、2过程，直到没有多边形顶点在AC射线上；④计算多边形的每一条边线是否与射线相交，并计量相交的次数；⑤如果相交的次数是偶数个则表示点在多边形外部，若为奇数个则表示点在多边形内部。找出所有被多边形包含的小正方形后，基于这组正方形的中心点集合进行迭代操作。小正方形的数目和抽样点的数目相关，设定其数目为抽样点数目的平方乘以系数10，下限为1万，上限为1 000万。

本研究设计的第二个算法力求抽样点分布更加均匀并考虑到内部存在无法抽样的孔洞的情况。冯洛诺伊图(Voronoi Diagram)，由一组连接两邻点直线的垂直平分线组成的连续多边形组成。冯洛诺伊图是对空间平面的一种剖分，多边形内的任何位置离该多边形样点的距离最近，离相邻多边形内样点的距离远，且每个多边形内含且仅包含一个样点。算法同样随机选取K个空间坐标点作为初始样点，计算K个点的冯洛诺伊图，图中每个划分的区域称为冯洛诺伊细胞。然后将每个冯洛诺伊细胞的样点移动到质心，重新计算新的K个点的冯洛诺伊图，并移动各新细胞的样点到新的质心。重复上述过程，直到每个冯洛诺伊细胞质心不再发生变化而达到收敛。最后各冯洛诺伊细胞的面积接近均等，其质心在冯洛诺伊细胞的中心，称为centroidal voronoi tessellation。CVT方法和K-means在算法上具有一定的相似性，在没有内部孔洞的凸多边形中，通过冯洛诺伊图算法得到的抽样点更连续，更均匀。但是，由于存在孔洞的多边形中，连接两邻点直线的垂直平分线可能被截断，导致无法有效地生成冯洛诺伊细胞，我们受到CVT方法的启发开发了一种近似的算法。首先在多边形内部密集产生均匀分布的离散点，并随机划分为K个包含尽可能相等数目离散点的区域。由于离散点分布均匀，我们使用每个区域包含的离散点数目表示该区域的面积(总的离散点数/K)。我们在每个区域内选择一个离散点，该点至区域内其他离散点的可达最短距离之和最短，这个离散点作为这个区域的质心。由于孔洞的存在，可达最短距离可能不再是直线，我们通过广度优先遍历得到近似最短的离散点。如果区域A中一个离散点a与被相邻区域B中的一个相邻离散点a交换，可降低A与B区域的可达最短距离之和，则交换点ab。交换之后各区域的面积保持稳定不变，但是形状有调整，重新计算两个区域质心。重复上述步骤，直到不存在可交换离散点ab。此时各区域的质心即为最终的取样点。

对于上述两种算法，我们编制了程序“一个不规则空间内均匀采样的工具软件V1.0”。程序读入记录于文件中的多边形顶点坐标，和需要抽样的数目。程序输出每一个抽样点的坐标，以及每个抽样点对应的监测区域。

表1 软件的输入输出格式

2 结果与分析

2.1 假设转基因出现率输出power表格

通过SAS程序获得不同零假设率、备择假设率、检验水平(α)和检验效能(1-β)梯度情况下，所需的样本量表。查表可知：如果预抽样估计转基因出现率为1.5%，在0.01的检验水平和90%的检验效能的情况下，需要不少于1 024个样本来进行统计，并做出转基因出现率显著高于正常可控水平的结论。同样，如果预抽样估计转基因出现率为7%，在0.05的检验水平和80%的检验效能的情况下，需要不少于80个样本来进行统计，并做出转基因出现率存在显著风险的结论。如果预抽样估计转基因出现率为29%，在0.05的检验水平和80%的检验效能的情况下，需要不少于134个样本来进行统计，并做出转基因林木扩散已经失控的结论(表2)。

表2 不同检验参数下的抽样数列表

2.2 K-means算法抽样点迭代结果

根据表2的结果，以多边形区域内至少需要80个抽样点为例，使用K-means算法对通过无人机获得的森林左侧区域进行80个离散点的均匀抽样。通过初始的随机抽样点130次迭代后，抽样点的坐标变化小于设定阈值迭代停止。在图2中列出了第一次随机抽样点和最终抽样点的位置，以及一次迭代中间过程的情况。可以看到最终的结果中抽样点均匀分布在图形中。程序记录每一次迭代的过程，并输出其中每一个抽样点的坐标，以及每个抽样点对应的监测区域。

图2 K-means算法中第一次随机抽样点、中途迭代抽样点和最终抽样点的位置

2.3 改进的CVT算法抽样点迭代结果

对于有孔洞的多边形，我们使用改进的CVT算法对森林的右侧区域同样进行80个离散点的均匀抽样。该算法计算量较大，通过310次迭代后迭代停止。可以看到最终的结果中抽样点均匀分布在图形中，且抽样点被排除在多边形中间的孔洞内。图3中列出了第一次随机抽样点和最终抽样点的位置，以及一次迭代中间过程的情况。

图3 改进的CVT算法中第一次随机抽样点、中途迭代抽样点和最终抽样点的位置

3 讨论

本研究开发了一套森林转基因林木抽样的方案，从获取真实抽样区地形、制定抽样策略到获得抽样点。在使用无人机获得森林区域影像的基础上，尝试把聚类和空间分析的思想引入到森林转基因林木抽样方案中来，开发适用于不规则抽样区域内随机抽样的通用方法。本研究提出的两种方法均是基于成熟算法的改进，提供了更好的抽样准确性和均匀性，也更适应森林区域的实际情况，可作为森林转基因林木抽样策略的有力参考。

传统的抽样方法大多只适用于比较规则的抽样区域，在实际运用中效果会打折扣。机器学习聚类算法和空间分析算法已在多个领域中成功应用，本研究首次把相关算法应用于森林转基因林木抽样策略的研究。我们对相关算法进行了改进，使得算法能够在不规则多边形和带有孔洞的多边形中完成均匀抽样，加强了抽样方法在森林地区应用的普适性。算法本身还可以做进一步优化，比如考虑到森林不同区域林木的密度和覆盖度差异，对重点或者林木的密集区域进行重点抽样，这些可作为下一步工作的方向。

本研究建立的抽样方法不仅可用于森林转基因林木的监测，也可以作为一种通用方法对森林的其他可疑情况进行监测，包括森林病虫害监测等。森林病虫害种类繁多，可导致林木生长不良，产量质量下降，甚至引起林木枯死。我国森林病虫害大约有8 000余种，在全国范围内常发生的病虫害有200多种[8]。我国每年发生病虫害的森林面积接近1 000万hm2，因森林病虫害所造成的直接和间接损失接近1 000亿元[9]。森林病虫害的发生造成经济损失的同时，也造成了生态条件的恶化，并对人的居住环境产生影响。森林病虫害以防为主，提前进行监测预报[10]，防止灾害发生发展，最大限度保持森林健康和生态平衡，避免产生不可估量的后果。需要指出的是，本研究建立的抽样方法仅适用于区域内的均匀抽样，并不可直接应用于具有特定传播路径事件的抽样。对于这样的事件，建议对传播区域进行提取后再进行抽样。

无人机作为一种新的获取高分辨率影像的手段，具有成本低廉、操作简单、灵活快捷等优点[11]。无人机低空遥感在资源调查、基础测绘中具有广阔的前景。本研究中，无人机的应用主要是通过正射影像图获得林区的二维平面图。通过多传感器、多角度观测获取三维点云数据进行三维模型的重建，无人机可用于更为全面和复杂的监测任务。另外，随着计算机技术和影像设备的不断发展，利用无人机遥感直接判定森林转基因的类型，以及林木病理变化等成为可能，这对于森林动态监测，尤其是人工检测不易进行的区域具有重要的意义[12-13]。