頤和園古樹內(nèi)部異常的雷達波檢測

李偉林, 文 劍, 肖中亮, 張 京, 韓紅巖

(1.北京林業(yè)大學工學院，北京 100083；2.北京市頤和園管理處，北京 100091)

頤和園古樹內(nèi)部異常的雷達波檢測

李偉林1, 文 劍1, 肖中亮1, 張 京2, 韓紅巖2

(1.北京林業(yè)大學工學院，北京 100083；2.北京市頤和園管理處，北京 100091)

采用雷達波樹木無損探測技術(shù)與樹木雷達分析軟件，對頤和園內(nèi)6棵古柳的內(nèi)部異常情況進行雷達波探測成像試驗.結(jié)果表明古柳內(nèi)部均存在不同程度的異常，其中3棵古柳異常程度較高.雷達波成像解析結(jié)果準確且與實際情況相符程度高，表明樹木雷達系統(tǒng)可以對樹木內(nèi)部異常分布情況進行高效無損的檢測、成像和解析.

樹木雷達; 無損檢測; 介電常數(shù); 樹木異常; 頤和園

樹木因自然缺陷或生物降解而影響生長質(zhì)量，主要異常表現(xiàn)為樹木內(nèi)部的開裂、結(jié)疤、腐朽和中空等損壞[1-2].健康樹木的內(nèi)部結(jié)構(gòu)可以分成心材、邊材和樹皮3個主要部分，不同結(jié)構(gòu)中的成分關(guān)系很復雜.有些異常是從樹木外部可觀測到的，而有些則隱藏在樹干內(nèi)部，需要通過有損或無損的檢測方法進行識別.樹干無損探測方法有很多種，包括目測法、應力波法[3-4]、電阻率法[5]、樹木針法、超聲波法[6-7]、X射線(CT)掃描法[8]、雷達波法等.

在對名木古樹進行探測時盡量減少對其造成損壞.雷達波技術(shù)能夠?qū)δ繕宋矬w進行完全無損的探測，與應力波、超聲波、X射線掃描等方法相比，其精度高、無污染、探測安全，因此越來越多地被用于木質(zhì)體和樹木的探測[9].雷達波是使用電磁波來探測介質(zhì)電特性的技術(shù)，當電磁波在介質(zhì)中傳播時，會穿過電特性有差異的平面，部分電磁波被反射，并被接收天線接收.Godio et al[10]和Sambuelli et al[11]將探地雷達層析成像技術(shù)與超聲波、電阻率技術(shù)進行了比較.Rodríguez et al[12]使用探地雷達分析了樹木的非均質(zhì)性，并指出它的介電常數(shù)變化依賴于紋路方向(木材纖維縱向排列).探地雷達也被用于木質(zhì)橫梁內(nèi)部結(jié)構(gòu)的檢測[13].Fu et al[14]將探地雷達用于活立木成像.Butnor et al[15]研究發(fā)現(xiàn)探地雷達對松柏類、針葉類樹木的樹皮表面孔洞、內(nèi)部空氣孔洞的檢測效果很好.雷達波在土木工程[16]、文化遺產(chǎn)保護、運輸工程[17]、土壤濕度測繪[18]等領(lǐng)域也得到了廣泛應用.

由于樹木內(nèi)部結(jié)構(gòu)不同，電磁特性會受含水率影響，而且樹木內(nèi)部的密度、溫度和進行檢測時所用的頻率不同，以及樹木截面輪廓具有的特異性，使得樹木雷達數(shù)據(jù)的分析變得很復雜[19].在對頤和園寅輝城關(guān)的一級紅牌古樹油松進行無損探測時發(fā)現(xiàn)，目前常用的Tree Win分析軟件雖然可以反映樹木內(nèi)部結(jié)構(gòu)情況，但其圖像中的樹木輪廓近似圓形，不能如實復現(xiàn)實際樹木截面輪廓，而且不能識別內(nèi)部異常形狀及大小.為了快速有效地對頤和園界河橋旁6棵古柳的樹干進行探測，本研究利用樹木雷達分析軟件，使用探地雷達采集樹木截面內(nèi)部數(shù)據(jù)，分析雷達波圖像中的可視化特征，在無損樹木的情況下對古柳內(nèi)部結(jié)構(gòu)情況進行測定；并且結(jié)合樹木截面輪廓信息呈現(xiàn)樹木內(nèi)部可視化圖像，準確預測樹木內(nèi)部結(jié)構(gòu)情況.

