李超 蘇耀文 涂文俊 張怡卓

(東北林業(yè)大學(xué)，哈爾濱，150040)

責(zé)任編輯:戴芳天。

根系在陸地生態(tài)系統(tǒng)能量轉(zhuǎn)換和物質(zhì)循環(huán)中扮演重要角色，不僅具有固定植株、吸收水分和礦物質(zhì)的功能，而且對于生態(tài)系統(tǒng)的碳循環(huán)也是至關(guān)重要的［1］。因此，研究植物根系的相關(guān)參數(shù)，如根徑大小、分布范圍及生物量，對生態(tài)系統(tǒng)循環(huán)具有重要意義。探地雷達(GPR)作為一種非侵入性探測方法，在根系探測領(lǐng)域正不斷改進與完善［2］。Raz－Yaseef N［3］等闡述了探地雷達方法在單株植物根徑大小和生物量的估測上具有實效性。郭立［4］等提出了粗根探地雷達信號正演模擬方法，分析了粗根探測精度的敏感因素，改進了探地雷達反演粗根生物量方法，建立了基于探地雷達的根系三維結(jié)構(gòu)模型。

由于探地雷達對根系進行高密度探測時，資源有限，而且采樣數(shù)據(jù)量大，無法實現(xiàn)實時壓縮，給傳輸信道造成巨大壓力［5］。屈樂樂［6］等提出了壓縮感知框架下的探地雷達技術(shù)，對可壓縮信號進行稀疏分解、觀測測量和重構(gòu)，有效減少了采集設(shè)備的壓力。正交匹配追蹤算法(OMP)是常用的壓縮感知重構(gòu)算法，但是，OMP 算法每一步分解都要在高維空間進行大量內(nèi)積計算，所需計算量巨大［7］。楊愚［8］等提出了基于粒子群算法(PSO)優(yōu)化的OMP重構(gòu)算法，通過粒子群快速搜索的特點，尋取最佳匹配原子，完成信號重構(gòu)。針對PSO 容易陷入局部最優(yōu)解，張建軍［9］等提出了基于混合粒子群的匹配追蹤算法，但是求出能逼近Hessian 矩陣的正定矩陣較為困難。趙知勁［10］等提出了量子粒子群的運動模式，該算法比基本粒子群算法速度快且精度高，但是由于粒子的中心位置參數(shù)設(shè)置為區(qū)間［0，1］的隨機數(shù)，造成粒子位置的隨機性，導(dǎo)致了算法的不穩(wěn)定。

模擬退火(SA)是全局尋優(yōu)算法，能以一定概率接受次優(yōu)解，能克服PSO 算法陷入局部最優(yōu)的缺陷［11］。于海平［12］等證明了PSO－SA 算法在多峰信號的應(yīng)用中，能有效增強PSO 算法的全局搜索能力和穩(wěn)定性。筆者提出PSO－SA－OMP 的優(yōu)化方法，針對探地雷達根系數(shù)據(jù)的快速重構(gòu)問題，通過PSO算法快速搜索，降低OMP 重構(gòu)算法運算時間;通過SA 算法，確定全局最優(yōu)解，提高重構(gòu)精度。

1 壓縮感知框架下的根系數(shù)據(jù)重構(gòu)

設(shè)m 為探地雷達時間采樣點數(shù)，n 為探地雷達采樣位置點，則可以表示成維度為m×n 的探地雷達數(shù)據(jù)，用N 維向量x(N=m×n)表示，x=［X1，1，X2，1，…，X1，2，X2，2，…，Xm，n］，其在某正交基組成的Ψ 域內(nèi)可表示為:

式中:Ψ=［Ψ1，Ψ2，Ψ3，…，ΨN］，N×1 的列向量s 是投影系數(shù)序列，如果s 中僅有K(K?N)個非零的大系數(shù)，則認(rèn)為x 在Ψ 域內(nèi)是K 稀疏的。

用M×N 的測量矩陣Φ(M?N)對x 進行觀測得到M 維的壓縮向量y:

Gabor 過完備原子庫在信號稀疏表示上具有靈活性，能更加準(zhǔn)確地表示稀疏信號，因此，采用Gabor 原子構(gòu)建過完備原子庫Θ。一個Gabor 原子由一個經(jīng)過調(diào)制的高斯窗函數(shù)組成，如式(3)所示。

