崔浩貴 劉 濤 蔣宇中 高 俊

1 引言

極化合成孔徑雷達(dá)(Polarimetric Synthetic Aperture Radar, PolSAR)具有全天時(shí)全天候的觀測能力和多維信息遙感能力，是重要的遙感信息源。雜波統(tǒng)計(jì)建模及其參數(shù)估計(jì)方法是PolSAR圖像目標(biāo)檢測和分類識別等圖像解譯手段的基礎(chǔ)課題[1,2]。在高分辨情形下，地面或者海面散射雜波一般服從某種非高斯分布。針對PolSAR圖像雜波非高斯建模問題，目前研究最多的是多變量乘積模型。該模型將目標(biāo)的散射表示為一個(gè)表征目標(biāo)雷達(dá)散射截面的紋理分量和一個(gè)表征相干斑的復(fù)高斯矢量的乘積，并且這兩個(gè)變量之間相互統(tǒng)計(jì)獨(dú)立。在多變量模型中，Gamma分布、逆 Gamma分布、Fisher分布和α穩(wěn)定分布等非高斯分布常常用來對紋理分量進(jìn)行建模[37]-。

K分布模型是目前應(yīng)用最廣泛的PolSAR圖像統(tǒng)計(jì)模型，它是在紋理分量為Gamma分布時(shí)根據(jù)多變量乘積模型得到的。尋找快速、準(zhǔn)確的參數(shù)估計(jì)方法是K分布研究的核心問題。傳統(tǒng)的PolSAR圖像K分布參數(shù)估計(jì)方法是利用單極化SAR的方法分別估計(jì)出每個(gè)極化通道強(qiáng)度圖像的參數(shù)，然后將各通道的參數(shù)求平均。常用的單通道K分布參數(shù)估計(jì)方法為矩估計(jì)法[8]和最大似然估計(jì)法[9]。文獻(xiàn)[10]的研究表明在矩估計(jì)法中采用分?jǐn)?shù)階矩能取得更好的估計(jì)效果。文獻(xiàn)[11]提出了基于 zrlg(z)的K分布參數(shù)估計(jì)方法，該方法在 z l g(z)情形下有解析的形式。文獻(xiàn)[12]對文獻(xiàn)[11]的方法在不同r值下的性能進(jìn)行了仿真比較，指出0

上述基于單通道的估計(jì)方法實(shí)際上只用了PolSAR圖像協(xié)方差矩陣的對角元素，沒有用到各極化通道間的相關(guān)信息。在單通道方法的基礎(chǔ)上，針對 K分布下 PolSAR圖像協(xié)方差矩陣的分布特性，以矩陣為變量對參數(shù)估計(jì)方法進(jìn)行了研究。文獻(xiàn)[13]將文獻(xiàn)[10]的方法擴(kuò)展到多維情形，給出了基于協(xié)方差矩陣行列式值二階矩的參數(shù)估計(jì)方法。文獻(xiàn)[14,15]提出了基于協(xié)方差矩陣的對數(shù)累積量(Matrix Log-Cumulants, MLC)的參數(shù)估計(jì)方法，該方法本質(zhì)上是協(xié)方差矩陣行列式值的Mellin變換。研究表明基于對數(shù)累積量的估計(jì)方法在已有方法中是估計(jì)精度最優(yōu)的，并已獲得廣泛應(yīng)用[13,16]。該方法的缺陷在于沒有解析的表達(dá)式，雖然根據(jù)表達(dá)式的單調(diào)性能較快地找到方程的解，但是在處理海量的PolSAR圖像數(shù)據(jù)時(shí)，其運(yùn)算速度仍有一定的局限性。本文將文獻(xiàn)[11]的估計(jì)方法擴(kuò)展到全極化多視數(shù)據(jù)的情形，提出了基于|z |rlg|z|混合矩(ZrLZ)的參數(shù)估計(jì)方法，該方法在r=1/d時(shí)有解析表達(dá)式。并且對ZrLZ方法在r取值不同值時(shí)的估計(jì)性能進(jìn)行了仿真分析。最后用仿真數(shù)據(jù)和實(shí)測數(shù)據(jù)對本文方法與已有方法進(jìn)行了比較，結(jié)果表明 ZrLZ方法對不同α值區(qū)域的適應(yīng)能力更好，并且在r=1/d時(shí)的計(jì)算速度明顯優(yōu)于對數(shù)累積量方法。