1 材料與方法

1.1 研究對象

研究區(qū)位于北京市海淀區(qū)頤和園管理區(qū)界湖橋旁，地理坐標：北緯39°59′N，東經(jīng)116°15′E，海拔52 m.該地區(qū)屬于溫帶季風性氣候.頤和園占地300.8 hm2，擁有古樹名木 1 600余株，界湖橋位于西堤北端，始建于乾隆年間，其沿岸古柳距今已有百余年歷史.

北京市頤和園管理處曾對園內(nèi)界湖橋旁6棵古柳進行過樹干修復作業(yè).以界河橋為起點，自北向南將古柳進行編號，從YHY201611-L0001到Y(jié)HY201611-L0006；平均樹齡在200 a左右；最大周長(YHY201611-L0003)達到5.21 m，最小周長(YHY201611-L0002)為2.90 m.這6棵古柳的樹齡是頤和園古樹中較大的，且生長狀況較差.

1.2 雷達波探測

雷達波以高頻帶短脈沖的形式由發(fā)射天線T發(fā)射，不同層位和異常反射的回波被接收天線接收.由于在傳遞過程中，雷達波會受到外部干擾，出現(xiàn)低頻漂移、水平道間干擾、雜波信號干擾等現(xiàn)象，需要對回波數(shù)據(jù)進行初始化與振幅增益、小波變換、零點校正、噪聲去除等預處理[20]，提高信噪比和探測精度.然后采用匹配濾波器算法對預處理后的信號進行層位識別，得到不同層位的回波時延，進而確定雷達波在介質(zhì)層內(nèi)的傳播時間.接著采用逐層反演的方法估計各個層位的介電常數(shù)，獲取當前介質(zhì)中的雷達波傳輸波速，由回波時延和傳播時間確定層位的厚度.

(1)

(2)

式中，Vp表示電磁波在介質(zhì)中傳播速度，c表示電磁波在真空中的傳播速度，ε′是相對介電常數(shù)，Z表示層位厚度，t表示電磁波雙程的時間.

1.3 樹木層位分析

在確定樹木內(nèi)部各層位厚度時，采用匹配濾波器算法計算回波時延，分析樹木內(nèi)部結(jié)構(gòu)的不同層位.匹配濾波是信號檢測中比較常用的一種方法，需要根據(jù)待識別目標的空間特性，設(shè)計與其幅值相匹配的濾波器參數(shù)，以便通過有效地增強信號能量而抑制干擾能量，從而使輸出的信噪比達到最大.采用匹配濾波器檢測算法[21]分析樹木內(nèi)部的層位結(jié)構(gòu)，匹配濾波器的脈沖響應h(t)可表示為：

h(t)=x(T-t)

(3)

式中，T表示信號x(t)持續(xù)的時間.

依據(jù)N-P準則，匹配濾波器中的閾值St取決于容許的最大虛警率Pf，即

(4)

式中，E是已知信號x(t)信號的能量，erfc-1(Pf)是虛警率Pf的逆互補誤差函數(shù).

采用式(3)計算匹配濾波器的脈沖響應h(t)，對檢測到的回波信號yr(t)進行匹配濾波并獲得濾波信號yMF(t)，找到y(tǒng)MF(t)絕對值的最大值ymax及其對應的時間延遲tmax；將yMF最大值與設(shè)定閾值St進行比較，根據(jù)被檢測信號與原始信號的極性，匹配濾波器會輸出最大值或最小值；通過最大時間延遲可推算出每個層位的時間延遲.

1.4 全內(nèi)反射現(xiàn)象

樹干內(nèi)部結(jié)構(gòu)可以近似看成圓柱體，要使用探地雷達觀測這類圓柱形物體，必須了解從不同介質(zhì)之間界面散射出的電磁波的內(nèi)部反射.

n1gsin(?i)=n2gsin(?r)

(5)

式中，n1表示出射介質(zhì)折射率，n2表示入射介質(zhì)折射率，?i表示入射角，?r表示折射或反射角.

