基于拉普拉斯算子的綜合孔徑輻射計成像算法

楊曉城，楊真乙，閻敬業(yè)，武 林，蔣明峰

（1. 浙江理工大學(xué)信息學(xué)院，浙江杭州 310018；2. 中國科學(xué)院國家空間科學(xué)中心，北京 100190）

1 引言

作為被動微波遙感領(lǐng)域的一項新技術(shù)，綜合孔徑微波輻射計利用小口徑天線干涉綜合成大口徑天線的方法，有效提高了對地觀測的空間分辨率，同時避免了大口徑天線帶來的折疊、機(jī)械掃描和重量等方面的困難. 目前，綜合孔徑微波輻射計逐漸成熟，被應(yīng)用于陸地、海洋和大氣遙感以及人體安檢等領(lǐng)域［1，2］.

成像反演是綜合孔徑微波輻射計的一項關(guān)鍵內(nèi)容［3，4］. 然而，綜合孔徑輻射計成像反演在數(shù)學(xué)上是病態(tài)的反問題，需要進(jìn)行正則化處理以獲得穩(wěn)定的唯一解. 目前常用的正則化方法可歸為2 類：第一類是采用數(shù)學(xué)的手段進(jìn)行約束，如Tikhonov 正則化、最小范數(shù)正則化、共軛梯度最小二乘（Conjugate Gradients Least Squares，CGLS）等［5～7］；第二類方法是通過綜合孔徑輻射計自身的物理特性進(jìn)行約束，即帶限正則化法［8］. 雖然上述正則化方法能夠有效克服其病態(tài)性，但依然存在較大的重構(gòu)誤差［9］.

值得注意的是，傳統(tǒng)的Tikhonov 正則化和最小范數(shù)正則化方法，將固定的約束信息添加到了所有的頻帶上，這會導(dǎo)致重構(gòu)得到的亮溫解過于平滑，從而無法保留原始亮溫分布的所有特征. 因此，借鑒雙參數(shù)正則化［10］的思想，本文通過在正則化過程中引入一項新的拉普拉斯算子項，提出了一種拉普拉斯混合正則化方法，并基于研制的L 波段FPIR 樣機(jī)進(jìn)行了仿真驗證.

2 綜合孔徑輻射計成像原理

與真實孔徑輻射計不同，綜合孔徑輻射計干涉測量結(jié)果是觀測場景亮溫分布在頻率域上的可見度函數(shù)采樣. 可見度函數(shù)V（u）與觀測場景亮溫分布TB（ξ）之間的關(guān)系可以表示為［11，12］

其中，ξ=（sinθcosφ，sinθsinφ）為在笛卡爾坐標(biāo)系下方向余弦矢量，ukl為與天線k和l對應(yīng)的基線矢量；Ωk和Ωl分別為天線k和l的等效立體角，F(xiàn)k和Fl分別為天線k和l歸一化電壓方向圖（表示Fl的復(fù)共軛）；Tr為接收機(jī)的物理溫度，kl(τ)為fringe washing 函數(shù)，τ=-ukl ξ/f0為時間延遲，f0為中心頻率.

將式（1）進(jìn)行離散化，并考慮實際測量過程中可見度函數(shù)測量誤差和噪聲的影響［13～15］，可得

其中，G為系統(tǒng)調(diào)制矩陣，V為可見度函數(shù)向量，T為觀測亮溫向量，n為噪聲向量. 由于亮溫采樣點(diǎn)數(shù)M一般大于可見度函數(shù)采樣點(diǎn)數(shù)N，因此式（2）是欠定的. 同時由于G矩陣的不適定性，噪聲n的微小變化都會使得反問題解產(chǎn)生很大的擾動，因此矩陣方程是病態(tài)的，其解既不唯一也不穩(wěn)定.

對于病態(tài)的線性方程，通常是進(jìn)行正則化處理. 正則化的作用是通過增加有關(guān)解空間的附加信息，并利用約束泛函獲得穩(wěn)定可行解，從而緩解反問題的病態(tài)性. 目前，綜合孔徑輻射計中廣泛使用的是最小范數(shù)正則化和帶限正則化方法.