2 PolSAR圖像K分布統(tǒng)計(jì)模型

PolSAR圖像中每個(gè)像素點(diǎn)可用目標(biāo)矢量表示，在互易條件下可寫為

多變量乘積模型是由一個(gè)表征目標(biāo)雷達(dá)散射截面積的紋理分量 t和一個(gè)表征相干斑分量的復(fù)高斯矢量y組成，它們是相互獨(dú)立的，即目標(biāo)矢量為

其中復(fù)高斯矢量y的概率密度函數(shù)為

這里d是復(fù)高斯矢量的維數(shù)， Γ = E (y yH)，E(?)表示隨機(jī)變量的數(shù)學(xué)期望。

實(shí)際中往往通過多視處理來抑制相干斑，該過程可由式(4)表示。

這里L(fēng)為視圖數(shù)，即進(jìn)行平均的像素點(diǎn)的個(gè)數(shù)，k代表第k個(gè)像素點(diǎn)。如果視圖數(shù)不是很多，往往假設(shè)多視處理的像素點(diǎn)具有相同的紋理分量，即 ()X k獨(dú)立于k，此時(shí)式(4)可簡化為

其中

式中Y為復(fù)高斯矢量y對應(yīng)的協(xié)方差矩陣，Y服從Wishart分布，即 Y ～ Wd(L,Γ)，其概率密度函數(shù)為

其中d是Y的維數(shù)， Γ = E (y yH)，函數(shù) Γd(L)為復(fù)數(shù)形式的多變量Gamma函數(shù)

當(dāng)式(5)中紋理分量服從均值為1的Gamma分布時(shí)，其乘積模型服從K分布，即 Z ～ K (L, Σ, α )，其概率密度函數(shù)為

其中L為視圖數(shù)， Σ = E (Z), α為Gamma分布的形狀參數(shù)，Kα-Ld(.)為第α-Ld階的第2類修正Bessel函數(shù)。

3 基于 |z |r lg|z |混合矩的 K 分布參數(shù)估計(jì)方法

假定高斯假設(shè)下協(xié)方差矩陣測量值Z和期望值Σ的行列式值分別為γ=Z,Υ=Σ , d是協(xié)方差矩陣的維數(shù)，L是視圖數(shù)。已知隨機(jī)變量η= ( 2L)dγ/Υ 為d個(gè)相互獨(dú)立的 χ2分布變量的乘積構(gòu)成，其自由度分別為2L, 2(L - 1),…,2(L-d+ 1 )[14,15]，考慮如式(5)所示的乘積模型，第 2類(Mellin變換)特征函數(shù)滿足：

根據(jù) χ2分布變量的Mellin變換[14]，可得行列式值隨機(jī)變量 η = ( 2L)dγ/Υ 的Mellin變換為

將式(11)進(jìn)行Mellin逆變換，可知高斯假設(shè)下協(xié)方差矩陣的行列式值的高階矩為

Gamma分布 t ～ γ(1,α)的矩特征為

對于式(5)所示的乘積模型，由于紋理分量與相干斑分量是相互獨(dú)立的，因此其矩特征存在如下關(guān)系：

因此可知K分布下協(xié)方差矩陣行列式值的矩特征為E(| |)r

Z 關(guān)于r求導(dǎo)，可得

將式(15)代入式(16)，可得

其中ψ(0)(x)為Digamma函數(shù)，并且ψ(0)(x)=Γ'( x ) /Γ( x )。對式(17)取不同的r值，并將兩個(gè)式子相減，可以得到α的參數(shù)估計(jì)方法如式(18)所示