部分折射和部分反射的電磁波強度會逐漸變?nèi)?，而全?nèi)反射幾乎保留了電磁波所有信號.實現(xiàn)全內(nèi)反射有2個必要條件：(1)電磁波從光密介質(zhì)向光疏介質(zhì)傳播；(2)入射角?i大于臨界角?c.

當折射角或反射角?r等于90°時，入射角?i即為對應的臨界角?c.這種情況下sin(?r)等于1，臨界角可由下式表示：

(6)

電磁波在2種不同介質(zhì)間的傳播方式可分為3種：(1)入射角小于臨界角，電磁波發(fā)生部分折射；(2)電磁波以臨界角入射時?r為90°，即沿臨界面?zhèn)鞑ィ?3)入射角大于臨界角，電磁波不能穿過臨界面，以等同于入射角的角度發(fā)生全反射.

在圓柱形物體中電磁波的傳播路徑可以有各種形式，而要記錄一種反射波，就必須能夠?qū)⑵浣邮盏?探地雷達圖像可以反映外表面、內(nèi)部目標和全內(nèi)反射的不同形式.Pinel et al[22]也描述了粗糙表面反射和折射的特征.

1.5 角相關(guān)函數(shù)方法

圖1 圓周探測原理圖Fig.1 The principle diagram of circumference detection

樹木雷達的收發(fā)天線共用，探測時天線沿截面的圓周測線運動，雷達波收發(fā)點的坐標為rsn=(xsn,ysn),n=1,2,K,N，樹木截面中的點為r0=(x0,y0)，如圖1所示.

在rsn處接收到的信號為sn(t)，sn(t)通過聚焦濾波器hn(t,r0)使信號聚焦在r0[23].聚焦濾波器傅里葉變換后的聚焦函數(shù)如式(7)所示.對各個掃描點的數(shù)據(jù)進行相位校正，當r0處有目標時，使各個掃描點的數(shù)據(jù)聚焦后在r0(x0,y0)相位相同.成像處理后聚焦濾波器的輸出如下式所示：

(7)

r0n=|rsn-r0|

(8)

(9)

為了提高圖像的信噪比，可以利用角相關(guān)技術(shù)對聚焦濾波器輸出的信號進行相關(guān)疊加，角相關(guān)函數(shù)表示為：

(10)

式中，arg(rsi,rsj)表示半徑軸rsi與rsj之間的夾角，△Φ稱為相關(guān)角，可取不同的夾角.

通過聚焦濾波器和角相關(guān)函數(shù)的方法，可以將樹木雷達采集到的B-scan數(shù)據(jù)轉(zhuǎn)換成極坐標形式；根據(jù)樹木截面輪廓的形狀進行對比，從而生成與實際內(nèi)部結(jié)構(gòu)接近的圖像，達到復現(xiàn)的效果.根據(jù)層位識別的結(jié)果標記正常與異常區(qū)域，生成易于解釋的可視化圖像，有助于古樹維護人員準確掌握樹木內(nèi)部結(jié)構(gòu)的變化，并及時采取相關(guān)措施.

1.6 檢測方法

樹木內(nèi)部探測系統(tǒng)使用的探地雷達是美國TRU(Tree Radar Unit)雷達系統(tǒng)，包括2個獨立設(shè)備，即頻率為900 MHz的探地雷達天線和數(shù)據(jù)管理器.分析成像軟件是采用Matlab編寫的基于雷達波回波數(shù)據(jù)的樹木雷達分析軟件.

采集數(shù)據(jù)時，使用TRU系統(tǒng)天線對每棵樹5個不同高度(0.6、0.9、1.2、1.5、1.8 m)部位(圖2)進行探測，先用細繩標記5個高度平面和正北起始軸線，并且測量各截面樹干周長.然后從正北方向開始順時針掃描，保持勻速，緊貼樹皮，掃描一圈到正北方向結(jié)束，并記錄數(shù)據(jù)編號.同一高度重復3次.

