森林轉(zhuǎn)基因林木抽樣策略的研究

趙毅強(qiáng), 蔡新宇, 秦朝秋, 高秉錚, 徐君怡*

1.中國農(nóng)業(yè)大學(xué)生物學(xué)院，北京 100193；2.大連海關(guān)技術(shù)中心，遼寧 大連 116001

林木轉(zhuǎn)基因技術(shù)的主要工作之一是對其性狀進(jìn)行改良。目前對林木基因性狀的改良主要集中在抗蟲、抗病原微生物、抗除草劑、抗非生物脅迫、降低木質(zhì)素合成以及增加纖維素合成等方面[1]。轉(zhuǎn)基因技術(shù)是林業(yè)前所未有的機(jī)遇，但是其潛在的負(fù)面影響同樣不容小覷。雖然轉(zhuǎn)基因林木的主要用途是在木材或生態(tài)景觀方面，不像轉(zhuǎn)基因農(nóng)作物一樣涉及人類食品健康與安全等相關(guān)問題，但是由于林木生產(chǎn)周期長、經(jīng)營管理粗放、風(fēng)媒傳粉等特點(diǎn)，可能存在轉(zhuǎn)基因林木對于林場以及周邊生態(tài)環(huán)境的環(huán)境釋放以及潛在的生態(tài)風(fēng)險(xiǎn)性[2]。一方面，如果轉(zhuǎn)基因林木的轉(zhuǎn)基因成分發(fā)生轉(zhuǎn)移并與天然林木的遺傳物質(zhì)進(jìn)行交流，則會(huì)使天然林木受到基因污染。另一方面，對轉(zhuǎn)基因林木抗蟲、抗病等優(yōu)勢性狀的自然選擇會(huì)增加相關(guān)抗性基因型在群體中的頻率，從而導(dǎo)致遺傳多態(tài)性減少，由此進(jìn)一步造成生態(tài)系統(tǒng)生物多樣性的破壞。此外，抗蟲抗病的轉(zhuǎn)基因林木可能會(huì)使害蟲、病菌進(jìn)行協(xié)同進(jìn)化變成“超級(jí)病蟲”[3]，進(jìn)一步打破生態(tài)平衡。因此，在這些問題被解決之前，建立一種在森林中轉(zhuǎn)基因林木監(jiān)測的抽樣策略十分必要。

由于森林面積遼闊，森林監(jiān)測中通常采用抽樣的辦法。傳統(tǒng)的抽樣方法，如五點(diǎn)取樣法、平行線取樣法、對角線取樣法、“Z”字形取樣法和棋盤式取樣法等，大多根據(jù)實(shí)驗(yàn)人員的經(jīng)驗(yàn)提出，更適合規(guī)則區(qū)域的田間抽樣[4]。這類方法無法滿足復(fù)雜條件下的抽樣均勻性，比如抽樣區(qū)域的不規(guī)則形狀、抽樣點(diǎn)的密度差異等。另一方面，在實(shí)際工作中抽樣點(diǎn)的確定需要配合森林的測繪信息。森林面積遼闊和抽樣區(qū)復(fù)雜的地質(zhì)特征導(dǎo)致測繪難度和人工成本增加[5]。無人機(jī)航測是傳統(tǒng)測量手段的補(bǔ)充，具有機(jī)動(dòng)靈活、作業(yè)成本低、適用范圍廣等優(yōu)點(diǎn)[6-7]。將無人機(jī)引入林業(yè)監(jiān)測可以大大節(jié)省人力成本，并對林區(qū)野外作業(yè)提供有力的支撐。

本研究期望開發(fā)一套森林轉(zhuǎn)基因林木抽樣方案，從而獲取真實(shí)抽樣區(qū)地形、制定抽樣策略到獲得抽樣點(diǎn)。在使用無人機(jī)獲得森林區(qū)域影像的基礎(chǔ)上，我們嘗試把聚類和空間分析的思想引入到森林轉(zhuǎn)基因林木抽樣方案中來，期望開發(fā)適用于不規(guī)則抽樣區(qū)域內(nèi)隨機(jī)抽樣的通用方法。

