師 君張曉玲 韋順軍 向 高 楊建宇

(電子科技大學(xué)電子工程學(xué)院 成都 611731)

1 引言

陣列3維SAR[1－8]是指利用線型天線陣列代替?zhèn)鹘y(tǒng)合成孔徑雷達(dá)的單天線結(jié)構(gòu)，通過控制天線陣列的運(yùn)動(dòng)，合成2維等效陣列，并結(jié)合脈沖壓縮技術(shù)獲得3維空間分辨能力，以實(shí)現(xiàn)對(duì)觀測(cè)目標(biāo)3維成像的技術(shù)。與傳統(tǒng)SAR相比，陣列3維SAR具有下視3維成像能力，在地形測(cè)繪、災(zāi)害監(jiān)測(cè)等領(lǐng)域具有廣泛的應(yīng)用前景。

但是，由于陣列3維SAR需要在機(jī)翼上布設(shè)線陣天線以獲得該方向的分辨率，載機(jī)平臺(tái)尺寸限制使得其陣列方向分辨率遠(yuǎn)遠(yuǎn)低于距離向和航跡向。該缺陷嚴(yán)重制約了陣列3維SAR系統(tǒng)整體性能的提升。因此，需要針對(duì)陣列3維SAR數(shù)據(jù)特點(diǎn)，從成像/圖像處理的角度研究提高陣列向分辨率的方法。

由于3維SAR圖像顯著的稀疏特性[9]，稀疏重建技術(shù)，如L1正則化方法，OMP，基追蹤等技術(shù)目前已經(jīng)被廣泛用于層析SAR、圓周SAR、軌道3維SAR的成像處理，在提高圖像分辨率、消除旁瓣等方面取得了一定的效果。但在深入研究后發(fā)現(xiàn)，直接將稀疏重建技術(shù)用于陣列3維SAR圖像增強(qiáng)問題仍存在一些問題。

首先，稀疏特性并不能完全描述陣列3維SAR圖像所具有的典型特征。在地形測(cè)繪任務(wù)中，目標(biāo)散射點(diǎn) 3維空間中沿著某個(gè)特定的曲面(數(shù)字高程圖/DEM)分布，圖像增強(qiáng)應(yīng)直接針對(duì)DEM圖進(jìn)行處理。而在采用稀疏重建技術(shù)時(shí)，DEM 圖特有的“單值性”、“連續(xù)性”等特征被重述為“稀疏性”，實(shí)際上導(dǎo)致了約束條件的放松。單值性的缺失可能導(dǎo)致構(gòu)造的高程圖存在一個(gè)x,y位置對(duì)應(yīng)多個(gè)高度的不合理情況，最終導(dǎo)致DEM提取誤差。

其次，由于陣列3維SAR圖像為3維數(shù)據(jù)，采用稀疏重建技術(shù)可能需要構(gòu)造上億維的線性方程組，并進(jìn)行最優(yōu)化求解，大大增加了算法的實(shí)現(xiàn)難度和可行性。而降維處理策略則會(huì)導(dǎo)致額外的模型誤差，進(jìn)一步影響DEM增強(qiáng)效果。

針對(duì)上述問題，本文討論了基于變分模型的陣列3維SAR最優(yōu)DEM重建方法，該方法直接將DEM 圖作為最優(yōu)化目標(biāo)，通過尋找最優(yōu)化 DEM圖和對(duì)應(yīng)的散射系數(shù)，實(shí)現(xiàn)最小二乘意義下的最優(yōu)DEM重建。第2節(jié)簡(jiǎn)要回顧了陣列3維SAR的原理和成像處理方法，第3節(jié)對(duì)陣列3維SAR最優(yōu)DEM重建的變分模型進(jìn)行了分析討論。第4節(jié)討論了變分模型的求解方法，并提出了一種基于變分模型的陣列3維SAR最優(yōu)DEM重建方法。第5節(jié)對(duì)該方法的性能進(jìn)行了分析比較，驗(yàn)證了該方法的有效性。

2 陣列3維SAR原理

陣列3維SAR系統(tǒng)在機(jī)翼上安裝線陣天線，通過陣列天線的相對(duì)運(yùn)動(dòng)合成2維虛擬面陣，以獲得2維空間分辨率，并結(jié)合SAR雷達(dá)的距離分辨率實(shí)現(xiàn)對(duì)觀測(cè)區(qū)域的3維成像，如圖1所示。