其中r1≠r2。當(dāng)r=0時(shí)，式(17)實(shí)際上就是一階對數(shù)累積量的表達(dá)式：

已有的研究表明將式(19)用于參數(shù)估計(jì)，其估計(jì)性能較好[7]，同時(shí)為了便于對估計(jì)方法進(jìn)行分析，這里將 r2的值取0，此時(shí)式(18)可寫為

我們將式(20)的估計(jì)方法稱為基于|z |rlg|z|混合矩的參數(shù)估計(jì)法(ZrLZ)。對式(20)求導(dǎo)可得求解方程具有單調(diào)遞減性，可以利用基于中值的搜索方法進(jìn)行快速求解。另外我們發(fā)現(xiàn)當(dāng)r=1/d時(shí)，由式(20)可得到ZrLZ(r=1/d)估計(jì)α的解析表達(dá)式

式(21)的推導(dǎo)應(yīng)用了等式

4 實(shí)驗(yàn)結(jié)果與分析

這里對以下幾種常用的K分布模型的參數(shù)估計(jì)方法性能進(jìn)行比較：

(1)文獻(xiàn)[15]提出的基于協(xié)方差矩陣行列式值二階對數(shù)累積量的參數(shù)估計(jì)方法(Second-order MLC,SMLC)；

(2)文獻(xiàn)[13]提出的基于二階矩特征(Second MOMent, SMOM)的參數(shù)估計(jì)方法；

(3)文獻(xiàn)[10]提出的基于單極化通道強(qiáng)度數(shù)據(jù)的分?jǐn)?shù)階矩的參數(shù)估計(jì)方法(Fractional MOMent,FMOM)；

(4)本文給出的ZrLZ方法。文獻(xiàn)[12]指出在單極化通道情形下 0

上述幾種方法中，F(xiàn)MOM為基于單極化通道的方法，其它的是基于協(xié)方差矩陣的估計(jì)方法。ZrLZ(r=1/d)和 SMOM 方法具有解析的表達(dá)式。SMLC和ZrLZ方法(r≠1/d)可以用基于中值的搜索方法進(jìn)行求解，這里搜索精度設(shè)為 0.001。FMOM方法的求解方程不具備單調(diào)性，因此用遍歷搜索的方法求解，選取的遍歷范圍為 α ∈ ( 0,100)，搜索精度為0.001。

4.1 仿真數(shù)據(jù)分析

K 分布仿真數(shù)據(jù)由如式(5)所示的乘積模型得到，視圖數(shù)取 L = 1 0。形狀參數(shù)α的取值范圍一般在0.1(極不均勻區(qū)域)到10(海洋、平原等均勻區(qū)域)之間。圖1給出了在不同的α取值下( α ∈ [ 0.1,10])，不同參數(shù)估計(jì)方法在樣本數(shù)為512時(shí)的相對估計(jì)偏差與相對估計(jì)方差，仿真次數(shù)為 N = 1 0000次。其中相對估計(jì)偏差定義為對 | α︿ -α|/α求平均，相對估計(jì)方差為 v ar(α︿ ) /α2。首先比較ZrLZ方法在不同參數(shù)r下的估計(jì)性能，可以看到r=0.60時(shí)的估計(jì)性能較差，而r=0.05, r=0.20以及r=1/d這幾種方法估計(jì)性能非常接近。因此下面在與常用的幾種K分布模型的參數(shù)估計(jì)方法相比較時(shí) ZrLZ方法只討論r=1/d的情形。

同時(shí)從圖1中可以看出當(dāng)α值較小時(shí)，SMOM方法相對估計(jì)偏差和方差都非常大，此時(shí)該方法幾乎失效。FMOM方法在α很小時(shí)效果最好，但是在α較大時(shí)，估計(jì)性能最差。另外可以看出本文的ZrLZ(r=1/d)方法在α值較小時(shí)，估計(jì)性能優(yōu)于SMLC方法。當(dāng)α值增大時(shí)，ZrLZ方法與 SMLC方法性能相當(dāng)。因此可以得出ZrLZ方法對于不同α值的估計(jì)性能要更穩(wěn)健。