最小范數(shù)正則化是指求解式（2）的最小范數(shù)解，即

為了求解式（3），先對矩陣G進(jìn)行奇異值分解，即

其中，si≥0為G矩陣按降序排列的奇異值，向量ui和vi分別為si的左奇異向量和右奇異向量. 然后，通過拋棄那些對誤差傳遞影響很大的較小奇異值來緩解矩陣G的病態(tài)性，即采用一個秩虧缺的良態(tài)矩陣Gk=來代替矩陣G. 因此，反演亮溫為

其中，k為正則化參數(shù).

對于綜合孔徑輻射計，由于天線單元數(shù)量是有限的，因此測量的可見度函數(shù)采樣值被限制在有限的帶寬H內(nèi). 帶限正則化是指求解下列約束最小二乘問題：

其中，T′=UT為T的傅里葉變換，U是離散傅里葉變換運(yùn)算符. 因此，反演亮溫為

其中，J+=（J*J）-1J*是矩陣J=GU*Z的廣義逆矩陣，Z是補(bǔ)零運(yùn)算符，*表示共軛轉(zhuǎn)置.

由于綜合孔徑輻射計覆蓋的基線范圍有限，無法獲取觀測目標(biāo)亮溫的高頻分量，因此存在固有的系統(tǒng)誤差. 為了降低系統(tǒng)誤差，引入一個先驗亮溫分布，對差分可見度函數(shù)ΔV=V-進(jìn)行反演成像，即求解最小二乘問題：

其中，ΔT=T-. 為了獲得先驗亮溫分布，首先利用陸地與海洋分界的先驗信息，將觀測場景劃分為陸地區(qū)域和海洋區(qū)域，并將每個區(qū)域的亮溫設(shè)為恒定值，然后根據(jù)最小二乘準(zhǔn)則求出每個區(qū)域的亮溫值，詳見文獻(xiàn)［16］. 因此，對于最小范數(shù)正則化，其反演亮溫變?yōu)?/p>

對于帶限正則化，其反演亮溫變?yōu)?/p>

為了方便區(qū)分，式（9）的方法稱為改進(jìn)的最小范數(shù)正則化，式（10）的方法稱為改進(jìn)的帶限正則化.

3 拉普拉斯混合正則化

對于雙參數(shù)混合正則化方法，王文娟等將其應(yīng)用于電導(dǎo)率的反演成像中，得到了更為精確穩(wěn)定的解［17］；劉婷等將其應(yīng)用于腦磁時序信號的逆問題求解中，重構(gòu)出更為準(zhǔn)確且時域平滑的腦內(nèi)神經(jīng)信號［18］. 受此啟發(fā)，為了進(jìn)一步降低綜合孔徑輻射計的重構(gòu)誤差，本文借鑒雙參數(shù)正則化的思想，提出了一種拉普拉斯混合正則化方法.

3.1 拉普拉斯混合正則化方法

為了獲取高精度的觀測場景亮溫，本文在標(biāo)準(zhǔn)的單參數(shù)Tikhonov 正則化的基礎(chǔ)上，通過引入新的拉普拉斯算子項，對差分亮溫添加2 種不同的約束，即求解如下極小值問題：

其中，‖ΔV-GΔT‖2為數(shù)據(jù)擬合項，即解的殘差范數(shù)；‖ΔT‖2為能量泛函約束項，使得解全局能量最??；‖LΔT‖2為拉普拉斯算子項，使得重構(gòu)解能夠保留原始亮溫分布的多尺度信息，其中L是離散化后與拉普拉斯算子對應(yīng)的正則化矩陣；λ1和λ2為大于0 的正則化參數(shù). 拉普拉斯算子可以分為一階拉普拉斯算子和二階拉普拉斯算子等，其中一階拉普拉斯算子可以反映圖像的邊緣特征，二階拉普拉斯算子可以反映圖像的紋理結(jié)構(gòu)特征. 本文采用二階拉普拉斯算子，則離散化后對應(yīng)的正則化矩陣L為

