極化SAR圖像K分布模型參數(shù)估計(jì)方法研究與仿真

晏慶,崔浩貴,易帆帆

(1.91469部隊(duì),北京 100841;2.海軍工程大學(xué) 電子工程學(xué)院,武漢 430033)

極化SAR圖像K分布模型參數(shù)估計(jì)方法研究與仿真

晏慶1,崔浩貴1,易帆帆2

(1.91469部隊(duì),北京 100841;2.海軍工程大學(xué) 電子工程學(xué)院,武漢 430033)

摘要：極化合成孔徑雷達(dá)(SAR)已成為一種不可或缺的地理遙感和軍事偵察信息源,但是其圖像解譯困難,需要借助精確的建模方法與分析手段。K分布模型是常用的極化SAR圖像模型,但是傳統(tǒng)參數(shù)估計(jì)方法準(zhǔn)確性較差的問題影響了模型的擬合精度．針對該問題,首先推導(dǎo)了K分布矩陣行列式值的概率密度函數(shù)及其最大似然(ML)估計(jì)方法,并用仿真數(shù)據(jù)驗(yàn)證了其有效性，然后將ML估計(jì)方法應(yīng)用于K分布模型的等效視圖數(shù)(ENL)估計(jì)中,比較了其與傳統(tǒng)的方法在ENL估計(jì)上的性能差異。仿真和實(shí)測數(shù)據(jù)的實(shí)驗(yàn)結(jié)果均表明，新方法具有更小的估計(jì)偏差和方差,對解決極化SAR圖像中K分布區(qū)域的參數(shù)準(zhǔn)確估計(jì)問題有重要的實(shí)用價(jià)值。

關(guān)鍵詞：K分布;雷達(dá)圖像;最大似然估計(jì);協(xié)方差矩陣

0引言

極化合成孔徑雷達(dá)(SAR)具有全天時(shí)、全天候的多維信息遙感能力和不同目標(biāo)極化特性的鑒別能力,但是相干成像技術(shù)都固有的相干斑增加了極化SAR圖像解譯的難度。為了從極化SAR圖像中提取出有用的信息,建立精確的相干斑模型并進(jìn)行準(zhǔn)確的數(shù)學(xué)分析是非常有效的手段。傳統(tǒng)的相干斑建模方法是假設(shè)散射回波服從高斯分布[1]。

但是隨著雷達(dá)分辨率的提高,需要借助非高斯模型來刻畫回波的特性。在極化SAR圖像非高斯建模領(lǐng)域中,最常用的是基于乘積模型推導(dǎo)得到的K分布[2]，該模型將協(xié)方差矩陣表示為服從伽馬分布的紋理變量與服從Wishart分布的相干斑變量的乘積。

參數(shù)估計(jì)的準(zhǔn)確性很大程度上決定了極化SAR圖像K分布模型的擬合精度。最常用的矩陣估計(jì)方法存在表達(dá)式非常復(fù)雜并且估計(jì)精度不高的問題[3]。Anfinsen等將Mellin變換應(yīng)用到了K分布的參數(shù)估計(jì)中,提出了基于一階對數(shù)累積量的等效視圖數(shù)(ENL)估計(jì)方法[3-5]。但是該方法在預(yù)先估計(jì)矩陣行列式值時(shí)采用了先求協(xié)方差矩陣樣本的均值,再取其行列式值的方法，這實(shí)際上是一種有偏估計(jì)方法,導(dǎo)致ENL的估計(jì)性能較差[6]。

本文推導(dǎo)了K分布矩陣行列式值的分布和其最大似然(ML)估計(jì)方法,用查表法解決了ML估計(jì)方法計(jì)算復(fù)雜度高的問題。仿真結(jié)果表明,相比傳統(tǒng)的估計(jì)方法,本文方法具有更小的偏差和方差。新方法為ENL這一關(guān)鍵參數(shù)的準(zhǔn)確估計(jì)提供了新的途徑,對后續(xù)的PolSAR圖像相干斑濾波和檢測等領(lǐng)域的應(yīng)用有著重要意義。

1極化SAR圖像協(xié)方差矩陣建模

