森林轉基因林木抽樣策略的研究

趙毅強, 蔡新宇, 秦朝秋, 高秉錚, 徐君怡*

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]。