圖1 陣列3維SAR工作幾何Fig.1 Geometry of LASAR

陣列3維SAR的高度向分辨率由脈沖壓縮技術(shù)獲得，根據(jù)脈沖壓縮理論，其分辨率為c/(2B)，其中，c為光速，B為發(fā)射信號(hào)帶寬。陣列3維SAR的沿航跡向分辨率通過平臺(tái)運(yùn)動(dòng)構(gòu)成的合成孔徑獲得，根據(jù)合成孔徑理論，其分辨率為λ/(2θ)，其中，λ為波長，θ為合成孔徑角度。陣列 3維 SAR的切航跡向分辨率通過在機(jī)翼上布設(shè)陣列天線獲得，其分辨率為 λR/(κA)，其中，R為雷達(dá)到目標(biāo)的距離，A 為陣列長度，κ為常數(shù)，1＜κ＜2。假設(shè)雷達(dá)工作于Ka波段，陣列長度10 m，飛行高度10 km，則可得到系統(tǒng)陣列向分辨率約為4～8 m，遠(yuǎn)遠(yuǎn)低于距離向和航跡向分辨率。

在成像處理算法方面，由于陣列3維SAR的等效2維天線陣列結(jié)構(gòu)復(fù)雜，RD算法、ωK算法適用性較差，無法滿足任意模式3維SAR成像處理要求。因此，陣列3維SAR成像處理主要采用后向投影算法等時(shí)域處理技術(shù)。

采用后向投影算法進(jìn)行成像處理后得到的3維SAR圖像如圖2所示，其中，x方向?yàn)殛嚵蟹较颍瑈方向?yàn)槠脚_(tái)運(yùn)動(dòng)方向，z方向?yàn)楦叨认颉？梢钥闯?，其能反映出目?biāo)的典型幾何特征。但是，觀察其單個(gè)高度切片的圖像可以看出，其圖像在陣列方向分辨率較低，且由于稀疏陣列的采用，導(dǎo)致圖像像素間的串?dāng)_較大。

3 最優(yōu)DEM重建的變分模型

3.1 變分模型

在地形測(cè)繪任務(wù)中，陣列3維SAR主要觀測(cè)目標(biāo)為起伏地形，可寫成散射系數(shù)與高程相對(duì)于水平面坐標(biāo)系的函數(shù)的形式，即：

其中，x,y表示 2維水平面的變量，z(x,y)反映了x,y處地形的高度值，σ(x,y)反映了x,y處的復(fù)值散射系數(shù)。

假設(shè)陣列 3維 SAR系統(tǒng)的 3維模糊函數(shù)為χ(x,y,z)，則采用后向投影算法得到的成像結(jié)果可由高程函數(shù)z(x,y)和散射系數(shù)函數(shù)σ(x,y)表示為：

圖2 陣列3維SAR成像處理結(jié)果Fig.2 Imaging results of LASAR

其中，u,v表示2維水平面的變量，I(x,y,z)表示利用系統(tǒng)模糊函數(shù)和散射系數(shù)估計(jì)得到的 3維 SAR圖像。

利用式(2)和實(shí)際系統(tǒng)成像結(jié)果，可建立式(3)所示最優(yōu)化問題模型：

其中，e(x,y,z)表示實(shí)際 3維圖像與估計(jì)圖像間的誤差，表示實(shí)際系統(tǒng)得到的3維SAR圖像。

式(3)表明，陣列3維SAR圖像增強(qiáng)問題的可描述為針對(duì)z(x,y)和σ(x,y)的最優(yōu)化問題：尋找高程函數(shù)z(x,y)和散射系數(shù)函數(shù)σ(x,y)，使得該高程函數(shù)和散射系數(shù)函數(shù)對(duì)應(yīng)的估計(jì)結(jié)果與實(shí)際測(cè)量結(jié)果的誤差最小。由于該最優(yōu)化問題的優(yōu)化變量為關(guān)于自由變量x和y的函數(shù)，其在數(shù)學(xué)上屬于典型的變分問題。

為了便于分析與算法設(shè)計(jì)，需要進(jìn)一步將3維模糊函數(shù) χ(x,y,z)寫成3個(gè)1維模糊函數(shù)乘積的形式，即：