1.1極化SAR圖像特性

極化SAR是一種主動(dòng)式的微波成像設(shè)備,具有分辨率不受雷達(dá)傳感器高度影響的優(yōu)點(diǎn),其雷達(dá)信號(hào)能穿透云層。另外,波長較長的L波段信號(hào)能穿透森林和地表植被覆蓋,在軍事上能夠用于發(fā)現(xiàn)叢林中的隱藏目標(biāo)和獲取地表以下的信息,這些是光學(xué)等其它傳感器無法實(shí)現(xiàn)的。

最常見的極化SAR為機(jī)載形式,美國NASA的JPL實(shí)驗(yàn)室研制的AIRSAR系統(tǒng)可以在P、L和C三個(gè)波段上進(jìn)行全極化測量成像。該系統(tǒng)自首飛后每年至少執(zhí)行一次測量任務(wù),成為驗(yàn)證新開發(fā)雷達(dá)技術(shù)的主要實(shí)測數(shù)據(jù)來源。AIRSAR系統(tǒng)安裝于DC-8飛機(jī)的底部,如圖1所示。

圖1　AIRSAR系統(tǒng)

Pauli分解是將目標(biāo)的散射矩陣表示為各Pauli矩陣的加權(quán)和,其中每個(gè)Pauli矩陣都對應(yīng)著一種特定的散射機(jī)制。在得到Pauli分解的RGB合成圖時(shí)往往只用其中三種散射機(jī)制進(jìn)行合成,分別為奇次散射、偶次散射和體散射。其中奇次反射一般發(fā)生在平面上,例如海洋或者光禿的地面;偶次散射經(jīng)常發(fā)生于城區(qū);體散射常見于植被區(qū)域。在進(jìn)行合成時(shí)這三種散射機(jī)制分別用藍(lán)色、紅色和綠色來表示,其加權(quán)值代表顏色的深度。

圖2示出了Flevoland地區(qū)和San Francisco海灣地區(qū)極化SAR圖像的Pauli分解RGB合成圖。這兩幅圖像都是由AIRSAR系統(tǒng)獲取的,其圖像大小分別為750×1024和900×1024,是常用的極化SAR實(shí)測數(shù)據(jù)。

圖2　極化SAR圖像RGB合成圖Pauli分解　　　(a)Flevoland； (b)San Francisco

1.2極化SAR圖像K分布建模

為了解譯如圖2所示的極化SAR圖像,首先要對其進(jìn)行建模。常用的非高斯建模方法是乘積模型,在該模型假設(shè)下,極化SAR圖像的協(xié)方差矩陣C可表示為

C=τY,

(1)

式中: τ為紋理變量; Y為相干斑變量,并且紋理變量和相干斑變量相互統(tǒng)計(jì)獨(dú)立。一般認(rèn)為相干斑變量Y服從Wishart分布,其PDF為[5]

(2)

式中:L為視圖數(shù);d表示矩陣的維數(shù);Γ為方差矩陣; Tr(·)表示取矩陣的跡;Γd(L)為多變量伽馬函數(shù)

(3)

式中,Γ(L)為標(biāo)準(zhǔn)歐拉函數(shù)。

當(dāng)紋理變量τ滿足伽馬分布時(shí),根據(jù)式(1)所示的乘積模型,可以得到極化SAR圖像的協(xié)方差矩陣服從K分布。一般假設(shè)紋理變量滿足單位均值的條件,歸一化伽馬分布的概率密度函數(shù)(PDF)為

(4)

式中,α為伽馬分布的形狀參數(shù)。

根據(jù)式(1)、式(2)和式(4),可以得到K分布的PDF為

(5)

式中,Σ=E(C);Kα-Ld(·)為第α-Ld階的第二類修正Bessel函數(shù)。

2K分布矩陣行列式值的最大似然估計(jì)

