極化SAR圖像K分布模型參數估計方法研究與仿真

晏慶,崔浩貴,易帆帆

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

極化SAR圖像K分布模型參數估計方法研究與仿真

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

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

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

關鍵詞：K分布;雷達圖像;最大似然估計;協方差矩陣

0引言

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

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

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

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

1極化SAR圖像協方差矩陣建模

1.1極化SAR圖像特性

極化SAR是一種主動式的微波成像設備,具有分辨率不受雷達傳感器高度影響的優點,其雷達信號能穿透云層。另外,波長較長的L波段信號能穿透森林和地表植被覆蓋,在軍事上能夠用于發現叢林中的隱藏目標和獲取地表以下的信息,這些是光學等其它傳感器無法實現的。

最常見的極化SAR為機載形式,美國NASA的JPL實驗室研制的AIRSAR系統可以在P、L和C三個波段上進行全極化測量成像。該系統自首飛后每年至少執行一次測量任務,成為驗證新開發雷達技術的主要實測數據來源。AIRSAR系統安裝于DC-8飛機的底部,如圖1所示。

圖1　AIRSAR系統

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

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

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

1.2極化SAR圖像K分布建模

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

C=τY,

(1)

式中: τ為紋理變量; Y為相干斑變量,并且紋理變量和相干斑變量相互統計獨立。一般認為相干斑變量Y服從Wishart分布,其PDF為[5]

(2)

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

(3)

式中,Γ(L)為標準歐拉函數。

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

(4)

式中,α為伽馬分布的形狀參數。

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

(5)

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

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

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

(10)

式中MeijerG函數的定義為

(11)

那么其最大似然估計可由對數似然比公式求取

(12)

求得

(13)

其中參數z為下式的解

(14)

3仿真實驗分析

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

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

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

(15)

則式(13)可表示為

(16)

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

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

圖3　最大似然修正系數仿真圖

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

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

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

3.2矩陣行列式值ML估計方法在極化SAR圖像ENL估計中的應用

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

(17)

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

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

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

2) 預先建立的是一個非均勻間隔索引的表,即L和α關于f(L;α)的二維索引表:當L≤4時,以間隔值ΔL=0.01建表;當L>4時,以ΔL=0.1建表;

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

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

3.3實測數據仿真分析

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

這里用無監督的滑窗估計方法來對實測數據的ENL值進行估計[4]。該方法首先用k×k的滑動窗口將整個圖像分解成很多個小區域,然后根據每個區域中的樣本數據計算出其對應的ENL值,最后對所有區域的ENL估計結果進行統計。由于上述兩幅極化SAR圖像中的大部分區域都可以用K分布來進行建模,這里取ENL統計分布圖的峰值作為整個圖像的ENL值。

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

圖6　實測數據下各估計方法對ENL的估計結果的比較(a) Flevoland; (b) San Francisco

4結束語

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

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

參考文獻

[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分布矩陣行列式值的統計特性及其在參數估計中的應用[J]. 電子學報, 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-),男,湖南人,碩士， 主要研究方向為通信技術與計算機應用。

崔浩貴(1987-),男,浙江人,博士,研究方向為雷達極化信息處理、電子戰建模與仿真。

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

中圖分類號：TN95

文獻標志碼：A

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

作者簡介

資助項目： 國家自然科學基金(批準號:61372165)

聯系人： 晏慶E-mail: 63934470@qq.com