式中:g(t)=e－πt2是高斯窗函數(shù);γ=(s，u，v，w)是時頻參數(shù)，一個原子由(s，u，v，w)決定;s 是伸縮因子;u 是位移原子;v 是原子的頻率;w 是原子相位。時頻參數(shù)可以按式(4)進行離散化。

式中:a、Δu 為系數(shù)參量;Δv、Δw 為相位參量;0＜j≤lbN，0≤p≤N2－j+1，0≤k＜2j+1，0≤i≤12。

當(dāng)Θ 滿足式(5)所示的約束等距性準(zhǔn)則(RIP)時，可通過求解式(6)問題進行恢復(fù)。

式中:0＜δk＜1;δk為等容常數(shù)。

顯然通過求解最小l0范數(shù)，從低維信號y 恢復(fù)出高維的稀疏信號s 或原始信號x 是NP－h(huán)ard問題。

OMP 從Θ 中選取構(gòu)成y 的列向量，重構(gòu)的關(guān)鍵步驟如下:

①在第p 次迭代中選出最佳原子ip，也就是使信號殘差rp－1與列g(shù)j相關(guān)性最大，即:

②將選中的原子構(gòu)成集合Hp，采用最小二乘:

③更新殘差:

在式(7)中，每一步分解都要在高維空間進行內(nèi)積計算，造成了其計算時間長?？梢赃x取OMP 的匹配函數(shù)|＜rp－1，gj＞|作為適應(yīng)度函數(shù)，將式(7)看作是求解最優(yōu)化問題，可以避免大量的內(nèi)積計算，明顯降低重構(gòu)算法的計算復(fù)雜度。

2 PSO-SA-OMP 重構(gòu)算法

PSO 算法是群體智能算法，源于鳥類捕食行為，算法中每個粒子都代表問題的一個潛在解，粒子在解空間中運動，不斷更新粒子的速度和位置，計算適應(yīng)度值，更新個體極值和群體極值，實現(xiàn)函數(shù)求解。SA 算法有加溫過程、等溫過程和冷卻過程，從某一較高初溫出發(fā)，伴隨溫度參數(shù)的不斷下降，結(jié)合概率突跳特性，在解空間中隨機尋找目標(biāo)函數(shù)的全局最優(yōu)解，能量變化就是目標(biāo)函數(shù)，最優(yōu)解就是能量最低態(tài)。

PSO 算法可以快速搜索匹配原子，但是在越來越接近最優(yōu)解時，粒子速度越來越小，導(dǎo)致粒子群出現(xiàn)“趨同和振蕩”現(xiàn)象，使算法不穩(wěn)定，易陷入局部最小值。而SA 算法能依據(jù)Metropolis 準(zhǔn)則，以一定概率接受次優(yōu)解，進行全局尋優(yōu)，結(jié)合兩者優(yōu)點，既減少了OMP 重構(gòu)過程的運算時間，又保證了重構(gòu)精度。

在d 維(M×N 維)搜索空間內(nèi)，Θ 匹配原子搜索優(yōu)化過程如下:

①初始化參數(shù):種群大小、慣性權(quán)重ω、學(xué)習(xí)因子c1和c2、初始溫度T0、迭代次數(shù)t 等。給每個粒子在Θ 中生成一隨機位置，每個位置對應(yīng)Θ 中的一個原子。

②第i 個粒子的速度vi=(vi，2vi，2…vi，d)、位置xi=(xi，1xi，2…xi，d)，粒子通過不斷計算匹配函數(shù)的適應(yīng)度值，選取每個粒子所經(jīng)歷過的具有最好適應(yīng)度值的位置作為個體最好位置，記為pi(pi1pi，2…pi，d);選取所有粒子所經(jīng)歷過的最好適應(yīng)度值的位置作為全局最好位置，記為pg(pg，1pg，2…pg，d)。速度和位置更新公式如式(10)、(11)所示。

粒子在解空間范圍的速度和位置限制在［－Vmdax，Vmdax］和［－Xmdax，Xmdax］。

③模擬退火算法在退溫搜索過程中，在解x 的區(qū)域中產(chǎn)生新的可行解x'時，計算f(x')和目標(biāo)函數(shù)f(x)的差值，依照概率min{1，exp(－Δf/Tt)}＞rand［0，1］接受解x'，其中rand［0，1］表示［0，1］的隨機數(shù)，因此能夠有效避免陷入局部極小值，Metropolis 準(zhǔn)則概率公式如式(12)所示。