1 實(shí)驗(yàn)方法

1.1 抽樣樣本量計(jì)算

樣本量指研究包含的樣本例數(shù)。由于大多數(shù)研究都是通過樣本統(tǒng)計(jì)量推斷總體參數(shù)，樣本量過小，會(huì)導(dǎo)致由于抽樣誤差的影響不能真實(shí)地反映總體規(guī)律；樣本量過大，會(huì)增加不必要的資源開銷。合理的樣本量對于做出有效、可靠的科學(xué)結(jié)論至關(guān)重要。

假設(shè)轉(zhuǎn)基因的出現(xiàn)率在0.5%的水平屬于正常可控，達(dá)到2%的水平則提示存在擴(kuò)散風(fēng)險(xiǎn)，到達(dá)20%提示風(fēng)險(xiǎn)失控，實(shí)驗(yàn)人員可以通過預(yù)抽樣或者其他途徑獲得當(dāng)前可能的轉(zhuǎn)基因出現(xiàn)率作為備擇檢驗(yàn)率，通過二項(xiàng)式分布的正態(tài)近似來估計(jì)正式檢驗(yàn)所需的樣本量，從而判斷當(dāng)前的風(fēng)險(xiǎn)等級(jí)。在本研究中，我們通過SAS中的power過程進(jìn)行所需樣本量的計(jì)算。樣本量取決于以下4個(gè)條件：①假設(shè)檢驗(yàn)的一類錯(cuò)誤率α；②假設(shè)檢驗(yàn)的二類錯(cuò)誤率β；③零假設(shè)的率；④備擇假設(shè)的率。1-β即為通常所說的檢驗(yàn)效能，定義為備擇假設(shè)為真的情況下拒絕零假設(shè)的概率。使用下面的SAS程序給出相關(guān)參數(shù)不同取值的情況下，所需樣本量的梯度表。

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為設(shè)置的零假設(shè)，proportion為設(shè)置的備擇假設(shè)，sides=u表示備擇假設(shè)高于零假設(shè)的單尾檢驗(yàn)，ntotal=表示要估計(jì)樣本量，alpha為一類錯(cuò)誤率，power為檢驗(yàn)效能，即1-β。實(shí)驗(yàn)人員也可通過SAS的power過程獲得現(xiàn)有抽樣量情況下的檢驗(yàn)效能，從而判斷是否需要增加抽樣量。

1.2 無人機(jī)林區(qū)平面圖獲取和預(yù)處理

本研究使用無人機(jī)拍照獲取森林區(qū)域的影像。由于本研究并不需要進(jìn)行高精度定位，根據(jù)經(jīng)驗(yàn)10米級(jí)定位精度即可滿足要求，因此采用機(jī)載GPS數(shù)據(jù)進(jìn)行定位。用于定位和拍照無人機(jī)為DJ Mavic 2專業(yè)版，掛載高分辨率相機(jī)，使用飛控軟件Pix4Dcapture編程自動(dòng)飛行和拍照。將獲取的照片導(dǎo)入Agisoft PhotoScan Professional V1.43軟件，使用Align Photos功能對齊照片，并人工檢查對齊效果。軟件生成點(diǎn)云數(shù)據(jù)，并輸出正射影像用于后期分析。

獲得影像后對影像的測繪參數(shù)進(jìn)行設(shè)置。以圖片左上角設(shè)置為笛卡爾坐標(biāo)的原點(diǎn)，以圖片的寬度方向設(shè)為笛卡爾坐標(biāo)的X軸，以圖片高度方向設(shè)置為笛卡爾坐標(biāo)的Y軸，以一個(gè)像素為基本單位構(gòu)建笛卡爾坐標(biāo)系。如果數(shù)據(jù)本身來自測繪，則無須進(jìn)行這一步處理。我們以坐標(biāo)點(diǎn)的形式標(biāo)注出森林區(qū)域多邊形的頂點(diǎn)，順時(shí)針或逆時(shí)針方向依次排列構(gòu)成多邊形外邊框，并記錄于文件中。同樣的方法，標(biāo)注出森林區(qū)域內(nèi)部無法抽樣的孔洞多邊形的頂點(diǎn)，順時(shí)針或逆時(shí)針方向依次排列構(gòu)成孔洞多邊形的外邊框，并記錄于文件中。

