崔浩貴，晏 慶，張 曉，楊 飛，范路芳

(中國(guó)人民解放軍91001部隊(duì)，北京 100841)

0 引言

極化合成孔徑雷達(dá)(Polarimetric Synthetic Aperture Radar，PolSAR)雜波統(tǒng)計(jì)建模及其參數(shù)估計(jì)方法是PolSAR圖像相干斑抑制[1]、分類識(shí)別[2-3]和目標(biāo)檢測(cè)[4]等圖像解譯手段的重要課題。在PolSAR圖像中，均勻區(qū)域一般用高斯模型進(jìn)行建模，而非均勻區(qū)域常用多變量乘積模型來進(jìn)行非高斯建模。多變量乘積模型將PolSAR復(fù)散射向量表示為一個(gè)紋理分量和一個(gè)服從復(fù)高斯分布的相干斑分量的乘積，并且二者之間相互統(tǒng)計(jì)獨(dú)立。G0分布模型[5-6]是應(yīng)用非常廣泛的PolSAR圖像非高斯統(tǒng)計(jì)模型，它是在假設(shè)紋理分量服從逆Gamma分布時(shí)根據(jù)多變量乘積模型得到的。G0分布模型特別適用于城區(qū)等極不均勻區(qū)域的統(tǒng)計(jì)建模。

尋找快速、準(zhǔn)確的參數(shù)估計(jì)是G0分布模型研究的核心問題。根據(jù)采用的數(shù)據(jù)源，可將G0分布參數(shù)估計(jì)方法分為單極化估計(jì)方法和全極化估計(jì)方法。其中傳統(tǒng)的單極化估計(jì)方法分別估計(jì)出每個(gè)極化通道的參數(shù)，然后將各通道的參數(shù)求平均。而全極化估計(jì)方法利用了各通道之間的極化信息，其估計(jì)性能優(yōu)于單極化的方法。目前最常用的全極化估計(jì)方法是Anfinsen等提出的基于協(xié)方差矩陣的對(duì)數(shù)累積量(Matrix Log Cumulants，MLC)的參數(shù)估計(jì)方法，該方法本質(zhì)上是協(xié)方差矩陣行列式值的Mellin變換[7-11]。最近，Khan等人[12]將基于多視極化白化濾波器(Multilook Polarimetric Whitening Filter，MPWF)的參數(shù)估計(jì)方法擴(kuò)展到分?jǐn)?shù)階矩(Fractional Moments，F(xiàn)M)的情形，提出了基于多視極化白化濾波器分?jǐn)?shù)階矩(MPWF-FM)的參數(shù)估計(jì)方法。實(shí)驗(yàn)結(jié)果表明，該方法比MLC方法具有更小估計(jì)方差和均方誤差(Mean Square Error，MSE),但是該方法存在計(jì)算復(fù)雜、最優(yōu)的分?jǐn)?shù)階矩?zé)o法理論推導(dǎo)得到和當(dāng)紋理參數(shù)較小時(shí)失效等問題。

本文提出了一種基于多視極化白化濾波器對(duì)數(shù)累積量(Log Cumulants of MPWF，MPWF-LC)的參數(shù)估計(jì)新方法。新方法解決了MPWF-FM方法存在的估計(jì)失效等問題，并且具有更優(yōu)的估計(jì)性能。通過G0分布數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)分別進(jìn)行了仿真分析，結(jié)果有效驗(yàn)證了新方法的參數(shù)估計(jì)結(jié)果準(zhǔn)確度高、計(jì)算速度快。

1 PolSAR圖像乘積模型

PolSAR系統(tǒng)在測(cè)量中獲取的信息可用復(fù)散射向量k表示：

k=[Shh,Shv,Svh,Svv]T，

(1)

式中，[·]T表示向量的轉(zhuǎn)置；Sxy表示發(fā)射極化方式為x，接收極化方式為y的復(fù)散射系數(shù)；h表示水平極化；v表示垂直極化。

多變量乘積模型常用來對(duì)PolSAR圖像進(jìn)行非高斯建模，它將復(fù)散射向量表示為相互統(tǒng)計(jì)獨(dú)立的紋理分量τ和相干斑分量y的乘積，即：

(2)

式中，y服從復(fù)高斯分布，其概率密度函數(shù)為：

(3)

式中，d是復(fù)高斯矢量的維數(shù)；Γ=E{yyH}，E{·}表示隨機(jī)變量的數(shù)學(xué)期望。

實(shí)際中，為了抑制相干斑，常對(duì)圖像進(jìn)行多視處理，其過程可表示為：

(4)

式中，L為視圖數(shù)，L≥d；上標(biāo)“H”表示共軛轉(zhuǎn)置；l代表第l個(gè)像素點(diǎn)。假設(shè)進(jìn)行多視處理的像素點(diǎn)具有相同紋理分量，此時(shí)式(4)可簡(jiǎn)化為：

C=τY，

(5)

其中，

(6)

式中，Y服從Wishart分布，其概率密度函數(shù)為：

(7)