特征函數(shù)(CF)是研究隨機(jī)變量PDF的一個(gè)重要工具,有些隨機(jī)變量要確定其分布可能比較困難,但是確定其特征函數(shù)卻比較簡單。平常所指的特征函數(shù)實(shí)際上是隨機(jī)變量PDF的逆Fourier變換。文獻(xiàn)[7]中Nicolas提出了用Mellin變換代替Fourier變換來分析隨機(jī)變量的分布,對于隨機(jī)變量X∈R+,定義其基于Mellin變換的特征函數(shù)(MCF)為[8]

(6)

式中,p(x)為隨機(jī)變量X的PDF,對應(yīng)的逆變換為

p(x)=M-1{φX(s)}(x)

(7)

對式(1)所示的乘積模型進(jìn)行Mellin變換,可得其MCF存在如下關(guān)系:

φC(s)=φτ(s)·φY(s).

(8)

由于x=(2L)d|C|/|Σ|服從d個(gè)χ2分布的乘積[9],可得其Mellin變換為

(9)

利用式(7)的定義對式(9)進(jìn)行Mellin逆變換,令Υ=|C|,可以求出協(xié)方差矩陣行列式值分布的PDF

(10)

式中MeijerG函數(shù)的定義為

(11)

那么其最大似然估計(jì)可由對數(shù)似然比公式求取

(12)

求得

(13)

其中參數(shù)z為下式的解

(14)

3仿真實(shí)驗(yàn)分析

3.1K分布矩陣行列式值ML估計(jì)的性能分析

在實(shí)際中利用式(13)對行列式值進(jìn)行最大似然估計(jì)時(shí),需要首先求解式(14)?？梢钥吹揭环矫媸?14)沒有解析的表達(dá)式,需要用搜索的方法求解。另一方面,該式中包含復(fù)雜的MeijerG函數(shù)。如在用樣本數(shù)據(jù)對行列式值進(jìn)行估計(jì)的過程中實(shí)時(shí)求解式(14)的話會(huì)碰到計(jì)算復(fù)雜度較大的問題,這降低了算法的實(shí)用性。這里采用查表法來解決ML方法計(jì)算復(fù)雜的問題。注意到式(13)中的d實(shí)際上為常數(shù),即d=3.令

f(L;α)=(Lα/d)d，

(15)

則式(13)可表示為

(16)

其中f(L;α)定義為行列式值最大似然修正系數(shù)。

在進(jìn)行參數(shù)估計(jì)之前,首先仿真得到二維表格形式的最大似然修正系數(shù)f(L;α)。由于α等于200時(shí),K分布已經(jīng)近似于高斯情形,而L要求大于d.因此在建表時(shí),參數(shù)選擇的范圍為α∈(0,200],L∈(3,1000]。圖3示出了K分布的形狀參數(shù)α分別為2,4,200時(shí),仿真得到的最大似然修正系數(shù)關(guān)于L的變化曲線。從圖中可以看出,隨著L的增加,最大似然修正系數(shù)趨向于一個(gè)固定值。其次,可以看出不同形狀參數(shù)下的修正系數(shù)只有在L較小時(shí)有稍微差別,隨著L的增加該差別基本可以忽略。另外,可以看到當(dāng)L<4時(shí),最大似然修正系數(shù)的變量速率較快,此時(shí)應(yīng)在建表時(shí)應(yīng)采用較小的間隔。因此在實(shí)際應(yīng)用中預(yù)先建立的是一個(gè)非均勻間隔索引的表:當(dāng)L≤4時(shí),以間隔值ΔL=0.01建表;當(dāng)L>4時(shí),以ΔL=0.1建表。

圖3　最大似然修正系數(shù)仿真圖

為了驗(yàn)證本文提出的估計(jì)方法的有效性,這里用K分布的仿真數(shù)據(jù)將本文方法與傳統(tǒng)的基于協(xié)方差矩陣均值的行列式值的估計(jì)方法(即|Σ|=|E{C}|)的結(jié)果進(jìn)行仿真比較。K分布的仿真數(shù)據(jù)采用Bootstrap的方法由式 (1)所示的乘積模型仿真得到。為了突出顯示數(shù)據(jù)分布圖的峰值,這里引入了一種基于核密度估計(jì)KDE的非參數(shù)估計(jì)方法[4]。其中h稱為帶寬,它決定圖像的平滑程度,在此采用Epanechnikov核函數(shù)。帶寬h的選擇對峰值的幅度有較大影響,但是對峰值的位置影響不大,本文選取帶寬h=0.3.