注：公路左右兩側(cè)的森林區(qū)域獨(dú)立區(qū)分。黃色線條為標(biāo)記出的森林區(qū)域邊緣，紅色線條為標(biāo)記出的森林內(nèi)部無法抽樣的區(qū)域(僅右側(cè)圖中有)。

1.3 抽樣點(diǎn)算法

傳統(tǒng)的抽樣方法大多設(shè)計(jì)在規(guī)則多邊形，比如正方形、矩形等抽樣區(qū)域中。為了提高抽樣方法在實(shí)際抽樣過程中的可用性，抽樣方法應(yīng)適用于各種不規(guī)則的抽樣區(qū)域，包括不規(guī)則多邊形，以及帶有不可抽樣孔洞的不規(guī)則多邊形。為此，我們開發(fā)了不同的算法來實(shí)現(xiàn)這個(gè)目標(biāo)。對于抽樣點(diǎn)，我們規(guī)定其要滿足如下條件：①不要在多邊形的邊上；②抽樣點(diǎn)在多邊形內(nèi)部足夠均勻；③避開多邊形內(nèi)部無法抽樣的區(qū)域。

K-means是一種廣泛使用的聚類算法，其算法的基本思想是隨機(jī)選取K個(gè)空間坐標(biāo)點(diǎn)作為初始中心，并將空間中的每個(gè)離散點(diǎn)按度量距離分配到對應(yīng)的簇中。重新計(jì)算每個(gè)簇的質(zhì)心，作為新的中心并再次將空間中的散點(diǎn)按度量距離分配到對應(yīng)的簇中。重復(fù)上述過程，直到每個(gè)簇類中心不再發(fā)生變化而達(dá)到收斂。本研究中，度量距離使用歐式距離，抽樣點(diǎn)的數(shù)目通過SAS的power過程查表獲得，初始中心的位置在多邊形內(nèi)部隨機(jī)產(chǎn)生。為了適配不規(guī)則多邊形，在讀入多邊形頂點(diǎn)坐標(biāo)后，遍歷多邊形的每一個(gè)頂點(diǎn)，找到多邊形的包圍盒。為了簡化后續(xù)的抽樣點(diǎn)迭代算法，我們使用一組小正方形平鋪包圍盒。如果正方形的中心點(diǎn)在多邊形內(nèi)部，則認(rèn)為這個(gè)多邊形包含這個(gè)正方形。我們使用下面的算法來判斷正方形的中心點(diǎn)A是否在多邊形內(nèi)部：①隨機(jī)任一方向C，做一條AC射線；②遍歷整個(gè)多邊形頂點(diǎn)，判斷是否有頂點(diǎn)在AC射線上；③若有多邊形頂點(diǎn)在AC上則重復(fù)1、2過程，直到?jīng)]有多邊形頂點(diǎn)在AC射線上；④計(jì)算多邊形的每一條邊線是否與射線相交，并計(jì)量相交的次數(shù)；⑤如果相交的次數(shù)是偶數(shù)個(gè)則表示點(diǎn)在多邊形外部，若為奇數(shù)個(gè)則表示點(diǎn)在多邊形內(nèi)部。找出所有被多邊形包含的小正方形后，基于這組正方形的中心點(diǎn)集合進(jìn)行迭代操作。小正方形的數(shù)目和抽樣點(diǎn)的數(shù)目相關(guān)，設(shè)定其數(shù)目為抽樣點(diǎn)數(shù)目的平方乘以系數(shù)10，下限為1萬，上限為1 000萬。