圖2為不同估計(jì)方法運(yùn)算時(shí)間的比較，這里數(shù)據(jù)源為α=6的K分布仿真數(shù)據(jù)。計(jì)算機(jī)CPU為Inter E5700，雙核3 GHz，內(nèi)存大小為2 G。從圖2可以看出在樣本數(shù)小于512時(shí)，ZrLZ(r=1/d)方法的運(yùn)算速度要明顯快于其它方法。實(shí)際中往往采用滑窗的方法來對實(shí)測圖像進(jìn)行估計(jì)，為了適應(yīng)不同的地物類型，每個(gè)滑窗包含的像素點(diǎn)一般都不會(huì)很多，因此該速度優(yōu)勢具有實(shí)用價(jià)值。

4.2 實(shí)測數(shù)據(jù)分析

傳統(tǒng)的PolSAR圖像統(tǒng)計(jì)模型的假設(shè)檢驗(yàn)方法是對各極化通道的幅度或強(qiáng)度數(shù)據(jù)分別進(jìn)行假設(shè)檢驗(yàn)。該方法的缺點(diǎn)在于只考慮了協(xié)方差矩陣的對角元素，沒有用到協(xié)方差矩陣中非對角元素的信息。最近，文獻(xiàn)[16]提出了基于矩陣對數(shù)累積量的假設(shè)檢驗(yàn)方法，實(shí)驗(yàn)結(jié)果表明該假設(shè)檢驗(yàn)方法有效并且直觀。對于模型參數(shù)未知需要從樣本值估計(jì)得到的復(fù)合假設(shè)檢驗(yàn)，統(tǒng)計(jì)量可由式(23)得到

假設(shè)檢驗(yàn)中概率值(probability value, p值)是指在由 H0所規(guī)定的總體中隨機(jī)抽樣，獲得等于及大于現(xiàn)有統(tǒng)計(jì)量的概率。根據(jù)式(23)所示的統(tǒng)計(jì)量可以通過蒙特卡羅仿真的方法得到p值[16]。進(jìn)行假設(shè)檢驗(yàn)時(shí)，向量k的維數(shù)要求大于m。對于K分布，其階數(shù)大于1的MLC與參數(shù)Σ無關(guān)。另外L采用名義視圖數(shù)代替，此時(shí)只剩下一個(gè)未知參數(shù)α，式(23)中取k=[,]T。

這里取的實(shí)測數(shù)據(jù)為PolSAR Pro軟件中常用的AIRSAR系統(tǒng)San Francisco地區(qū)的L波段PolSAR圖像(10 m×10 m)和 ESAR 系統(tǒng) Oberpfaffenhofen地區(qū)的L波段PolSAR圖像(3 m×3 m)，其中多視視圖數(shù)假定為 L = 4 (2×2)。在兩幅圖像中各選取了4個(gè)大小為20×20的區(qū)域進(jìn)行參數(shù)估計(jì)與假設(shè)檢驗(yàn)。

兩幅圖像中 4個(gè)區(qū)域的選擇及其（k2,k3）散點(diǎn)圖如圖3和圖4所示。該散點(diǎn)圖是在400個(gè)樣本中隨機(jī)選取200個(gè)計(jì)算（k2,k3），并且重復(fù)該過程50次得到的。圖中的兩條線分別為K分布和G0分布的理論（k2, k3）曲線[16]。San Francisco圖像中選取的4個(gè)區(qū)域分別為城區(qū)、海洋、植被A和植被B，由圖3可知城區(qū)的MLC靠近G0分布，海洋的MLC靠近Wishart分布，而植被A和植被B靠近K分布。Oberpfaffenhofen圖像中選取的4個(gè)區(qū)域分別為城區(qū)、森林、植被A和植被B，由圖4可知城區(qū)MLC在K分布和G0分布之間，實(shí)際上這是Fisher分布的區(qū)域[16]，而其它的3個(gè)區(qū)域都靠近K分布。