樹木雷達分析軟件的流程如圖3所示，數(shù)據(jù)采集后進行數(shù)據(jù)讀取及振幅增益、零點校正、噪聲去除等預處理，獲取B-scan數(shù)據(jù).采用匹配濾波器法獲取異常層面的時延信息，進行定位；采用逐層反演的方法估計各個層位的介電常數(shù)；并精確估計層位的厚度，識別結(jié)構(gòu)中的異常；結(jié)合全反射、角相關(guān)函數(shù)等方法以及外輪廓數(shù)據(jù)[24]，得到木質(zhì)體內(nèi)部異常的圖像.

圖2 樹木雷達掃描示意圖Fig.2 Sketch map of tree radar

2 結(jié)果與分析

2.1 樹干內(nèi)部的異常層位

根據(jù)雷達波異常層位分析方法可知，由于介電常數(shù)的差異性，雷達波會在木質(zhì)體與異常處產(chǎn)生明顯的回波信號，并出現(xiàn)在B-scan圖中.使用樹木雷達系統(tǒng)對頤和園界湖橋旁古柳進行探測的現(xiàn)場如圖4所示.圖5為編號YHY201611-L0001(高度0.6 m)古柳的探地雷達B-scan圖，通過灰度級來描述接收信號的幅度，橫向軸表示探測長度，深度軸表示探測深度.從圖5可以觀察到幾次多重反射，在反射面呈現(xiàn)出相同的形狀(上部綠線)；其他反射與樹木內(nèi)部存在的結(jié)構(gòu)特征相對應.由于不同介質(zhì)間的介電常數(shù)存在差異，雷達波在穿越介質(zhì)界面時會產(chǎn)生相應的反射特征，這些特征描述了存在異常與樹皮表面的相對位置，如圖5中紅色線所示.為了更準確直觀地呈現(xiàn)內(nèi)部異常的分布，需要使用角相關(guān)函數(shù)的方法將雷達層位掃描信息轉(zhuǎn)換為內(nèi)部絕對位置，并進行分布成像.

圖4 頤和園檢測現(xiàn)場Fig.4 On-site detection with tree radar

2.2 基于輪廓的角相關(guān)函數(shù)成像掃描

對于樹干這類圓柱型被測物體，需要通過圖像變換方法獲得更直觀的成像效果.TRU系統(tǒng)的分析軟件Tree Win可使用極坐標的形式生成雷達波數(shù)據(jù)分析圖像.如在頤和園的前期試驗中，對一級紅牌古樹油松進行了無損探測，并使用Tree Win軟件進行數(shù)據(jù)分析，結(jié)果如圖6所示.從圖6可以看出，在大多數(shù)樹木截面近似原型的情況下，可較為直觀地呈現(xiàn)樹木內(nèi)部缺陷的分布，也可以反映樹木內(nèi)部結(jié)構(gòu)缺陷的分布情況.采用自行研制的樹木雷達分析軟件以及角相關(guān)函數(shù)方法，對雷達波數(shù)據(jù)進行變換，并與輪廓形狀數(shù)據(jù)進行對比，生成分析圖像(圖7).從圖7可以看出，與Tree Win軟件結(jié)果相比，角相關(guān)函數(shù)結(jié)合輪廓的圖像變換方法,是利用輪廓掃描數(shù)據(jù)以及雷達波深度層位分析信息將雷達波深度的相對坐標轉(zhuǎn)換為實際二維掃描截面的絕對坐標，可更準確、直觀地復現(xiàn)樹干內(nèi)部截面圖像.

a.實物；b.B-scan;c.Tree Win分析.
圖6頤和園古樹分析圖
Fig.6 Abnormality analysis of trees by B-scan

a.實物；b.Tree Win分析;c.樹木雷達軟件分析.

圖7TreeWin與樹木雷達軟件分析對比圖
Fig.7 Comparison of Tree Win and tree radar

2.3 頤和園現(xiàn)場樹干內(nèi)部的異常成像

分別使用樹木雷達系統(tǒng)和樹木雷達分析軟件對頤和園界湖橋旁古柳進行探測，并將探測分析結(jié)果與Tree Win軟件分析結(jié)果進行對比.圖8a～8f為編號YHY201611-L0001到Y(jié)HY201611-L0006古柳內(nèi)部分析圖.圖8各個小圖中，左圖的藍色線表示5個不同截面高度，紅色線表示正北方向基準線；中間圖為樹木雷達軟件分析結(jié)果；右圖為經(jīng)Tree Win軟件分析出來的結(jié)果.通過對比可以觀察到，樹木雷達系統(tǒng)得到的探測結(jié)果比Tree Win系統(tǒng)更接近實際情況，而且可以更直觀地看出利用采集到的雷達波數(shù)據(jù)對樹木內(nèi)部腐朽等異常的定位識別和成像效果更好.