求解式（11），可得反演亮溫為

3.2 正則化參數(shù)的選擇

正則化參數(shù)的選取在綜合孔徑輻射計成像反問題研究中起著非常重要的作用. 然而，正則化參數(shù)的選取策略較為復(fù)雜而且沒有普適的方法. 目前，主要采用的2 種方法，分別是L 曲線和廣義交叉驗證（Generalized Cross-Validation，GCV）準(zhǔn)則.

對于最小范數(shù)正則化，將采用GCV 準(zhǔn)則選取正則化參數(shù).GCV 準(zhǔn)則的基本思想［19］是當(dāng)原始數(shù)據(jù)項V中去掉一項Vi時，由此產(chǎn)生的新方程可以預(yù)測出去掉的分量，即通過最小化GCV函數(shù)來確定正則化參數(shù)，其中f(k)=GG+k，tr表示方程的跡.

對于拉普拉斯混合正則化，其有2個正則化參數(shù)需要確定，通過傳統(tǒng)的GCV準(zhǔn)則無法同時得到最優(yōu)的2個參數(shù)，所以本文采用多維擴(kuò)展后的GCV準(zhǔn)則. 多維擴(kuò)展后的GCV 準(zhǔn)則［20］將GCV 函數(shù)視為向量λ=（λ1，λ2）T的多元函數(shù)，即

其中f（λ）=G（G*G+λ1I+λ2LTL）-1G*. 然而，通過最小化式（15）求解正則化參數(shù)計算代價太大，難以同時處理矩陣G和所有正則化矩陣. 文獻(xiàn)［20］證明了可以利用矩陣Z(λ)=(vi(λi))2來代替v（λ），通過求解目標(biāo)函數(shù)Z（λ）的極小值來獲取近似解. 其求解過程為分別對單參數(shù)的GCV 函數(shù)vi（λi）進(jìn)行最小化獲得正則化參數(shù)λ1和λ2，然后組成向量λ.

最小范數(shù)正則化和帶限正則化等傳統(tǒng)正則化方法對成像反問題添加單一的約束和限制，屬于單一正則化方法. 然而，綜合孔徑輻射計成像反問題是非常復(fù)雜的，單一正則化方法都有一定的側(cè)重性，具有相應(yīng)的優(yōu)點(diǎn)和缺點(diǎn)，這意味著即使在選擇了最優(yōu)的正則化參數(shù)情況下，單一正則化也無法確保反問題解的穩(wěn)定性和健全性. 因此，與傳統(tǒng)的單一正則化不同，本文提出了一種拉普拉斯混合正則化方法，其優(yōu)點(diǎn)在于通過引入拉普拉斯算子項來反映原始亮溫分布的其他尺度特征，如邊緣特征、紋理結(jié)構(gòu)特征等，從而提高反演亮溫的精度.

為了降低頻譜截斷所導(dǎo)致的吉布斯振蕩現(xiàn)象，一般需要利用窗函數(shù)對可見度函數(shù)進(jìn)行加窗平滑，可得加窗后重構(gòu)亮溫為

其中，W為窗函數(shù).

為了定量地分析重構(gòu)誤差，計算重構(gòu)結(jié)果與真實亮溫分布之間的均方根誤差RMSE 與峰值信噪比PSNR，即

其中，m為無混疊視場內(nèi)的像素個數(shù).

4 仿真結(jié)果

為了驗證上述拉普拉斯混合正則化方法的有效性，基于L 波段（1.4 GHz）FPIR 樣機(jī)系統(tǒng)進(jìn)行了仿真分析.FPIR 樣機(jī)是一套L/X 波段一維全極化綜合孔徑輻射計，主要用于探測海表面鹽度［21］. 圖1為FPIR樣機(jī)雙頻天線陣列布局，其中L 波段的陣列天線放置在兩邊，X 波段的陣列天線放置在中間區(qū)域. 仿真系統(tǒng)為工作在H 極化狀態(tài)的L 波段FPIR 樣機(jī)，最短基線為d=0.589λ0，6 根觀測天線a1，a2，…，a6 稀疏地排列在不同的位置，其間隔依次序分別為d、d、2d、5d、d，所以可得最長基線為10d.