式中:Tt表示第t 次迭代的溫度;Δf 表示計算目標(biāo)函數(shù)f(x')和目標(biāo)函數(shù)f(x)的差值。

根據(jù)以上闡述，設(shè)計PSO－SA－OMP 算法流程圖(見圖1)。

圖1 PSO－SA－OMP 重構(gòu)算法流程圖

3 仿真實驗及結(jié)果分析

3.1 重構(gòu)參數(shù)設(shè)計

探地雷達模型如圖2所示。XYZ 為探測空間，XOY 代表地表，XOZ 代表探測的地下剖面，發(fā)射天線(T)和接收天線(R)的距離為D，高度為l，根系上表面距地表的高度為h。當(dāng)探地雷達在某一位置探測時，此時采集到的回波數(shù)據(jù)為A－scan 數(shù)據(jù)，回波中的唯一變量為時間;當(dāng)探地雷達沿某一直線探測多次時，多道A－scan 數(shù)據(jù)組成B－scan 數(shù)據(jù)。

在XOZ 平面中，構(gòu)造一個40 cm×40 cm 的二維目標(biāo)空間，記作A，將其空間等分為L=1 600 個網(wǎng)格，記作A=［A1，…Ai，…AL］T，用b=［b1，…bi，…bL］T表示A 的加權(quán)系數(shù)向量，b 中的元素取布爾量，0 表示沒有目標(biāo)，1 表示有目標(biāo)。

系統(tǒng)采用收發(fā)分置天線，天線距離地表高度l=15 cm，收發(fā)天線間隔D=3 cm，天線移動步長1 cm，系統(tǒng)對接收信號采樣256 點，采樣間隔10 ps，采用壓縮感知的高斯隨機矩陣滿足均值為0，方差為1/，觀測維數(shù)為128 維。粒子群和模擬退火中主要參數(shù)選擇如下:種群大小為50，最大迭代次數(shù)為30，學(xué)習(xí)因子c1=c2=2.1，模擬退火初始溫度T0=－fitness(gbest)/log(0.2)，退火常數(shù)λ=0.2，選取50個原子重構(gòu)信號。

圖2 根系探地雷達模型圖

3.2 A-scan 數(shù)據(jù)重構(gòu)

按照參數(shù)設(shè)置，首先選取其中一組A－scan 數(shù)據(jù)進行實驗。當(dāng)發(fā)射波進入地面遇有目標(biāo)時，信號幅度會發(fā)生變化。壓縮感知信號重構(gòu)效果，可以通過式(13)所示的均方誤差(EMS)進行度量。

式中:N=256;f'(x)表示重構(gòu)信號;f(x)表示原始信號。

圖3是各類算法對A－scan 數(shù)據(jù)稀疏重構(gòu)的對比。由圖3a 得出，在信號長度50～256，原始數(shù)據(jù)受到雜波影響;而由圖3b、圖3c、圖3d 可知，采用壓縮感知后，僅保留了具有稀疏特性的目標(biāo)信號，很大程度上去除了雜波的干擾。

圖3 A－scan 數(shù)據(jù)重構(gòu)方法比較

表1為各類算法的MSE 和運算時間。在Ascan 數(shù)據(jù)重構(gòu)上，PSO－SA－OMP 算法的MSE 較小，但是PSO－SA－OMP 算法比傳統(tǒng)OMP 算法運算速度明顯提高，運算時間減少了10.471 s。

表1 A－scan 數(shù)據(jù)比較

3.3 B-scan 數(shù)據(jù)重構(gòu)

多個A－scan 數(shù)據(jù)組成B－scan 數(shù)據(jù)。定義重構(gòu)B－scan 數(shù)據(jù)圖像的峰值信噪比RPSN，如式(14)所示。

圖4 B－scan 數(shù)據(jù)重構(gòu)算法比較

重構(gòu)圖像表明，PSO－SA－OMP 比傳統(tǒng)OMP和PSO－OMP 重構(gòu)的B－scan 數(shù)據(jù)圖像更為精細，表2給出了PSNR 和時間的比較。PSO－SAOMP 重構(gòu)算法不僅在PSNR 方面比傳統(tǒng)OMP 算法提高了5.539 dB，而且運算時間減少了20.260 s;雖然PSO－OMP 算法運算速度更快，但是信噪比不高，本研究算法計算效率明顯提高。

表2 B－scan 數(shù)據(jù)比較

3.4 復(fù)雜探地雷達數(shù)據(jù)重構(gòu)

采用美國GSSI 實際采集的根系探地雷達數(shù)據(jù)進行測試，由上述實驗已經(jīng)得出，PSO－SA－OMP 重構(gòu)算法在理想的仿真數(shù)據(jù)下不僅效果好，而且對雜波具有抑制作用。對于復(fù)雜的實際根系探地雷達數(shù)據(jù)，實驗結(jié)果如圖5所示。

圖5a 和圖5d 為實際數(shù)據(jù)圖像，圖5b 和圖5e采用壓縮感知理論的PSO－SA－OMP 算法對實際根系探地雷達數(shù)據(jù)進行重構(gòu)。圖5c 和圖5f 為雜波抑制并進行圖像二值化分割后得到的目標(biāo)信息。分割圖像表明，雜波能得到很好抑制，目標(biāo)清晰。

4 結(jié)論

提出一種快速有效的根系探地雷達數(shù)據(jù)壓縮感知重構(gòu)算法，粒子群算法可以快速尋找OMP 過程中每一步分解的匹配原子，降低了運算復(fù)雜度，提高了重構(gòu)的速度。模擬退火算法避免了收斂到局部最優(yōu)，提高了信號重構(gòu)精度。實驗結(jié)果表明，PSO－SAOMP 算法能對探地雷達數(shù)據(jù)進行快速重構(gòu)，降低了OMP 算法的計算復(fù)雜度與提高了重構(gòu)精度，同時，優(yōu)化算法具有雜波抑制作用，對于實際采集的根系探地雷達數(shù)據(jù)，能得到更加清晰的目標(biāo)信息。

圖5 PSO－SA－OMP 對復(fù)雜探地雷達數(shù)據(jù)測試

［1］ Wu Y，Guo L，Cui X，et al.Ground-penetrating radar-based automatic reconstruction of three-dimensional coarse root system architecture［J］.Plant and Soil，2014，383(1/2):155－172.

［2］ Salvador J M，Jimenez V G，Lopez R，et al.Platform for research and education on ground penetrating radar［J］.Radar Sensor Technology XVII，2013，8714(4):813－831.

［3］ Raz-Yaseef N，Koteen L，Baldocchi D D.Coarse root distribution of a semi-arid oak savanna estimated with ground penetrating radar［J］.Journal of Geophysical Research:Biogeosciences，2013，118(1):135－147.

［4］ 郭立，范碧航，吳淵，等.探地雷達應(yīng)用于植物粗根探測的研究進展［J］.中國科技論文，2014(4):494－498.

［5］ 盧策吾，劉小軍，方廣有.基于感知壓縮的探地雷達數(shù)據(jù)壓縮采集［J］.電子學(xué)報，2011，39(9):2204－2206.

［6］ 屈樂樂，黃瓊，方廣有，等.基于壓縮感知的頻率步進探地雷達成像算法［J］.系統(tǒng)工程與電子技術(shù)，2010，32(2):295－297.

［7］ Liu J，Xu S，Gao X，et al.Novel imaging methods of stepped frequency radar based on compressed sensing［J］.Journal of Systems Engineering and Electronics，2012，23(1):47－56.

［8］ 楊愚.利用粒子群算法實現(xiàn)信號OMP 稀疏分解［J］.微計算機信息，2008，24(43):178－201.

［9］ 張建軍，王仲生，余匯.采用混合粒子群算法實現(xiàn)匹配追蹤算法［J］.振動與沖擊，2010，29(1):143－147.

［10］ 趙知勁，馬春暉.一種基于量子粒子群的二次匹配OMP 重構(gòu)算法［J］.計算機工程與應(yīng)用，2012，48(29):157－161.

［11］ Bahmani-Firouzi B.A new hybrid algorithm based on PSO，SA，and k-means for cluster analysis［J］.International Journal of Innovative Computing，Information and Control，2010，6(7):3177－3192.

［12］ 于海平，劉會超，吳志健.基于模擬退火的自適應(yīng)粒子群優(yōu)化算法的改進策略［J］.計算機應(yīng)用研究，2012，29(12):4448－4450.