圖4示出了K分布仿真數(shù)據(jù)下,取不同的形狀參數(shù)時(shí),傳統(tǒng)的基于均值的行列式值估計(jì)方法和本文提出的ML估計(jì)方法的估計(jì)偏差和方差。這里仿真數(shù)據(jù)服從L=4,|Σ|=1的K分布,仿真次數(shù)M=10000．圖中的橫軸為樣本數(shù),圖4(a)的形狀參數(shù)為α=8,圖4(b)為α=2．從圖中首先可以看出隨著樣本數(shù)的增加,估計(jì)偏差和方差都在減小。其次,可以看出相比傳統(tǒng)的方法,本文方法具有更小的估計(jì)偏差和方差,其性能更優(yōu)。另外,對比圖4(a)和圖4(b),可以看出當(dāng)形狀參數(shù)較大時(shí),參數(shù)估計(jì)的偏差和方差相對較小，這是因?yàn)殡S著形狀參數(shù)的增大,K分布數(shù)據(jù)會(huì)變均勻,并趨向于高斯情形。而形狀參數(shù)較小時(shí),K分布的數(shù)據(jù)比較奇異,導(dǎo)致參數(shù)估計(jì)性能變差。

圖4　不同算法下K分布矩陣行列式值的估計(jì)偏差與方差 (a) α=8 ; (b) α=2

3.2矩陣行列式值ML估計(jì)方法在極化SAR圖像ENL估計(jì)中的應(yīng)用

ENL是極化SAR圖像中的一個(gè)重要參數(shù),它表征了參與多視的樣本的數(shù)目,一般采用基于一階對數(shù)累積量的方法進(jìn)行估計(jì)[3]

(17)

這里不詳細(xì)討論α的估計(jì),用基于二階矩特征的方法預(yù)先估計(jì)[3],并在處理過程中將其視為已知。

由于傳統(tǒng)的基于矩陣均值的行列式值的估計(jì)方法是有偏估計(jì),為了提高ENL估計(jì)的準(zhǔn)確性,我們可以聯(lián)合式(13)與式(17)來估計(jì)ENL。首先用名義視圖數(shù)代入式(13),估計(jì)出|Σ|;然后將|Σ|的估計(jì)值代入式(17)估計(jì)出ENL。該方法的步驟可總結(jié)如下:

1) 根據(jù)式(13)和(14),用數(shù)值仿真方法計(jì)算出最大似然修正系數(shù);

2) 預(yù)先建立的是一個(gè)非均勻間隔索引的表,即L和α關(guān)于f(L;α)的二維索引表:當(dāng)L≤4時(shí),以間隔值ΔL=0.01建表;當(dāng)L>4時(shí),以ΔL=0.1建表;

圖5示出了仿真數(shù)據(jù)情形下,分別用本文ML方法和傳統(tǒng)的基于均值的方法估計(jì)出行列式值后,由一階對數(shù)累積量方法估計(jì)得到的ENL結(jié)果。其中仿真數(shù)據(jù)服從參數(shù)為L=4,|Σ|=1的K分布,仿真次數(shù)為M=10000.圖5(a)的形狀參數(shù)為α=8,圖5(b)為α=2.從圖5可以看出隨著樣本數(shù)的增加,兩種方法的估計(jì)偏差和方差之間的差異在減小，其次,本文的ML估計(jì)方法在對K分布矩陣的行列式值進(jìn)行修正后,相比原有算法能有效改善ENL的估計(jì)精度。最后,對比圖5(a)和圖5(b)可以看到當(dāng)形狀參數(shù)較大時(shí),ENL估計(jì)結(jié)果的偏差和方差都相對較小。

圖5　不同算法下L的估計(jì)偏差與方差(a) α=8; (b) α=2

3.3實(shí)測數(shù)據(jù)仿真分析