本研究設(shè)計(jì)的第二個(gè)算法力求抽樣點(diǎn)分布更加均勻并考慮到內(nèi)部存在無法抽樣的孔洞的情況。馮洛諾伊圖(Voronoi Diagram)，由一組連接兩鄰點(diǎn)直線的垂直平分線組成的連續(xù)多邊形組成。馮洛諾伊圖是對空間平面的一種剖分，多邊形內(nèi)的任何位置離該多邊形樣點(diǎn)的距離最近，離相鄰多邊形內(nèi)樣點(diǎn)的距離遠(yuǎn)，且每個(gè)多邊形內(nèi)含且僅包含一個(gè)樣點(diǎn)。算法同樣隨機(jī)選取K個(gè)空間坐標(biāo)點(diǎn)作為初始樣點(diǎn)，計(jì)算K個(gè)點(diǎn)的馮洛諾伊圖，圖中每個(gè)劃分的區(qū)域稱為馮洛諾伊細(xì)胞。然后將每個(gè)馮洛諾伊細(xì)胞的樣點(diǎn)移動(dòng)到質(zhì)心，重新計(jì)算新的K個(gè)點(diǎn)的馮洛諾伊圖，并移動(dòng)各新細(xì)胞的樣點(diǎn)到新的質(zhì)心。重復(fù)上述過程，直到每個(gè)馮洛諾伊細(xì)胞質(zhì)心不再發(fā)生變化而達(dá)到收斂。最后各馮洛諾伊細(xì)胞的面積接近均等，其質(zhì)心在馮洛諾伊細(xì)胞的中心，稱為centroidal voronoi tessellation。CVT方法和K-means在算法上具有一定的相似性，在沒有內(nèi)部孔洞的凸多邊形中，通過馮洛諾伊圖算法得到的抽樣點(diǎn)更連續(xù)，更均勻。但是，由于存在孔洞的多邊形中，連接兩鄰點(diǎn)直線的垂直平分線可能被截?cái)?，?dǎo)致無法有效地生成馮洛諾伊細(xì)胞，我們受到CVT方法的啟發(fā)開發(fā)了一種近似的算法。首先在多邊形內(nèi)部密集產(chǎn)生均勻分布的離散點(diǎn)，并隨機(jī)劃分為K個(gè)包含盡可能相等數(shù)目離散點(diǎn)的區(qū)域。由于離散點(diǎn)分布均勻，我們使用每個(gè)區(qū)域包含的離散點(diǎn)數(shù)目表示該區(qū)域的面積(總的離散點(diǎn)數(shù)/K)。我們在每個(gè)區(qū)域內(nèi)選擇一個(gè)離散點(diǎn)，該點(diǎn)至區(qū)域內(nèi)其他離散點(diǎn)的可達(dá)最短距離之和最短，這個(gè)離散點(diǎn)作為這個(gè)區(qū)域的質(zhì)心。由于孔洞的存在，可達(dá)最短距離可能不再是直線，我們通過廣度優(yōu)先遍歷得到近似最短的離散點(diǎn)。如果區(qū)域A中一個(gè)離散點(diǎn)a與被相鄰區(qū)域B中的一個(gè)相鄰離散點(diǎn)a交換，可降低A與B區(qū)域的可達(dá)最短距離之和，則交換點(diǎn)ab。交換之后各區(qū)域的面積保持穩(wěn)定不變，但是形狀有調(diào)整，重新計(jì)算兩個(gè)區(qū)域質(zhì)心。重復(fù)上述步驟，直到不存在可交換離散點(diǎn)ab。此時(shí)各區(qū)域的質(zhì)心即為最終的取樣點(diǎn)。

對于上述兩種算法，我們編制了程序“一個(gè)不規(guī)則空間內(nèi)均勻采樣的工具軟件V1.0”。程序讀入記錄于文件中的多邊形頂點(diǎn)坐標(biāo)，和需要抽樣的數(shù)目。程序輸出每一個(gè)抽樣點(diǎn)的坐標(biāo)，以及每個(gè)抽樣點(diǎn)對應(yīng)的監(jiān)測區(qū)域。

表1 軟件的輸入輸出格式