圖1 在不同的α取值下，不同參數(shù)估計(jì)方法在樣本數(shù)為512時(shí)的相對估計(jì)偏差與相對方差

圖2 本文方法與原有方法在不同樣本數(shù)下的運(yùn)算時(shí)間比較

圖3 San Francisco區(qū)域選取及MLC散點(diǎn)圖

圖4 Oberpfaffenhofen區(qū)域選取及MLC散點(diǎn)圖

由于 FMOM 估計(jì)性能較差，這里不做比較。由仿真比較可知 ZrLZ無解析解的估計(jì)方法中r=0.60時(shí)估計(jì)性能較差，r=0.05與r=0.20估計(jì)性能相差不大。因此這里參數(shù)估計(jì)方法選取ZrLZ(r=0.20和r=1/d), SMLC和SMOM方法。下面我們在顯著性水平α下，對各個(gè)區(qū)域在不同的參數(shù)估計(jì)方法下檢驗(yàn)假設(shè) H0：區(qū)域內(nèi)樣本值服從K(L = 4 ,Σ, α︿ )。

表1給出了San Francisco 4個(gè)區(qū)域不同算法下的K分布模型的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn)p值，從表中可以看出，對于城區(qū)，所有估計(jì)方法得到的 p值都無法通過顯著性水平α=0.05的假設(shè)檢驗(yàn)。從圖3中散點(diǎn)圖中可看到這是因?yàn)樵搮^(qū)域的樣本MLC更靠近G0分布，離K分布較遠(yuǎn)。對于海洋區(qū)域，所有的估計(jì)方法都能通過α=0.05的假設(shè)檢驗(yàn)，但是相比較來說SMLC方法的p值最大，SMOM方法p值最小，而其它兩種方法p值相近。在植被A區(qū)域，4種方法的p值都較大。而在植被B區(qū)域，SMOM方法的p值特別小，無法通過α=0.05的假設(shè)檢驗(yàn)，而其它3種方法的p值都很大。

表2給出了Oberpfaffenhofen 4個(gè)區(qū)域不同算法下的 K分布模型的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn) p值，從表中可以看出，對于城區(qū)，同樣地所有估計(jì)方法得到的p值都無法通過顯著性水平α=0.05的假設(shè)檢驗(yàn)。對于森林區(qū)域，所有估計(jì)方法的p值都較大，但是相比較來說SMOM方法最小，而其它3種方法結(jié)果相差不大。對于植被A區(qū)域，SMOM方法無法通過α=0.05的假設(shè)檢驗(yàn)，而其它 3種方法的p值都很大。對于植被B區(qū)域，SMOM方法無法通過α=0.05的假設(shè)檢驗(yàn)，而其它方法都能通過α=0.05的假設(shè)檢驗(yàn)。

表1 San Francisco 4個(gè)區(qū)域不同算法下的K分布模型的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn)p值

總結(jié)來說對于樣本分布本身就離K分布較遠(yuǎn)的區(qū)域，例如兩幅圖中的城區(qū)，在α=0.05下四種方法都拒絕假設(shè)0H。對于接近 K分布的實(shí)測數(shù)據(jù)區(qū)域，SMOM方法的穩(wěn)健性較差，其它3種方法的都能通過α=0.05的假設(shè)檢驗(yàn)，并且p值相差不大。但是另一方面從算法復(fù)雜度來說，ZrLZ(r=1/d)具有解析的表達(dá)式，其計(jì)算速度要明顯優(yōu)于其它方法(見圖 2)。

表2 Oberpfaffenhofen 4個(gè)區(qū)域不同算法下的K分布模型的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn)p值

5 結(jié)論