圖1 FPIR雙頻天線陣列布局

仿真系統(tǒng)其他參數(shù)如下：帶寬為B=20 MHz，fringe washing 函數(shù)設(shè)為sinc（Bτ），積分時間為1 s. 包含零基線并考慮共軛對稱性后，總共可獲得31 個可見度函數(shù)采樣點(diǎn). 在天線暗室中，對L 波段的陣列中6 根觀測天線的方向圖進(jìn)行測量，綜合孔徑方向上測量的歸一化天線電壓方向圖如圖2所示.

圖2 綜合孔徑方向歸一化天線方向圖

原始亮溫分布如圖3 所示，包括2 種典型的觀測目標(biāo)區(qū)域，其中圖3（a）為純海洋區(qū)域，圖3（b）為海洋陸地交界區(qū)域. 原始亮溫分布來源于Aquarius 衛(wèi)星2013 年8 月19 日觀測的某區(qū)域L 波段H 極化輻射亮溫數(shù)據(jù).作為一維綜合孔徑輻射計，F(xiàn)PIR 系統(tǒng)單次測量只能獲取一個條帶的亮溫數(shù)據(jù)，通過載荷平臺的移動進(jìn)行推掃成像，以獲取二維亮溫圖像.

圖3 原始亮溫分布

基于上述仿真系統(tǒng)，利用式（1）可仿真得到可見度函數(shù). 在實際測量過程中，不可避免地存在測量噪聲，所以利用均值為零方差為0.01×max（Vi）的高斯白噪聲來模擬測量誤差或者噪聲. 然后，分別利用改進(jìn)的最小范數(shù)正則化、改進(jìn)的帶限正則化和拉普拉斯混合正則化對含噪聲的可見度函數(shù)進(jìn)行重構(gòu)，得到無混疊視場內(nèi)重構(gòu)結(jié)果，如圖4 所示，其中圖4（a）為純海洋區(qū)域重構(gòu)結(jié)果，圖4（b）為海陸交界區(qū)域重構(gòu)結(jié)果. 在圖4 中使用的窗函數(shù)為漢寧窗，原始亮溫分布也進(jìn)行了加窗處理以便于比較. 從圖4 中可以看出，與改進(jìn)的最小范數(shù)正則化和帶限正則化相比，拉普拉斯混合正則化重構(gòu)誤差明顯降低，更加接近于真實亮溫分布.

此外，計算了重構(gòu)結(jié)果與原始圖像之間的RMSE和PSNR. 表1 為圖4 中純海洋區(qū)域和海陸交界區(qū)域亮溫分布經(jīng)過不同成像算法后重構(gòu)誤差對比結(jié)果.

表1 不同成像算法的重構(gòu)誤差

圖4 重構(gòu)結(jié)果

為了進(jìn)一步驗證拉普拉斯混合正則化的魯棒性，分析了噪聲對重構(gòu)結(jié)果的影響. 在式（1）可見度函數(shù)基礎(chǔ)上加入不同方差的高斯白噪聲，然后分別利用改進(jìn)的最小范數(shù)正則化、改進(jìn)的帶限正則化和拉普拉斯混合正則化進(jìn)行成像反演. 圖5 所示為不同噪聲情況下純海洋區(qū)域重構(gòu)結(jié)果性能. 在圖5 中，當(dāng)噪聲為零時，改進(jìn)的最小范數(shù)、改進(jìn)的帶限正則化和拉普拉斯混合正則化的RMSE 分別為0.05 K，0.06 K 和0.03 K. 從圖5中可以看出，雖然測量噪聲大小不同，但是拉普拉斯混合正則化比改進(jìn)的最小范數(shù)和帶限正則化重構(gòu)結(jié)果RMSE 降低了40%以上，PSNR 提高了5 dB 以上. 值得注意的是，仿真中帶限正則化的性能明顯差于最小范數(shù)正則化，這是由于FPIR 樣機(jī)天線單元數(shù)量僅有6 根，導(dǎo)致頻率上覆蓋的基線范圍很小，從而導(dǎo)致帶寬受限的約束性能不佳.