其中，χR(?),χL(?)和χA(?)分別為距離向、陣列方向和航跡向的模糊函數(shù)。一般情況下，距離向和航跡向的系統(tǒng)分辨率優(yōu)于或等于3維圖像像素間隔，χR(?)和χA(?)可近似看作沖激函數(shù)，則式(2)可近似表示為：

此時(shí)，式(3)的最優(yōu)化問題可沿航跡向劃分為一組相互獨(dú)立的最優(yōu)化問題：

其中，“/y”表示“給定航跡向距離單元y”，χ(x)=χL(x)。式(6)將二元變分優(yōu)化問題轉(zhuǎn)化為一組一元變分問題，大大降低了算法實(shí)現(xiàn)難度。

觀察式(6)可以發(fā)現(xiàn)，一方面，由于最優(yōu)化變量z(x)與模糊函數(shù)復(fù)合，該最優(yōu)化問題屬于典型的非線性變分，采用歐拉-拉格朗日方程法求解較為困難。另一方面，在z(x)已知條件下，關(guān)于散射系數(shù)σ(x,y)的最優(yōu)化問題屬于典型的最小二乘問題，可以采用偽逆、稀疏重建等方法進(jìn)行求解。

因此，針對(duì)式(6)的最優(yōu)化問題可分成兩個(gè)階段：首先，給定z(x)，求解最小二乘問題，得到對(duì)應(yīng)的最優(yōu)化散射系數(shù)和最優(yōu)化代價(jià)函數(shù)值；其次，搜索所有可能z(x)，比較最優(yōu)化代價(jià)函數(shù)值，選擇最小的一個(gè)作為最優(yōu)化z(x)其對(duì)應(yīng)的最優(yōu)化散射系數(shù)作為期望的散射系數(shù)。

3.2 局部化近似

為了降低求解該變分問題的難度，首先需要對(duì)其進(jìn)行局部化近似處理，其核心是忽略陣列向模糊函數(shù)的旁瓣(在超分辨問題中，像素間隔遠(yuǎn)遠(yuǎn)小于陣列向分辨率，故散射點(diǎn)間的旁瓣串?dāng)_影響遠(yuǎn)遠(yuǎn)小于臨近區(qū)域散射點(diǎn)間的主瓣串?dāng)_)。在此近似條件下，某個(gè)x位置處高度z和其對(duì)應(yīng)的散射系數(shù)的改變只對(duì)其臨近區(qū)域的數(shù)據(jù)產(chǎn)生影響，如圖3所示。

圖3 DEM增強(qiáng)問題局部化Fig.3 Localization of the DEM enhancement problem

假設(shè)主瓣寬度為 2L+1，當(dāng)前狀態(tài)為 x，則改變當(dāng)前位置的值將導(dǎo)致x?L到x+L段數(shù)據(jù)值(包含所有高度)發(fā)生改變。相應(yīng)地，x?L到x+L數(shù)據(jù)段值的改變則會(huì)對(duì)x?2L到x+2L段散射系數(shù)產(chǎn)生影響。基于此特征，滿足最優(yōu)化問題式(6)的x處最優(yōu)高度z，同樣應(yīng)滿足式(7)所示最優(yōu)化問題：

式(7)的物理意義為：尋找以x為中心的4L+1點(diǎn)最優(yōu)路徑使得以x為中心的2L+1點(diǎn)數(shù)據(jù)對(duì)應(yīng)的均方誤差最小。通過進(jìn)一步分析可以發(fā)現(xiàn)，在無誤差條件下，式(7)局部化問題的最優(yōu)路徑同樣也是最優(yōu)化問題式(6)的最優(yōu)解。

為了便于計(jì)算最小二乘問題，需要通過數(shù)據(jù)重排技術(shù)，進(jìn)一步將式(7)表示為矩陣形式，該過程可以通過簡(jiǎn)單的代數(shù)知識(shí)解決，此處不再贅述。當(dāng)L=2時(shí)，給定高度h其對(duì)應(yīng)的矩陣公式為：

需要注意的是，此處的散射系數(shù)σ(u,h)為 x-z切片內(nèi)的散射系數(shù)，其為x和高度h的函數(shù)，而非式(1)中x與y的函數(shù)。此時(shí)，σ(u,h)具有稀疏性，即 ：σ(x,h)=0,h ≠ z(x)。 一 般 情 況 下 ，φ 為(2L+1)×(4L+1)矩陣，σ[h]為 4L+1 點(diǎn)列向量，d[h]為2L+1 點(diǎn)列向量。