2 結(jié)果與分析

2.1 假設(shè)轉(zhuǎn)基因出現(xiàn)率輸出power表格

通過SAS程序獲得不同零假設(shè)率、備擇假設(shè)率、檢驗(yàn)水平(α)和檢驗(yàn)效能(1-β)梯度情況下，所需的樣本量表。查表可知：如果預(yù)抽樣估計(jì)轉(zhuǎn)基因出現(xiàn)率為1.5%，在0.01的檢驗(yàn)水平和90%的檢驗(yàn)效能的情況下，需要不少于1 024個(gè)樣本來進(jìn)行統(tǒng)計(jì)，并做出轉(zhuǎn)基因出現(xiàn)率顯著高于正常可控水平的結(jié)論。同樣，如果預(yù)抽樣估計(jì)轉(zhuǎn)基因出現(xiàn)率為7%，在0.05的檢驗(yàn)水平和80%的檢驗(yàn)效能的情況下，需要不少于80個(gè)樣本來進(jìn)行統(tǒng)計(jì)，并做出轉(zhuǎn)基因出現(xiàn)率存在顯著風(fēng)險(xiǎn)的結(jié)論。如果預(yù)抽樣估計(jì)轉(zhuǎn)基因出現(xiàn)率為29%，在0.05的檢驗(yàn)水平和80%的檢驗(yàn)效能的情況下，需要不少于134個(gè)樣本來進(jìn)行統(tǒng)計(jì)，并做出轉(zhuǎn)基因林木擴(kuò)散已經(jīng)失控的結(jié)論(表2)。

表2 不同檢驗(yàn)參數(shù)下的抽樣數(shù)列表

2.2 K-means算法抽樣點(diǎn)迭代結(jié)果

根據(jù)表2的結(jié)果，以多邊形區(qū)域內(nèi)至少需要80個(gè)抽樣點(diǎn)為例，使用K-means算法對通過無人機(jī)獲得的森林左側(cè)區(qū)域進(jìn)行80個(gè)離散點(diǎn)的均勻抽樣。通過初始的隨機(jī)抽樣點(diǎn)130次迭代后，抽樣點(diǎn)的坐標(biāo)變化小于設(shè)定閾值迭代停止。在圖2中列出了第一次隨機(jī)抽樣點(diǎn)和最終抽樣點(diǎn)的位置，以及一次迭代中間過程的情況?？梢钥吹阶罱K的結(jié)果中抽樣點(diǎn)均勻分布在圖形中。程序記錄每一次迭代的過程，并輸出其中每一個(gè)抽樣點(diǎn)的坐標(biāo)，以及每個(gè)抽樣點(diǎn)對應(yīng)的監(jiān)測區(qū)域。

圖2 K-means算法中第一次隨機(jī)抽樣點(diǎn)、中途迭代抽樣點(diǎn)和最終抽樣點(diǎn)的位置

2.3 改進(jìn)的CVT算法抽樣點(diǎn)迭代結(jié)果

對于有孔洞的多邊形，我們使用改進(jìn)的CVT算法對森林的右側(cè)區(qū)域同樣進(jìn)行80個(gè)離散點(diǎn)的均勻抽樣。該算法計(jì)算量較大，通過310次迭代后迭代停止。可以看到最終的結(jié)果中抽樣點(diǎn)均勻分布在圖形中，且抽樣點(diǎn)被排除在多邊形中間的孔洞內(nèi)。圖3中列出了第一次隨機(jī)抽樣點(diǎn)和最終抽樣點(diǎn)的位置，以及一次迭代中間過程的情況。

圖3 改進(jìn)的CVT算法中第一次隨機(jī)抽樣點(diǎn)、中途迭代抽樣點(diǎn)和最終抽樣點(diǎn)的位置

3 討論