圖5 不同噪聲情況下純海洋區(qū)域重構(gòu)結(jié)果性能

類似地，不同噪聲情況下海陸交界區(qū)域重構(gòu)結(jié)果性能如圖6 所示. 從圖6 中可以看出，與改進(jìn)的最小范數(shù)和帶限正則化相比，在不同噪聲情況下，拉普拉斯混合正則化重構(gòu)結(jié)果RMSE 有效降低，PSNR 明顯上升.此外，在圖6中，當(dāng)噪聲為零時，改進(jìn)的最小范數(shù)和帶限正則化以及拉普拉斯混合正則化的RMSE 分別為1.46 K，1.64 K 和1.18 K. 對于海陸交界區(qū)域，由于觀測亮溫分布跨度大，其吉布斯振蕩現(xiàn)象更加明顯，殘留的重構(gòu)誤差較大. 因此，對于海陸交界等亮溫變化劇烈的場景，拉普拉斯混合正則化性能比改進(jìn)的最小范數(shù)和帶限正則化有一定的提升.

圖6 不同噪聲情況下海陸交界區(qū)域重構(gòu)結(jié)果性能

綜上所述，本文提出的拉普拉斯混合正則化方法對于測量誤差或者噪聲干擾具有良好的魯棒性，對于不同的觀測場景都能有效地降低其重構(gòu)誤差，尤其是對于純海洋區(qū)域.

拉普拉斯混合正則化的計算時間與矩陣G的大小密切相關(guān). 當(dāng)矩陣G大小為31×91 時，在2.4 GHz Intel i5-4210U 處理器和12 GB 內(nèi)存的筆記本電腦上，改進(jìn)的最小范數(shù)正則化、改進(jìn)的帶限正則化和混合正則化的Matlab 計算時間分別為0.159 6 s，0.112 6 s 和0.318 3 s.結(jié)果表明，雖然拉普拉斯混合正則化能有效地提高反演的精度，但是需要花費(fèi)更長的計算時間，其計算時間大約是改進(jìn)的最小范數(shù)正則化的兩倍. 對于大型的綜合孔徑天線陣列，可利用GPU 平臺來計算拉普拉斯混合正則化反演結(jié)果以加快計算速度.

5 結(jié)論

高效的成像反演算法是綜合孔徑輻射計的關(guān)鍵技術(shù)之一. 綜合孔徑輻射計干涉輸出的是可見度函數(shù)，成像反演算法的目的是將可見度函數(shù)重構(gòu)為精確的輻射亮溫圖像. 然而，綜合孔徑輻射計成像反問題是病態(tài)的，雖然傳統(tǒng)的正則化方法能有效克服其病態(tài)性，但依然存在較大的重構(gòu)誤差. 針對上述問題，提出了一種基于拉普拉斯算子的綜合孔徑輻射計成像算法. 該方法通過引入新的拉普拉斯算子項，對差分亮溫構(gòu)造拉普拉斯混合正則化，并利用多維擴(kuò)展的廣義交叉驗證準(zhǔn)則選取2 個正則化參數(shù). 與傳統(tǒng)的單一正則化相比，其優(yōu)點(diǎn)在于通過添加2 個約束項來構(gòu)造拉普拉斯混合正則化，以反映原始亮溫分布的多尺度特征，從而得到更精確的反演亮溫結(jié)果. 利用研制的L 波段FPIR 樣機(jī)進(jìn)行了仿真分析，仿真結(jié)果表明：與改進(jìn)的最小范數(shù)正則化和帶限正則化相比，本文方法重構(gòu)結(jié)果更加接近真實亮溫分布，能明顯降低重構(gòu)誤差，并且具有良好的魯棒性.