將6棵古柳樹干探測結(jié)果進行量化，內(nèi)部異常區(qū)域占比(表1)顯示頤和園界河橋旁的這6棵古柳全都有不同程度的異常.YHY201611-L0004古柳的異常占比最高，達到76.92%；YHY201611-L0002古柳的異常占比最低，為54.37%.所有測量對象異常占比平均值高達67.76%，其中YHY201611-L0004、YHY201611-L0005、YHY201611-L0001存在異常程度極高，易發(fā)生傾倒等.

表1 古柳樹干內(nèi)部異常區(qū)域占比Table 1 The proportion of abnormal area within tree trunks

a.YHY201611-L0001;b.YHY201611-L0002;c.YHY201611-L0003;d.YHY201611-L0004;e.YHY201611-L0005;d.YHY201611-L0005;e.YHY201611-L0006.

3 小結(jié)

本文采用基于雷達波的樹木內(nèi)部層位結(jié)構(gòu)及異常探測成像系統(tǒng)，對頤和園界河橋旁6棵古柳進行探測.針對樹木內(nèi)部層位多、雷達回波在不同介質(zhì)層的傳遞特性差異大的特點，采用匹配濾波器算法計算回波時延，確定層位的厚度，結(jié)果表明樹木雷達分析軟件可以提高樹木內(nèi)部異常區(qū)域的分層探測精度.頤和園探測結(jié)果表明6棵古柳全部都有不同程度的異常.本文采用的環(huán)繞切向測量方法，可得到樹木截面信息，并生成便于直接觀測的識別圖像，有助于名木古樹的快速無損檢測.

[1] 王岳,王海燕,楊小娟,等.基于土壤因子的近天然落葉松云冷杉林枯損模型研究[J].福建農(nóng)林大學學報(自然科學版),2015,44(4):379-383.

[2] 岳小泉,王立海,葛曉雯.風荷載作用下樹木力學特性研究進展[J].森林工程,2015,31(6):33-36.

[3] 廖春暉,張厚江,王喜平,等.明十三陵裕陵修繕工程中的木構(gòu)件無損檢測[J].北京林業(yè)大學學報,2014,36(1):132-137.

[4] 房友盼,劉英,徐兆軍,等.運用靜態(tài)法與應力波法測定寬度對硬材彈性模量的影響[J].森林工程,2015,31(4):60-63.

[5] 岳小泉,王立海,劉澤旭,等.不同含水率下利用電阻與應力波斷層成像定量檢測木材腐朽[J].福建農(nóng)林大學學報(自然科學版),2016,45(5):593-598.

[6] 曾祥,周曉軍,楊辰龍,等.基于經(jīng)驗模態(tài)分解和S變換的缺陷超聲回波檢測方法[J].農(nóng)業(yè)機械學報,2016,47(11):414-420.

[7] 梅萍,夏兆鵬,張耀麗,等.檀香紫檀和染料紫檀快速微損鑒別方法[J].福建農(nóng)林大學學報(自然科學版),2017,46(2):154-158.

[8] 葛浙東,侯曉鵬,魯守銀,等.基于反投影坐標快速算法的木材CT檢測系統(tǒng)研究[J].農(nóng)業(yè)機械學報,2016,47(3):327-341.

[9] GAO L, LIU J, XUE H. Nondestructive detection of standing trees and radar wave detection[C]∥IEEE. International Conference on Computing, Control and Industrial Engineering. Wnhan: IEEE Computer Society, 2010:304-307.

[10] GODIO A, SAMBUELLI L, SOCCO L. Application and comparison of three tomographic techniques for detection of decay in trees[J]. Journal of Arboriculture, 2003,29(2):66-78.

[11] SAMBUELLI L, SOCCO L V, GODIO A. Ultrasonic, electric and radar measurements for living trees assessment[J]. Bollettino Di Geofisica Teorica Ed Applicata, 2003,44(3):253-279.