用圖2所示的實(shí)測數(shù)據(jù)來對兩種ENL估計(jì)方法進(jìn)行仿真分析。首先從圖2(a)中可以看出,Flevoland圖像中大部分區(qū)域?yàn)檗r(nóng)田和植被等符合K分布的數(shù)據(jù),右上角的河流可認(rèn)為服從形狀參數(shù)較大的K分布。其次從圖2(b)可以看出,SanFrancisco主要由海洋、城區(qū)和植被3個(gè)部分構(gòu)成。其中,植被區(qū)域一般認(rèn)為服從K分布。城區(qū)較為奇異,可能滿足形狀參數(shù)較小的K分布。而海洋接近高斯分布,可認(rèn)為是形狀參數(shù)較大的K分布。

這里用無監(jiān)督的滑窗估計(jì)方法來對實(shí)測數(shù)據(jù)的ENL值進(jìn)行估計(jì)[4]。該方法首先用k×k的滑動(dòng)窗口將整個(gè)圖像分解成很多個(gè)小區(qū)域,然后根據(jù)每個(gè)區(qū)域中的樣本數(shù)據(jù)計(jì)算出其對應(yīng)的ENL值,最后對所有區(qū)域的ENL估計(jì)結(jié)果進(jìn)行統(tǒng)計(jì)。由于上述兩幅極化SAR圖像中的大部分區(qū)域都可以用K分布來進(jìn)行建模,這里取ENL統(tǒng)計(jì)分布圖的峰值作為整個(gè)圖像的ENL值。

圖6示出了實(shí)測數(shù)據(jù)下本文給出的ML估計(jì)方法和傳統(tǒng)的估計(jì)方法在不同的滑窗值k下ENL的估計(jì)結(jié)果。其中滑窗值的大小分別選取為k={2,3,5,7,11}。對比不同窗口k下的結(jié)果可以看出,當(dāng)窗口值較小時(shí)本文方法的優(yōu)勢比較明顯,但是隨著窗口值的增大兩種方法的差別在減小。該結(jié)論也與圖3和圖4中通過仿真數(shù)據(jù)得到的結(jié)果相符。另外,可以看到本文的ML估計(jì)方法的在不同窗口值下的表現(xiàn)比較穩(wěn)定,而傳統(tǒng)的方法在k較小時(shí)其分布的峰值出現(xiàn)偏差。

圖6　實(shí)測數(shù)據(jù)下各估計(jì)方法對ENL的估計(jì)結(jié)果的比較(a) Flevoland; (b) San Francisco

4結(jié)束語

本文研究了在三維情形下,K分布矩陣的行列式值的統(tǒng)計(jì)分布特性,推導(dǎo)了其PDF和ML估計(jì)方法。針對ML估計(jì)方法的表達(dá)式計(jì)算復(fù)雜度高的問題,提出了先仿真出ML修正系數(shù)再進(jìn)行查表的解決方法。仿真數(shù)據(jù)的實(shí)現(xiàn)表明,相比原有的基于矩陣均值的行列式值的估計(jì)方法,該方法的優(yōu)勢明顯。然后,將矩陣行列式值的ML估計(jì)應(yīng)用于極化SAR圖像的基于一階對數(shù)累積量的ENL估計(jì)方法中。最后用仿真數(shù)據(jù)和實(shí)測數(shù)據(jù)驗(yàn)證了該方法在對ENL估計(jì)的準(zhǔn)確性和有效性,實(shí)現(xiàn)解決現(xiàn)實(shí)其估計(jì)效果改善明顯。

但是,在用查表法獲得最大似然修正系數(shù)的過程中,由于表的間隔以及數(shù)據(jù)位數(shù)的限制,會(huì)使估計(jì)結(jié)果存在一定的誤差,該誤差的影響還有待分析。另外,四維或者更多維情形K矩陣行列式值的PDF及其ML估計(jì)方法也有待后續(xù)的研究。

參考文獻(xiàn)

[1]LEEJS,HOPPELKW,MANGOSA, et al.IntensityandphasestatisticsofmultilookpolarimetricandinterferometricSARimagery[J].IEEETransactionsonGeoscienceandRemoteSensing, 1994, 32(5):1017-1028.