對(duì)于整個(gè)2維相關(guān)數(shù)據(jù)段區(qū)域，其對(duì)應(yīng)的矩陣表示為：

其中，Nz為3維圖像沿z方向的像素點(diǎn)數(shù)。需要注意的是，根據(jù)其物理意義，σ只包含 4L+1 個(gè)非零分量。另外，如果事先已知σ的某些分量為0，則可以去除式(9b)中的對(duì)應(yīng)行，以簡(jiǎn)化最小二乘問題。

4 基于OMP算法的變分模型求解

為了求解式(6)的全局最優(yōu)化問題，需要消減其可行解集以達(dá)到降低運(yùn)算量的目的。本文方法的基本策略是將全局最優(yōu)化問題按照陣列方向分解為一系列的局部最優(yōu)化問題，每個(gè)階段只決定一個(gè)x位置對(duì)應(yīng)的最優(yōu)高度，并通過滑窗迭代的結(jié)構(gòu)實(shí)現(xiàn)對(duì)整個(gè)x-z切片的最優(yōu)化搜索。

通過觀察可以發(fā)現(xiàn)，由于滑窗結(jié)構(gòu)的采用，當(dāng)優(yōu)化x位置時(shí)，其先前位置的最優(yōu)化高程已經(jīng)獲得，而無需進(jìn)行重復(fù)搜索。另一方面，x位置的后續(xù)位置只用于估計(jì)當(dāng)前局部最優(yōu)化問題的最小均方誤差，其各位置的高程將在后續(xù)迭代中計(jì)算，因此可以放松后續(xù)位置的“單值映射條件”。

基于上述考慮，式(7)可轉(zhuǎn)化為改進(jìn)的稀疏重建問題，即：

其中，第1個(gè)約束條件表示位置x?2L到x?1對(duì)應(yīng)的高程為先前獲得的最優(yōu)化路徑； Θ[x]為一個(gè)較大的復(fù)數(shù)，表示x狀態(tài)處存在一個(gè)“顯著的”散射系數(shù)；Sp(σ)表示σ的稀疏度，第3個(gè)約束條件表示整個(gè)優(yōu)化問題的稀疏度為4L+1。

該最優(yōu)化問題可以通過對(duì)OMP算法進(jìn)行改進(jìn)求解。第1個(gè)約束條件表明在進(jìn)行OMP算法迭代初始，需要將x?2L 到x?1位置處已知高度對(duì)應(yīng)于矩陣Ψ的指標(biāo)集作為初始支撐集傳遞給OMP算法；第3個(gè)約束條件表明，OMP算法的迭代次數(shù)應(yīng)為2L次(前 2L個(gè)高程已知，x處對(duì)應(yīng)的高程需要遍歷，在單次OMP調(diào)用時(shí)也相當(dāng)于已知)。由于Θ[x]未知，且其取值隨著x的變化而變化，第2個(gè)約束條件處理起來較為困難，本文采用式(11)估計(jì)式(10)的最優(yōu)代價(jià)函數(shù)(詳見附錄)。

至此，基于OMP算法的變分模型求解算法如下：

目標(biāo)：求解式(7)定義的局部最優(yōu)化問題。

輸入：給定矩陣Ψ，先前高程z(u)，觀測(cè)數(shù)據(jù)d,L和當(dāng)前狀態(tài)x。

步驟1 高程遍歷循環(huán)：遍歷當(dāng)前位置的所有可能高程，并進(jìn)行如下操作：

(1) 選擇z(x)，先前高程以及所有后續(xù)位置所對(duì)應(yīng)的矩陣Ψ的行向量，獲得Ψ的子矩陣；

(2) 記錄z(x)和先前高程所對(duì)應(yīng)的Ψ的子矩陣的行序號(hào)，作為OMP算法的初始支撐集；

(3) 調(diào)用OMP算法；

(4) 記錄OMP算法的殘留誤差和z(x)對(duì)應(yīng)的散射系數(shù)σ(z)。

步驟 2 尋找σ(z)的最大值，并利用式(11)計(jì)算J～[x](z)。

步驟 3 尋找J～[x](z)的最小值對(duì)應(yīng)的位置作為當(dāng)前位置x的最優(yōu)高程。