式中，d是Y的維數(shù)；Γ=E(Y)；Tr(·)表示矩陣的跡；函數(shù)Γd(L)為復(fù)數(shù)形式的多變量伽馬函數(shù)：

(8)

式中，Γ(·)為標(biāo)準(zhǔn)Eular伽馬函數(shù)。由于這里假設(shè)在多視窗口內(nèi)紋理分量恒定，樣本的協(xié)方差矩陣可表示為：

Σ=E{C}=E{τ}E{Y}=E{τ}Γ。

(9)

當(dāng)式中紋理分量服τ為逆伽馬分布，即τ～γ-1(u,α)時(shí)，其乘積模型服從G0分布。

2 基于MPWF對(duì)數(shù)累積量的參數(shù)估計(jì)方法

MPWF濾波器通過最優(yōu)地組合極化協(xié)方差矩陣中的所有元素，得到一幅降斑的圖像[13-14]。進(jìn)行MPWF濾波后的數(shù)據(jù)可表示為：

(10)

將式(5)和式(9)代入式(10)，MPWF濾波后的數(shù)據(jù)可分解為：

(11)

x～γ(d,Ld)，

(12)

其中，伽馬分布γ(u,d)的概率密度分布函數(shù)為：

(13)

基于Mellin變換的對(duì)數(shù)累積量定義為：

(14)

式中，kr{z}表示隨機(jī)變量z的第r階對(duì)數(shù)累積量；E{zs-1}為關(guān)于隨機(jī)變量z的數(shù)學(xué)期望。當(dāng)隨機(jī)變量服從如式(11)所示的乘積模型時(shí)，其對(duì)應(yīng)的對(duì)數(shù)累積量為加性模型[8]：

(15)

由于隨機(jī)變量x服從如式(12)所示的伽馬分布，可得其對(duì)數(shù)累積量為：

(16)

式中，ψ(x)=Γ′(x)/Γ(x)為digamma函數(shù)；Γ(x)為伽馬函數(shù)。ψ(m,x)表示m階Polygamma函數(shù)：

(17)

(18)

將式(16)、式(18)代入式(15)，可得基于MPWF對(duì)數(shù)累積量(MPWF-LC)的G0分布參數(shù)估計(jì)方法為：

(19)

(20)

3 仿真結(jié)果與分析

用G0分布仿真數(shù)據(jù)比較以下幾種方法的相對(duì)估計(jì)偏差、相對(duì)估計(jì)方差、相對(duì)MSE和計(jì)算時(shí)間：

① MPWF-FM[12]；

② MLC方法，取二階MLC進(jìn)行參數(shù)估計(jì)[8]；

③ 本文提出的MPWF-LC參數(shù)估計(jì)方法，取二階對(duì)數(shù)累積量的表達(dá)式，即在式中取r=2。

(a)相對(duì)估計(jì)偏差

圖1 G0分布情形，不同參數(shù)λ下3種估計(jì)方法的相對(duì)估計(jì)偏差、相對(duì)估計(jì)方差與相對(duì)MSEFig.1 The relative estimation bias, relative estimation variance and relative MSE for three estimation methods with different parameters of λ in the case of G0 distribution

G0分布下各估計(jì)方法在不同樣本數(shù)下的計(jì)算時(shí)間如圖2所示。

圖2 各估計(jì)方法的計(jì)算時(shí)間比較Fig.2 Comparison of calculation time of different estimation methods

仿真數(shù)據(jù)參數(shù)L=10，λ=10。采用Matlab進(jìn)行編程計(jì)算，計(jì)算機(jī)CPU為Intel E5700，雙核3 GHz，內(nèi)存大小2 GB。由圖2可以看出，當(dāng)樣本數(shù)較小時(shí)，MPWF-LC計(jì)算所花的時(shí)間最短，計(jì)算復(fù)雜度最低；隨著樣本數(shù)的增加，當(dāng)樣本大于100后，MLC具有最短的計(jì)算時(shí)間；另外，在任意樣本數(shù)下，MPWF-FM所需的計(jì)算時(shí)間最長(zhǎng)。

4 實(shí)測(cè)數(shù)據(jù)分析

用實(shí)測(cè)數(shù)據(jù)對(duì)MPWF-LC，MPWF-FM和MLC三種方法進(jìn)行分析。選取的實(shí)測(cè)數(shù)據(jù)為AIRSAR系統(tǒng)于1989年獲得的Flevoland地區(qū)L波段PolSAR圖像(12.1 m×6.7 m)，以及1988年獲得的San Francisco地區(qū)的L波段PolSAR圖像(10 m×10 m)。對(duì)原始數(shù)據(jù)進(jìn)行了L=4(2×2)的多視處理，由于Flevoland地區(qū)的數(shù)據(jù)源為4視數(shù)據(jù)，因此其視圖數(shù)變?yōu)長(zhǎng)=16。

采用基于對(duì)數(shù)累積量的假設(shè)檢驗(yàn)方法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行擬合分析。簡(jiǎn)單假設(shè)檢驗(yàn)情形下，統(tǒng)計(jì)量Qp可表示為：