本研究開發(fā)了一套森林轉(zhuǎn)基因林木抽樣的方案，從獲取真實(shí)抽樣區(qū)地形、制定抽樣策略到獲得抽樣點(diǎn)。在使用無人機(jī)獲得森林區(qū)域影像的基礎(chǔ)上，嘗試把聚類和空間分析的思想引入到森林轉(zhuǎn)基因林木抽樣方案中來，開發(fā)適用于不規(guī)則抽樣區(qū)域內(nèi)隨機(jī)抽樣的通用方法。本研究提出的兩種方法均是基于成熟算法的改進(jìn)，提供了更好的抽樣準(zhǔn)確性和均勻性，也更適應(yīng)森林區(qū)域的實(shí)際情況，可作為森林轉(zhuǎn)基因林木抽樣策略的有力參考。

傳統(tǒng)的抽樣方法大多只適用于比較規(guī)則的抽樣區(qū)域，在實(shí)際運(yùn)用中效果會(huì)打折扣。機(jī)器學(xué)習(xí)聚類算法和空間分析算法已在多個(gè)領(lǐng)域中成功應(yīng)用，本研究首次把相關(guān)算法應(yīng)用于森林轉(zhuǎn)基因林木抽樣策略的研究。我們對相關(guān)算法進(jìn)行了改進(jìn)，使得算法能夠在不規(guī)則多邊形和帶有孔洞的多邊形中完成均勻抽樣，加強(qiáng)了抽樣方法在森林地區(qū)應(yīng)用的普適性。算法本身還可以做進(jìn)一步優(yōu)化，比如考慮到森林不同區(qū)域林木的密度和覆蓋度差異，對重點(diǎn)或者林木的密集區(qū)域進(jìn)行重點(diǎn)抽樣，這些可作為下一步工作的方向。

本研究建立的抽樣方法不僅可用于森林轉(zhuǎn)基因林木的監(jiān)測，也可以作為一種通用方法對森林的其他可疑情況進(jìn)行監(jiān)測，包括森林病蟲害監(jiān)測等。森林病蟲害種類繁多，可導(dǎo)致林木生長不良，產(chǎn)量質(zhì)量下降，甚至引起林木枯死。我國森林病蟲害大約有8 000余種，在全國范圍內(nèi)常發(fā)生的病蟲害有200多種[8]。我國每年發(fā)生病蟲害的森林面積接近1 000萬hm2，因森林病蟲害所造成的直接和間接損失接近1 000億元[9]。森林病蟲害的發(fā)生造成經(jīng)濟(jì)損失的同時(shí)，也造成了生態(tài)條件的惡化，并對人的居住環(huán)境產(chǎn)生影響。森林病蟲害以防為主，提前進(jìn)行監(jiān)測預(yù)報(bào)[10]，防止災(zāi)害發(fā)生發(fā)展，最大限度保持森林健康和生態(tài)平衡，避免產(chǎn)生不可估量的后果。需要指出的是，本研究建立的抽樣方法僅適用于區(qū)域內(nèi)的均勻抽樣，并不可直接應(yīng)用于具有特定傳播路徑事件的抽樣。對于這樣的事件，建議對傳播區(qū)域進(jìn)行提取后再進(jìn)行抽樣。

無人機(jī)作為一種新的獲取高分辨率影像的手段，具有成本低廉、操作簡單、靈活快捷等優(yōu)點(diǎn)[11]。無人機(jī)低空遙感在資源調(diào)查、基礎(chǔ)測繪中具有廣闊的前景。本研究中，無人機(jī)的應(yīng)用主要是通過正射影像圖獲得林區(qū)的二維平面圖。通過多傳感器、多角度觀測獲取三維點(diǎn)云數(shù)據(jù)進(jìn)行三維模型的重建，無人機(jī)可用于更為全面和復(fù)雜的監(jiān)測任務(wù)。另外，隨著計(jì)算機(jī)技術(shù)和影像設(shè)備的不斷發(fā)展，利用無人機(jī)遙感直接判定森林轉(zhuǎn)基因的類型，以及林木病理變化等成為可能，這對于森林動(dòng)態(tài)監(jiān)測，尤其是人工檢測不易進(jìn)行的區(qū)域具有重要的意義[12-13]。