[12] RODRIGUEZ-ABAD I, MARTINEZ-SALA R, GARCIA-GARCIA F, et al. Non-destructive methodologies for the evaluation of moisture content in sawn timber structures: ground-penetrating radar and ultrasound techniques[J]. Near Surface Geophysics, 2010,8(6):475-482.

[13] PEREZ-GRACIA V, SANTOS-ASSUNCAO S, CASELLES O, et al. Study of wood beams in buildings with ground penetrating radar[C]∥IEEE. International Conference on Ground Penetrating Radar. Brussels: IEEE, 2014:31-35.

[14] FU L, LIU S, LIU L. Internal structure characterization of living tree trunk cross-section using GPR: numerical examples and field data analysis[C]∥IEEE. 15th International Conference on Ground Penetrating Radar. Brussels:IEEE, 2014:155-160.

[15] BUTNOR J R, PRUYN M L, SHAW D C, et al. Detecting defects in conifers with ground penetrating radar: applications and challenges[J]. Forest Pathology, 2009,39(5):309-322.

[16] NUNEZ N X, SOLLA M, NOVO A, et al. Three-dimensional ground-penetrating radar methodologies for the characterization and volumetric reconstruction of underground tunneling[J]. Construction & Building Materials, 2014,71:551-560.

[17] BOLD R D, O′CONNOR G, MORRISSEY J, et al. Benchmarking large scale GPR experiments on railway ballast[J]. Construction & Building Materials, 2015,92:31-42.

[18] MINET J, BOGAERT P, VANCLOOSTER M, et al. Validation of ground penetrating radar full-waveform inversion for field scale soil moisture mapping[J]. Journal of Hydrology, 2012,425(5):112-123.

[19] TOMPPO L, TIITTA M, LAAKSO T, et al. Dielectric spectroscopy of Scots pine[J]. Wood Science & Technology, 2009,43(7-8):653-667.

[20] XIAO Z L, WEN J, GAO L, et al. Method of tree radar signal processing based on Curvelet transform[J]. Revista Tecnica de la Facultad de Ingenieria Universidad del Zulia, 2016,39(7):243-250.

[21] LAHOUR S. Development of Data Analysis Algorithms for Interpretation of Ground Penetrating Radar Data[D]. Blacksburrg: Virginia Polytechnic Institute and State University, 2003.

[22] PINEL N, BOURLIER C, PINEL N, et al. Electromagnetic Wave Scattering from Random Rough Surfaces: Basics[M]. USA: Wiley, 2013:1-39.

[23] 張安學,蔣延生,汪文秉.圓周探地雷達測量和成像方法的研究[J].電子學報,2002,30(6):853-856.

[24] 林耀海,景林,王長纓,等.基于圓弧的原木輪廓的識別與驗證[J].福建農(nóng)林大學學報(自然科學版),2016,45(6):649-654.

DetectionofinternalanomalyfortreesintheSummerPalacebasedonradar

LI Weilin1, WEN Jian1, XIAO Zhongliang1, ZHANG Jing2, HAN Hongyan2

(1.School of Technology, Beijing Forestry University, Beijing 100083, China; 2.The Summer Palace Management Office, Beijing 100091, China)

The internal abnormal situation of 6 ancient willows in the Summer Palace was detected based on tree radar system and radar analysis software. The results showed that all 6 willows had anomalies, 3 of which have high proportions of anomalies, indicating that the tree radar system is applicable to predict the internal distribution of tree and contour shape non-invasively.

tree radar; nondestructive testing; permittivity; trunk abnormality; Summer Palace

2017-02-26

2017-05-15

國家自然科學基金資助項目(31600589)；中央高校基本科研業(yè)務費專項資金資助項目(2016ZCQ08).

李偉林(1993-)，男，碩士研究生.研究方向：現(xiàn)代傳感與檢測技術(shù).Email：liweilin@bjfu.edu.cn.通訊作者文劍(1981-)，男，講師,博士.研究方向：林業(yè)無損檢測技術(shù).Email：wenjian@bjfu.edu.cn.

S781; S792.12

A

1671-5470(2017)06-0665-07

10.13323/j.cnki.j.fafu(nat.sci.).2017.06.011

(責任編輯:葉濟蓉)