輸出：當(dāng)前位置x的最優(yōu)高程。

通過遍歷所有沿陣列方向的位置x，即可得到當(dāng)前x-z切片內(nèi)的最優(yōu)高程函數(shù)。

根據(jù)上述分析，基于OMP算法的變分模型求解方法的運(yùn)算量為：

5 性能分析

由于超分辨率條件下圖像增強(qiáng)問題屬于典型的病態(tài)問題，其性能與信噪比密切相關(guān)。當(dāng)無噪聲條件下，很多算法都能獲得較好的增強(qiáng)效果；但隨著噪聲的增加，不同算法的性能差異很大。本節(jié)仿真實(shí)驗(yàn)在固定超分辨倍數(shù)的條件下，分析不同算法的重建誤差，并對(duì)基于變分模型的DEM重建算法的噪聲性能進(jìn)行分析比較。

仿真地形包括起伏山區(qū)地形和城市區(qū)域兩種，圖像尺寸為150像素×200像素。根據(jù)前面的分析可知，由于不同沿航跡單元內(nèi)距離-陣列切片分的問題相互獨(dú)立。本仿真也可以看作在相同信噪比條件下進(jìn)行的不同地形200次蒙特卡洛實(shí)驗(yàn)。為了便于對(duì)比，標(biāo)準(zhǔn)的OMP算法[10－12]和LASSO算法[13－15]也被用于陣列3維SAR圖像增強(qiáng)，其基本過程是抽取3維SAR圖像中每個(gè)航跡向和距離向柵格的陣列向數(shù)據(jù)，調(diào)用稀疏重建方法得到處理結(jié)果，然后對(duì)選擇每個(gè)水平面柵格節(jié)點(diǎn)處最大值作為該處高程的估計(jì)。

假設(shè)L=4，散射系數(shù)幅度為1(隨機(jī)相位)，噪聲標(biāo)準(zhǔn)差為0.1，變分模型方法、OMP和LASSO算法對(duì)起伏地形的處理結(jié)果如圖4所示。對(duì)比可以看出，采用變分模型的處理方法可以在較大噪聲條件下實(shí)現(xiàn)起伏地形的精確 DEM 增強(qiáng)，其對(duì)應(yīng)的DEM殘差圖如圖4(d)所示，其中只包含155個(gè)誤差點(diǎn)。與之相比，OMP算法和LASSO算法所得到的處理結(jié)果則存在較大的噪聲，其誤差點(diǎn)分別為12691個(gè)和5926個(gè)。

圖4 起伏地形的最優(yōu)DEM重建結(jié)果(L = 4，噪聲標(biāo)準(zhǔn)差0.1)Fig.4 Enhancement results of a mountain area with L = 4 and the standard deviation 0.1

圖5為仿真城市區(qū)域的DEM增強(qiáng)結(jié)果。觀察圖5(a)所示的觀測(cè)地形可以發(fā)現(xiàn)，由于城市區(qū)域房頂?shù)冉ㄖY(jié)構(gòu)與地面平行，其沿陣列方向抽取的數(shù)據(jù)表現(xiàn)出了“連續(xù)性”而非“稀疏性”。此時(shí)采用稀疏重建模型的算法處理將更為困難。圖5(b)為變分模型的處理結(jié)果，可以看出其能夠較好地重建城市區(qū)域地形。與之相反，圖5(e)和圖 5(f)的OMP算法、LASSO算法得到的結(jié)果則由于稀疏性的破壞導(dǎo)致處理效果嚴(yán)重下降。通過進(jìn)一步統(tǒng)計(jì)發(fā)現(xiàn)，變分模型法、OMP和LASSO算法得到的DEM圖像的殘留誤差分別為119,24473和14940。

通過上述分析比較可以看出，由于利用了DEM圖的單值性特征，基于變分模型的DEM重建方法性能要遠(yuǎn)高于基于稀疏信號(hào)模型的重建方法，且在處理起伏地形、城市區(qū)域時(shí)都能得到較好的處理結(jié)果。

6 結(jié)論

本文提出了基于變分模型的陣列3維SAR最優(yōu)DEM重建方法，該方法直接將DEM圖作為最優(yōu)化目標(biāo)，通過尋找最優(yōu)化DEM圖和對(duì)應(yīng)的散射系數(shù)，實(shí)現(xiàn)最小二乘意義下的最優(yōu)DEM重建。仿真分析表明，該方法可以實(shí)現(xiàn)各種地形(山區(qū)、城市)等的穩(wěn)健DEM增強(qiáng)，其性能遠(yuǎn)優(yōu)于OMP算法和正則化方法。

