李超 蘇耀文 涂文俊 張怡卓

(東北林業(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)具有重要意義。探地雷達(dá)(GPR)作為一種非侵入性探測方法，在根系探測領(lǐng)域正不斷改進(jìn)與完善［2］。Raz－Yaseef N［3］等闡述了探地雷達(dá)方法在單株植物根徑大小和生物量的估測上具有實(shí)效性。郭立［4］等提出了粗根探地雷達(dá)信號正演模擬方法，分析了粗根探測精度的敏感因素，改進(jìn)了探地雷達(dá)反演粗根生物量方法，建立了基于探地雷達(dá)的根系三維結(jié)構(gòu)模型。

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

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

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

設(shè)m 為探地雷達(dá)時(shí)間采樣點(diǎn)數(shù)，n 為探地雷達(dá)采樣位置點(diǎn)，則可以表示成維度為m×n 的探地雷達(dá)數(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)個(gè)非零的大系數(shù)，則認(rèn)為x 在Ψ 域內(nèi)是K 稀疏的。

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

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

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

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

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

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

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

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

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

②第i 個(gè)粒子的速度vi=(vi，2vi，2…vi，d)、位置xi=(xi，1xi，2…xi，d)，粒子通過不斷計(jì)算匹配函數(shù)的適應(yīng)度值，選取每個(gè)粒子所經(jīng)歷過的具有最好適應(yīng)度值的位置作為個(gè)體最好位置，記為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'時(shí)，計(jì)算f(x')和目標(biāo)函數(shù)f(x)的差值，依照概率min{1，exp(－Δf/Tt)}＞rand［0，1］接受解x'，其中rand［0，1］表示［0，1］的隨機(jī)數(shù)，因此能夠有效避免陷入局部極小值，Metropolis 準(zhǔn)則概率公式如式(12)所示。

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

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

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

3 仿真實(shí)驗(yàn)及結(jié)果分析

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

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

在XOZ 平面中，構(gòu)造一個(gè)40 cm×40 cm 的二維目標(biāo)空間，記作A，將其空間等分為L=1 600 個(gè)網(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，天線移動(dòng)步長1 cm，系統(tǒng)對接收信號采樣256 點(diǎn)，采樣間隔10 ps，采用壓縮感知的高斯隨機(jī)矩陣滿足均值為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è)原子重構(gòu)信號。

圖2 根系探地雷達(dá)模型圖

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

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

式中: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 和運(yùn)算時(shí)間。在Ascan 數(shù)據(jù)重構(gòu)上，PSO－SA－OMP 算法的MSE 較小，但是PSO－SA－OMP 算法比傳統(tǒng)OMP 算法運(yùn)算速度明顯提高，運(yùn)算時(shí)間減少了10.471 s。

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

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

多個(gè)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ù)圖像更為精細(xì)，表2給出了PSNR 和時(shí)間的比較。PSO－SAOMP 重構(gòu)算法不僅在PSNR 方面比傳統(tǒng)OMP 算法提高了5.539 dB，而且運(yùn)算時(shí)間減少了20.260 s;雖然PSO－OMP 算法運(yùn)算速度更快，但是信噪比不高，本研究算法計(jì)算效率明顯提高。

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

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

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

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

4 結(jié)論

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

圖5 PSO－SA－OMP 對復(fù)雜探地雷達(dá)數(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］ 郭立，范碧航，吳淵，等.探地雷達(dá)應(yīng)用于植物粗根探測的研究進(jìn)展［J］.中國科技論文，2014(4):494－498.

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

［6］ 屈樂樂，黃瓊，方廣有，等.基于壓縮感知的頻率步進(jìn)探地雷達(dá)成像算法［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］ 楊愚.利用粒子群算法實(shí)現(xiàn)信號OMP 稀疏分解［J］.微計(jì)算機(jī)信息，2008，24(43):178－201.

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

［10］ 趙知?jiǎng)?，馬春暉.一種基于量子粒子群的二次匹配OMP 重構(gòu)算法［J］.計(jì)算機(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ìn)策略［J］.計(jì)算機(jī)應(yīng)用研究，2012，29(12):4448－4450.