[2]AKBARIV,DOULGERISAMOSERG, et al.ATextural-contextualmodelforunsupervisedsegmentationofmultipolarizationsyntheticapertureradarimages[J].IEEETransactionsonGeoscienceandRemoteSensing, 2013, 51(4):2442-2453.

[3]ANFINSENSN,ELTOFTT.Applicationofthematrix-variateMellintransformtoanalysisofpolarimetricradarimages[J].IEEETransactionsonGeoscienceandRemoteSensing, 2011, 49(6):2281-2295.

[4]ANFINSENSN,DOULGERISAP,ELTOFTT.Estimationoftheequivalentnumberoflooksinpolarimetricsyntheticapertureradarimagery[J].IEEETransactionsonGeoscienceandRemoteSensing, 2009, 47(11):3795-3809.

[5]ANFINSENSN,DOULGERISAP.OULGERISP,etal.Goodness-of-fittestsformultilookpolarimetricradardatabasedontheMellintransform[J].IEEETransactionsonGeoscienceandRemoteSensing, 2011, 49(7):2764-2781.

[6]劉濤, 崔浩貴, 高俊.Wishart分布矩陣行列式值的統(tǒng)計(jì)特性及其在參數(shù)估計(jì)中的應(yīng)用[J]. 電子學(xué)報(bào), 2013, 41(6):1231-1237.

[7]NICOLASJM.Introductiontosecondkindstatistics:applicationoflog-momentsandlog-cumulantstoSARimagelawanalysis[J].TraitementduSignal, 2002, 19(3):139-168.

[8]POULARIKASAD.Transformsandapplicationshandbook[D].CRCPress, 2010.

[9]GOODMANN.StatisticalanalysisbasedonacertainmultivariatecomplexGaussiandistribution(anintroduction)[J].theAnnalsofMathematicalStatistics, 1963, 34(1):152-177.

晏慶(1979-),男,湖南人,碩士， 主要研究方向?yàn)橥ㄐ偶夹g(shù)與計(jì)算機(jī)應(yīng)用。

崔浩貴(1987-),男,浙江人,博士,研究方向?yàn)槔走_(dá)極化信息處理、電子戰(zhàn)建模與仿真。

ResearchandSimulationonParameterEstimationMethodofKDistributionforPolarimetricSARImagery

YANQing1,CUIHaogui1,YIFanfan2

(1.Unit 91469, Beijing 100841, China;2.School of Electronic Engineering,Naval University of Engineering, Wuhan 430033, China)

Abstract:Thepolarimetricsyntheticapertureradar(SAR)isakindofindispensablemeansofremotesensingandmilitaryreconnaissance,buttheimagesneedaprecisemodelingandanalysistointerpret.TheKdistributionisusuallyusedtomodelpolarimetricsyntheticapertureradar(SAR)imageinthenon-Gaussiancase,butthebiasestimationoftheparametersaffectthefittingprecisionofthemodel.Inthispaper,theprobabilitydensityfunctionandthemaximumlikelihood(ML)estimationforthedeterminantofKdistributionmatrixarederived.Itsvalidityisalsoconfirmedbythesimulateddata.Then,theMLmethodisappliedtotheestimationofequivalentnumberoflooks(ENL)inPolSARimages.ComparesaremadebetweentheMLmethodandtheformermethodusingthedeterminantofthemeanofmatrix.Theresultsofexperimentsbasedonsimulatedandrealdatashowthattheproposedalgorithmobviouslypromotestheestimationbiasandvariance,whichisofpracticalvalueinsolvingtheaccurateestimationproblemforKdistributionregionsinpolarimetricSARimagery.

Keywords:K distribution; radar image; maximum likelihood estimation; covariance matrix

doi:10.13442/j.gnss.1008-9268.2016.02.013

收稿日期：2016-03-12

中圖分類號(hào)：TN95

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1008-9268(2016)02-0071-06

作者簡介

資助項(xiàng)目： 國家自然科學(xué)基金(批準(zhǔn)號(hào):61372165)

聯(lián)系人： 晏慶E-mail: 63934470@qq.com