本文以|z |rlg|z|混合矩為基礎(chǔ)推導(dǎo)了K分布的參數(shù)估計(jì)方法，并且該方法在r=1/d時(shí)有解析表達(dá)式。研究了最優(yōu)r值的選擇問題，仿真結(jié)果表明在選定的幾個(gè)r值中，r取0.20時(shí)最優(yōu)。最后用仿真數(shù)據(jù)和實(shí)測數(shù)據(jù)對本文方法與原有方法進(jìn)行了比較，實(shí)驗(yàn)結(jié)果表明本文方法對不同α值的適應(yīng)性最好，并且ZrLZ方法在r=1/d時(shí)運(yùn)算速度的優(yōu)勢明顯。但是最優(yōu)r值的推導(dǎo)問題還有待研究。另外可以將基于|z |rlg|z|混合矩的參數(shù)估計(jì)方法的推導(dǎo)應(yīng)用于G0分布和U分布等PolSAR圖像模型，其推導(dǎo)結(jié)果及估計(jì)效果也值得研究。

[1] Collins M J, Denbina M, and Atteia G. On the reconstruction of quad-pol SAR data from compact polarimetry data for ocean target detection[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 591-600.

[2] Akbari V, Doulgeris A P, Moser G, et al.. A textural–contextual model for unsupervised segmentation of multipolarization synthetic aperture radar images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(4): 2442-2453.

[3] Lee J S, Schuler D, Lang R, et al.. K-distribution for multi-look processed polarimetric SAR imagery[C].Preceedings of the IEEE Geoscience and Remote Sensing Symposium, Pasadena, 1994: 2179-2181.

[4] Khan S and Guida R. On fractional moments of multilook polarimetric whitening filter for polarimetric SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing, 2014,52(6): 3502-3512.

[5] Cheng J, Gao G, Ding W, et al.. An improved scheme for parameter estimation of G0 distribution model in highresolution SAR images[J]. Progress in Electromagnetics Research, 2013, 134: 23-46.

[6] Sei T, Shibata H, Takemura A, et al.. Properties and applications of Fisher distribution on the rotation group[J].Journal of Multivariate Analysis, 2013, 116: 440-445.

[7] Liu T, Cui H G, Mao T, et al.. Modeling multilook polarimetric SAR images with heavy-tailed rayleigh distribution and novel estimation based on matrix log-cumulants[J]. Science China Information Sciences, 2013, 56(6): 1-14.

[8] Iskander D R and Zoubir A M. Estimating the parameters of K-distribution using higher-order and fractional moments[J].IEEE Transactions on Aerospace and Electronic Systems,1999, 35(4): 1453-1457.

[9] Blacknell D. Comparison of parameter estimators for K-distribution[J]. IEE Proceedings-Radar, Sonar and Navigation, 1994, 141(1): 45-52.

[10] Frery A C, Correia A H, and Da Freitas C. Classifying multifrequency fully polarimetric imagery with multiple sources of statistical evidence and contextual information[J].IEEE Transactions on Geoscience and Remote Sensing, 2007,45(10): 3098-3109.

[11] Blacknell D and Tough R. Parameter estimation for the K-distribution based on [zlog(z)][J]. IEE Proceedings-Radar,Sonar and Navigation, 2001, 148(6): 309-312.

[12] 胡文琳, 王永良, 王首勇. 基于 zrlog(z)期望的K分布參數(shù)估計(jì)[J]. 電子與信息學(xué)報(bào), 2008, 30(1): 203-205.Hu Wen-lin, Wang Yong-liang, and Wang Shou-yong.Estimation of the parameters of K-distribution based on zrlog(z)[J]. Jounal of Electronics & Information Technology,2008, 30(1): 203-205.

[13] Doulgeris A P, Anfinsen S N, and Eltoft T. Classification with a non-Gaussian model for PolSAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008,46(10): 2999-3009.

[14] Anfinsen S N, Doulgeris A P, and Eltoft T. Estimation of the equivalent number of looks in polarimetric synthetic aperture radar imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(11): 3795-3809.

[15] Anfinsen S N and Eltoft T. Application of the matrix-variate Mellin transform to analysis of polarimetric radar images[J].IEEE Transactions on Geoscience and Remote Sensing, 2011,49(6): 2281-2295.

[16] Anfinsen S N, Doulgeris A P, and Eltoft T. Goodness-of-fit tests for multilook polarimetric radar data based on the Mellin transform[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(7): 2764-2781