圖5 城市區(qū)域的最優(yōu)DEM重建結(jié)果(L = 4，噪聲標(biāo)準(zhǔn)差0.1)Fig.5 Enhancement results of an urban area with L = 4 and the standard deviation 0.1

附錄

本附錄將推導(dǎo)式(11)，該問題可以重述為：

給定線性方程組Ax=y和其OMP算法解*x，計(jì)算該線性方程組在約束條件xi=c下的殘留誤差，其中，i屬于支撐集。

為了分析上述問題，可將 OMP算法的殘留誤差寫為：

另一方面，Ax=y且xi=c可以寫為：

其中，A↓i表示剔除第i列后A矩陣的子矩陣，x↓i表示剔除第i個(gè)分量后x向量的子向量，αi表示矩陣A的第i列。

由于i屬于支撐集，式(A-3)可表示為：

根據(jù)OMP算法的正交性，則有：

即式(11)。

[1]Matthias W and Markus G.Initial ARTINO radar experiments[C].20108th European Conference on Synthetic Aperture Radar (EUSAR),Aachen,Germany,2010: 1-4.

[2]Han Kuo-ye,Wang Yan-ping,Tan Wei-xian,et al..Efficient pseudopolarformat algorithm for down-looking linear-array SAR 3-D imaging[J].IEEE Geoscience and Remote Sensing Letters,2015,12(3): 572-576.

[3]Peng Xue-ming,Hong Wen,Wang Yan-ping,et al..Downward looking linear array 3D SAR sparse imaging with wave-front curvature compensation[C]. International Conference on Signal Processing,Communication and Computing (ICSPCC),Kunming,2013: 1-4.

[4]Zhang S,Zhu Y,and Kuang G.Imaging of downward-looking linear array three-dimensional SAR based on FFT-MUSIC[J].IEEE Geoscience and Remote Sensing Letters,2015,12(4):885-889.

[5]Peng Xue-ming,Hong Wen,Wang Yan-ping,et al..Polar format imaging algorithm with wave-front curvature phase error compensation for airborne DLSLA three-dimensional SAR[J].IEEE Geoscience and Remote Sensing Letters,2014,11(6): 1-5.

[6]Wei Shun-jun,Zhang Xiao-ling,and Shi Jun.Sparse autofocus via bayesian learning iterative maximum and applied for LASAR 3-D imaging[C].2014 IEEE Radar Conference,Cincinnat,USA,2014: 666-669.

[7]Zhang Ying-jie,Han Kuo-ye,Wang Yan-ping,et al..Study on motion compensation for airborne forward looking array SAR by time division multiplexing receiving[C].2013 Asia-Pacific Synthetic Aperture Radar (APSAR),Tsukuba,Japan,2013:392-395.

[8]Zhuge Xiao-dong and Yarovoy A G.A sparse aperture MIMO-SAR-Based UWB imaging system for concealed weapon detection[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(1): 509-518.

[9]Shi Jun,Zhang Xiao-ling,Xiang Gao,et al..Signal processing for microwave array imaging: TDC and sparse recovery[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(11): 4584-4598.

[10]Needell D and Tropp J A.CoSaMP: iterative signal recovery from incomplete and inaccurate samples[J].Applied and Computational Harmonic Analysis,2009,26(3): 301-321.

[11]Wang Jian,Kwon S,and Shim B.Generalized orthogonal matching pursuit[J].IEEE Transactions on Signal Processing,2012,60(12): 6202-6216.

[12]Wang Jian and Shim B.On the recovery limit of sparse signals using orthogonal matching pursuit[J].IEEE Transactions on Signal Processing,2012,60(9): 4973-4976.

[13]Tibshirani R.Regression shrinkage and selection via the lasso[J].Journal Of The Royal Statistical Society Series B-Methodological,1996,58(1): 267-288.

[14]EfronB,Hastie T,Johnstone I,et al..Least angle regression[J].Annals of Statistics,2004,32(2): 407-499.

[15]Rosset S and Zhu Ji.Piecewise linear regularized solution paths[J].Annals of Statistics,2007,35(3): 1012-1030.