Qp=n(〈k〉-k)TK-1(〈k〉-k)，

(21)

式中，〈k〉為樣本的MLC，k為特定模型下MLC的理論值；n為樣本數(shù)；下標(biāo)p為〈k〉或k的維數(shù)，K=nE{(〈k〉-k)(〈k〉-k)T}。可知Qp服從自由度為p的χ2分布[15]：

Qp→χ2(p)。

(22)

假設(shè)檢驗(yàn)中的概率值(Probability Value，p值)是指在由H0所規(guī)定的總體中隨機(jī)抽樣，獲得等于及大于現(xiàn)有統(tǒng)計(jì)量的概率。由上述分析可知，p值可由統(tǒng)計(jì)量Qp根據(jù)χ2(p)分布計(jì)算得到。在參數(shù)估計(jì)中視圖數(shù)L用名義視圖數(shù)代替，向量〈k〉和k采用2階和3階MLC，即〈k〉=[〈k2〉,〈k3〉]T，k=[k2,k3]T。

在2幅圖像中各選取了4個(gè)大小為20 pixel×20 pixel的區(qū)域進(jìn)行參數(shù)估計(jì)與假設(shè)檢驗(yàn)，比較各參數(shù)估計(jì)方法在G0分布假設(shè)下的p值。圖3和圖4分別給出了2幅圖像中4個(gè)區(qū)域的選擇以及各區(qū)域的樣本對(duì)數(shù)累積量〈k〉。如圖3(a)所示，F(xiàn)levoland圖像中選取的4個(gè)區(qū)域分別為河流、森林、小麥和油菜籽區(qū)域[16]。由圖3(b)可知，河流的MLC靠近G0分布，森林的MLC靠近K分布，而小麥和油菜籽區(qū)域與Wishart、K和G0分布都非常接近。如圖4(a)所示，San Francisco圖像中選取的4個(gè)區(qū)域分別為海洋、城區(qū)、植被1和植被2。由圖4(b)可知，海洋的MLC靠近Wishart分布，城區(qū)的MLC靠近G0分布，而植被1和植被2靠近K分布。

(a)4個(gè)區(qū)域的示意(Pauli RGB圖像)

(a)4個(gè)區(qū)域的示意(Pauli RGB圖像)

Flevoland四個(gè)區(qū)域在G0分布模型下的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn)p值如表1所示。由表1可以看出，對(duì)于河流區(qū)域，所有估計(jì)方法得到的p值都能通過顯著性水平α=0.05的假設(shè)檢驗(yàn)，但是MLC方法的p值最大。對(duì)于小麥和油菜籽區(qū)域，所有估計(jì)方法的p值都較大。森林區(qū)域同樣地?zé)o法通過顯著性水平α=0.05的假設(shè)檢驗(yàn)。

表1 Flevoland四個(gè)區(qū)域在G0分布模型下的參數(shù)估計(jì)結(jié)果及p值

San Francisco四個(gè)區(qū)域在G0分布模型下的參數(shù)估計(jì)結(jié)果及假設(shè)檢驗(yàn)p值如表2所示。由表2可以看出，對(duì)于城區(qū)，只有MLC方法的p值都能通過顯著性水平α=0.05的假設(shè)檢驗(yàn)，其他2種方法p值很小。對(duì)于海洋區(qū)域，3種方法的結(jié)果都能通過顯著性水平α=0.05的假設(shè)檢驗(yàn)，但是MLC方法的p值略小。植被1和植被2區(qū)域，所有估計(jì)方法的p值都無法通過顯著性水平α=0.05的假設(shè)檢驗(yàn)。

表2 San Francisco四個(gè)區(qū)域在G0分布模型下的參數(shù)估計(jì)結(jié)果及p值

從上述分析可以看出，MPWF-LC和MPWF-FM方法的參數(shù)估計(jì)結(jié)果和p值都非常接近。而MLC方法在某些區(qū)域表現(xiàn)較好，例如Flevoland圖像中的河流和San Francisco圖像中的植被2區(qū)域，這是因?yàn)樵搮?shù)估計(jì)方法和假設(shè)檢驗(yàn)都是基于MLC。在其他區(qū)域，MLC方法和另外2種方法的結(jié)果也基本上相差不大。從該結(jié)果可以看出，本文提出的MPWF-LC方法在實(shí)測(cè)數(shù)據(jù)參數(shù)估計(jì)中的有效性。

5 結(jié)束語

以多視極化白化濾波器和對(duì)數(shù)累積量為基礎(chǔ)，提出了基于多視極化白化濾波器對(duì)數(shù)累積量的G0分布模型的參數(shù)估計(jì)方法。用仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)對(duì)本文方法與MLC方法以及MPWF-FM方法進(jìn)行了比較。實(shí)驗(yàn)結(jié)果表明，本文方法具有最小的估計(jì)方差和MSE，并且在樣本數(shù)較少時(shí)，計(jì)算時(shí)間最短。另外，可以將本文方法應(yīng)用于Fisher分布、U分布等PolSAR圖像模型，其推導(dǎo)結(jié)果及估計(jì)效果值得研究。