毛滔 刘涛
(1.海军参谋部电子对抗雷达声纳局,北京 100841;2.海军工程大学电子工程学院,武汉 430033)
基于Fisher分布的极化合成孔径雷达统计建模及其参数估计方法
毛滔1刘涛2
(1.海军参谋部电子对抗雷达声纳局,北京 100841;2.海军工程大学电子工程学院,武汉 430033)
为更准确地描述高分辨极化合成孔径雷达(Synthetic Aperture Radar,SAR)图像的尖峰和长拖尾等统计分布特性,提出了基于Fisher分布的极化图像多变量乘积模型,并研究了其参数估计方法.首先给出了柯西分布相干斑噪声等效纹理分量的概率密度函数及其低阶矩特征;然后利用散射因子服从F分布的等效纹理变量与高斯散斑变量相乘形成的多变量乘积统计模型,得到了Fisher分布模型的概率密度函数,并推导了其多视协方差矩阵的概率密度函数和矩阵行列式值的低阶矩特征;最后提出了基于矩阵行列式值的矩估计和基于Mellin变换的对数累积量估计等两种参数估计方法,并进行了对比,同时通过仿真数据和实测数据验证了理论模型和新参数估计方法的有效性.这为高分辨极化SAR图像建模、目标检测和识别等领域的理论研究和工程实现提供了新途径.
Fisher分布;Kummer-U函数;极化SAR;Mellin变换;多视;建模估计
DOI 10.13443/j.cjors.2016030102
在极化合成孔径雷达(Synthetic Aperture Radar,SAR)图像解译中杂波模型的准确性起着重要的作用,是雷达目标检测、分类识别以及相干斑抑制的基础.雷达回波的统计分布一般服从满足有限方差中心极限定理的高斯分布(幅度瑞利分布)[1-2],但是随着现代雷达频率和分辨率的不断提高,需要借助非高斯分布函数来刻画地杂波和海杂波的统计特性.这种分辨单元内散射点数目较少或者有显著散射点的雷达回波往往具有长拖尾的分布特性,对数正态分布、Weibull分布、K分布常常用来描述这种海杂波和地杂波的起伏特性[3-4].在高分辨以及低掠角条件下Weibull杂波模型对回波的描述更加准确[5-6],同时在低分辨条件下平方根伽玛分布能较准确地描述杂波的幅度统计特性[7].对数正态乘积模型常常比观测数据有更强的拖尾,伽玛和Weibull分布的拖尾要比实际观测数据弱[8].K分布等基于纹理分量和散斑分量的多变量乘积模型在描述杂波模型上也备受关注,已经被广泛应用在长拖尾杂波的建模、信号处理以及数据分析等领域.乘积模型表明,在散射中心数目足够的条件下[9-10],雷达回波等价于满足有限方差中心极限定理的散斑高斯噪声与地面后向散射系数的乘积,目前存在的问题是相干斑噪声分量服从满足中心极限定理的高斯分布的这个假设未必总是成立.一般情况下如果被观测区域是均匀的,平方根伽玛乘积模型(也就是K分布模型)能较好地描述这一区域的随机分布[11].K分布等多变量模型是包含形状参数和尺度参数的两参数模型(强度归一化后为单参数模型),为了使多变量模型能够适用的场景更加广泛,更具一般性的G分布和Fisher分布等三参数模型(强度归一化后为双参数模型)备受关注[12-17],这实际上是以增加自由度和运算量来提高模型的适用度,但是存在的问题是某些统计模型的物理意义并不明确,参数估计复杂且不准确.同时,在高分辨雷达照射下由于陆海特别是海面散射回波随时间和空间的起伏变化较快,回波一般具有尖峰和长拖尾特性,尖峰特性的出现使得常规的K分布等多变量模型在信号处理中的可靠性下降[18-19].具有尖峰和长拖尾特性的稳定分布是满足广义中心极限定理的唯一的一类分布,并且能够保持自然噪声过程的产生机制和传播条件,从而在统计模型领域备受关注[19].R.D.Pierce已经将对称稳定分布模型用于描述海面舰船和杂波的散射截面积,并给出了微弱舰船目标的检测算法,取得了较好的效果[17].Kuruoglu等的研究结果表明基于稳定分布的重尾瑞利分布(heavy-tailed Rayleigh distribution)非常适合从物理原理上描述雷达波散射特性、回波长拖尾及其尖峰特性[20-21],但是稳定分布不具有封闭的概率密度函数解析形式,造成数值计算量大,也不利于进行深入的性能分析,造成了其实际应用的困难.柯西分布是具有封闭表达式的仅有的两种对称稳定分布之一(另一种是高斯分布),同时柯西分布的信号处理算法对对称稳定模型误差具有很好的鲁棒性和优于高斯算法的性能,使得柯西分布模型能够较好地近似其他对称稳定分布[19].借鉴不同杂波分布的各种统计模型,在文献[17, 22]的理论基础上,结合乘积相干斑模型各分量的实际物理意义,本文推导了服从伽玛分布和柯西分布的双变量乘积统计模型的概率密度函数的解析表达式.然后根据等效纹理分量服从Fisher分布的一般情况,将极化SAR杂波的柯西乘积模型拓展为Fisher分布的多变量乘积统计模型,并分析了其协方差矩阵行列式值的矩特征.
与此同时,快速准确的参数估计是保证统计模型有效的关键技术.统计模型参数的估计问题是SAR图像分割、分类等统计分析和信号处理的关键步骤和经典问题.参数估计常用的方法是矩估计和最大似然估计.最大似然估计在求解复杂的统计函数时计算可能非常复杂,而不得不采用数值计算的方法,运算量很大不利于快速求解;矩估计方法局限于高阶矩或者分数阶矩存在的情况下,同时矩估计方法受噪声的影响较大,较小的噪声扰动也会引起参数估计较大的偏差,并且和最大似然估计一样,矩估计也有可能产生较为复杂的数学表达式而造成计算复杂度的增加.最近的研究表明基于Mellin变换的对数累积量参数估计方法与矩估计和最大似然估计相比具有计算简单、精度高的优点,它在SAR图像参数估计领域的突出表现也引起了研究人员的兴趣[22-23].S N Anfinsen将研究方法成功应用到极化SAR图像等效视图数估计、模型拟合度分析和目标检测等领域,取得了相当理想的结果[22-24].因此本文提出了基于Mellin变换的对数累积量估计方法.新的处理方法有效提高了极化Fisher分布参数估计的准确性和高效性.
复高斯矢量Y的幅度服从平方根伽马分布,其概率密度函数为
y≥0, v>0, a>0.
(1)
(2)
(3)
a>0, t>0, v>0.
(4)
经计算等效后向散射幅度的r阶矩为
(5)
(6)
此分布函数是F分布函数的线性变形,由式(5)可知其具有(0,1/2)阶矩.可见基于柯西分布的乘积模型可以用F分布函数来表征.实际上,F分布也包含了高斯分布和柯西分布的强度分布形式,可以看作是对称稳定分布的一个简单替代.在噪声分量服从F分布时,采用上述方法可以推导得到等效后向散射强度的概率密度函数为
(7)
式中,u=n/2,n为逆伽玛分布的自由度, 直接将此分布作为乘积模型的纹理分量在单极化SAR图像里面已经有学者提出[15],L Bombrun将其应用到全极化图像并进行图像分割,取得了较好的效果,尽管他们并没有说明模型的物理意义.如何快速准确地获得模型参数仍然是目前极化Fisher分布实际应用的难点问题,也是本文后面需要解决的问题.
基于Mellin变换的对数累积量参数估计方法是目前广泛关注的有效参数估计方法,特别适用于乘积模型,具体方法描述详见文献[22,26].可以得到F分布的矩特征为[15]
(8)
基于Mellin变换的F分布的对数累积量为[15]
(9)
(10)
式中:ψ为普西函数;r=2,3,….
极化SAR中,复散射矩阵用S表示,记为
(11)
在互易条件下,SHV=SVH,这样散射矩阵S就可以用目标矢量来代替:
(12)
(13)
复高斯矢量服从以下分布:
(14)
d是复高斯矢量的维数.
为了抑制相干斑的影响,极化数据往往通过平均方差矩阵的方法来进行多视处理.从式(13)对方差矩阵求均值得到
(15)
式中:L为视图数;Z(k),Y(k),T(k)分别是第k个Z,Y,T的单视采样点;上标“H”表示共轭转置.如果视图数不是很多的话,在多视处理中T(k)独立于k[12], 在此假设下,式(15)转化为
(16)
(17)
式中:C=E(yyH);d是复矢量ZC的维数;tr(·)为行列式的迹;函数Γ(L,d)为
(18)
(19)
(20)
将式(7)和式(20)代入到式(19),交换积分次序,利用拉普拉斯变换[29]以及Kummer-U合流超几何函数定义[30]
exp(-pz)dz,
(21)
式(19)可化简为
U(Ld+u,Ld-v+1,Ltr(C-1z)),
(22)
则Fisher分布乘积模型方差矩阵概率密度函数为
U(Ld+u,Ld-v+1,Ltr(C-1z)).
(23)
高斯分布行列式值的高阶矩特征为[28]
(24)
那么易得Fisher分布极化协方差矩阵行列式值的高阶矩为
(25)
在只考虑同极化参数的强度分布情况下,即在式(25)中令d=1,通过交换积分次序易得其一维多视强度的分数阶矩为
r (26) 式(26)利用贝塔函数可重新表示为 r (27) k1{Z}=ψd(0,L)+ln|C|+d(ψ(v)- ψ(u)-ln L) , (28) kr{Z}=ψd(r-1,L)+dr(ψ(r-1,v)+ (-1)rψ(r-1,u)) (29) 所示,这是后续利用对数累积量进行参数估计的核心公式之一. (30) 多通道矩估计平均值的方法本质上只利用了共极化通道数据,并没有充分利用极化协方差矩阵的相关信息,必定会导致估计误差的增大,同时上述估计方法受阶数的选取影响较大(阶数选取对估计精度的影响依据见图1).因此我们提出了基于协方差矩阵行列式值的矩估计方法. 由式(25)可得L视强度的r/2阶矩特征为 (31) 那么易得 (32) 上述估计方法受阶数的选取影响也较大.为了充分利用极化协方差矩阵信息并避免阶数对结果的影响,可以利用式(28)~(29)中的二三阶累积量进行参数估计,即 (33) (34) 对数累积量的计算形式如下(四阶以下对数累积量等于其对应对数中心矩): (35) (36) (37) Goodman,Kuruoglu发现柯西分布非常适合描述雷达高分辨散射和相干斑特性,及其拖尾和尖峰特性[10,22],本节将从仿真数据和实测数据两个方面分析验证极化F分布理论推导的正确性和有效性. 4.1 仿真数据分析 为了仿真简单易行,假设极化SAR图像HH通道数据服从变形t分布模型,仿真算法借鉴文献[14]所述方法.设形状参数u=v=0.5,极化SAR的HH通道尺度参数|CHH|=0.1.通过仿真数据,可以验证概率密度和矩特征理论推导的正确性,如图1所示. 从图1可以看出,理论推导结果完全符合仿真数据的统计,但是由于其长拖尾无限均值方差的特性,导致其幅度矩特征在1附近产生振荡,与理论值相差较大,但是随着样本数的增加,仿真结果与理论值逐渐接近,这与仿真数据的无限均值、无限方差的特性是相符合的.另外仿真直方图尾部的高直特点恰恰反映了该分布具有的重尾特性.其低阶矩特征也说明,在样本数目一定和计算精度许可的条件下,低阶矩阶数的减小有助于参数矩估计精度的提高.对本文提出的三种估计方法,即单通道平均矩估计方法、基于矩阵行列式值的矩估计方法和基于矩阵行列式值的对数累积量估计方法,在L=4,u=2.5,v=3.5参数假定下,进行了性能的仿真比较,结果如表1所示.可以看出对数累积量的估计方法要比矩估计方法准确快速.矩估计随着矩阶数的减小,估计精度总体有所提高. 图1 仿真数据与理论结果的比较 表1 三种估计方法性能比较(样本数1 024,仿真次数1 000) 4.2 实测数据分析 这里采用的极化SAR图像(图2)是Pol-SARPro软件中常用的由ESAR/OP-AIRFIELD-提供的Oberpfaffenhofen全极化雷达图像(3m×3m),其中名义视图数为L=12,以及Pol-SARPro软件中AIRSAR的SanFrancisco全极化雷达图像分辨率(10m×10m)4视数据.为了验证极化F分布理论与实测数据的拟合程度,采取K-S非参数检验方法[5]和最小方差(MinimumSquareError,MSE)准则[15]计算二者的差异, 结果如表2所示.表3为显著性水平为α时的K-S检验门限允许值. 区域uv极化通道MSEK-S检验HH0.54071.3113A(城区)2.960010.3900HV0.05160.7052VV0.05010.5474HH0.22581.2133B(植被)+∞+∞HV0.62021.3527VV0.11830.7724HH0.05680.6762C(森林)+∞8.6200HV0.78581.3236VV0.06200.6853HH0.19571.0004D(城区)6.80003.7000HV0.57611.3038VV0.20641.1625 表3 显著性水平为α时的K-S检验门限值 从表2和表3可以看出:A和D城区可以方便地用极化Fisher分布模型进行拟合;B区植被的两个参数都为无穷大,说明B区植被协方差矩阵服从Wishart分布,也就是相干斑服从复高斯分布,这也是符合以前的理论结果的;C区森林服从K分布,这也与其典型分布相吻合.从表中也可以看出不管在森林、农作物,还是城区,极化Fisher分布都具有很小的K-S值和MSE值,表明极化Fisher分布模型能够很好地拟合大部分区域,包括高斯分布和K分布以及描述不均匀区域的G0分布(v→∞). 从表4可以看出:在San Francisco的海洋区域和森林区域,尺度参数数值都非常大,符合高斯分布或者K分布模型假设,这与实际情况比较相符;而在城区,形状参数数值非常大,这与非均匀区域G0的分布假设非常吻合,因此本文提出的Fisher分布模型具有一定的普适性,可以较好地对均匀和非均匀区域进行合理的建模与参数估计,为后续的目标检测与识别奠定了基础. 表4 San Francisco四个区域的参数估计与检验分析 本文提出了适用范围广泛的基于多变量乘积的极化Fisher分布数学模型,它具有闭合的解析表达式.通过研究极化F分布的矩特征和对数累积量特征,提出了模型参数基于矩阵行列式值的矩估计方法及基于Mellin变换的对数累积量估计方法.在多参数估计的情况下给出了基于对数累积量的极化Fisher统计分布数据估计步骤,部分避免了多元方程遍历的低效解法.研究成果为高分辨极化SAR图像建模、目标检测和识别等领域的理论研究和工程实现提供了新思路.由于基于高阶对数累积量方法的参数估计的准确性与样本数目相关性较大,容易引起小样本时参数估计误差较大甚至无解,如何解决小样本条件下的参数准确有效估计是下一步的研究重点. [1]刘涛, 黄高明, 王雪松, 等.基于H/A目标极化分解理论的统计分析[J].中国科学(F辑), 2010, 40(1):102-114. LIU T, HUANG G M, WANG X S, et al.Statistical analysis of the H/A target decomposition theory[J].Scientia sinica informationis, 2010,40(1):102-114.(in Chinese) [2]LIU T, WANG X S, XIAO S P.Statistical characteristics of the normalized Stokes parameters [J].Science China information sciences, 2008, 51(10):1594-1606. [3]黄培康, 殷红成, 许小剑.雷达目标特性[M].北京:电子工业出版社, 2005. [4]YACOUB M D, FRAIDENRAICH G.Joint statistics for two correlated Weibull variates[J].IEEE antennas &wireless propagation letters, 2005,4:129-132. [5]LIU T, HUANG G M, WANG X S, et al.Statistics of the polarimetric Weibull-distributed electromagnetic wave [J].IEEE transactions on antennas and propagation, 2009, 57(10):3232-3248. [6]FARINA A, RUSSO A, SCANNAPIECO F, et al.Theory of radar detection in coherent Weibull clutter [J].IEEE proceedings, 1987,134(2):174-190. [7]MEYER D P, MAYER H A.Radar target detection-handbook of theory and practice[M].Academic Press, 1973, 64-82. [8]IRVING W W, OWIRKA G J, NOVAK L M.Adaptive processing of polarimetric SAR imagery[C]//24th Asilomar Conference on Signals, Systems, and Computers.Pacific Grove, CA, 1990. [9]TUR M, CHIN K C, GOODMAN J W.When is speckle noise multiplicative?[J].Applied Optics, 1982, 21(7):1157-1159. [10]GOODMAN J W.Statistical optics[M].New York:John Wiley&Sons, 2000. [11]FRERY A C, MULLER H J, YANASSE C C F.A model for extremely heterogeneous clutter[J].IEEE transactions on geoscience and remote sensing, 1997, 35(3):648-659. [12]刘涛, 王雪松, 肖顺平.极化SAR图像处理中L分布杂波统计特性分析[J].自然科学进展,2009, 19(1):89-96. LIU T, WANG X S, XIAO S P.Statistical analysis of L-distributed clutter in Pol-SAR imagery [J].Progress in nature and science, 2009, 19(1):89-96.(in Chinese) [13]FREITAS C C, FRERY A C, CORREIA A H.The polarimetric G distribution for SAR data analysis[J].Environmetrics, 2005, 16:13-31. [14]MARIER L J, Jr.Correlated K-distributed clutter generation for radar detection and track[J].IEEE transcations on aerospace and electronic systems, 1995,31(2):568-580. [15]TISON C, NICOLAS J M, TUPIN F, et al.A new statistical model for Markovian classication of urban areas in high resolution SAR images[J].IEEE transactions on geoscience and remote sensing, 2004, 42:2046-2057. [16]ACHIM A, TSAKALIDES P, BEZERIANOS A.SAR image denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[J].IEEE transactions on geoscience and remote sensing, 2003, 41:1773-1784. [17]PIERCE R D.RCS characterization using the alpha-stable distribution[C]//Proceeding of IEEE National Radar Conference.Ann Arbor, MI, 13-16 May, 1996:154-159. [18]PANAGIOTIS T, CHRYSOSTOMOS L N.Maximum likehood localization of source in noise modeled as a stable process[J].IEEE transactions on signal processing, 1995, 43 (11):2700-2713. [19]邱天爽, 张旭秀, 李小兵, 等.统计信号处理-非高斯信号处理及其应用[M].北京:电子工业出版社, 2004. [20]ACHIM A, KURUOGLU E E, ZERUBIA J.SAR image filtering based on the heavy-tailed Rayleigh model [J].IEEE transactions image processing, 2006, 15(9):2686-2693. [21]KURUOGLU E E, ZERUBIA J.Modeling SAR images with a generalization of the Rayleigh distribution[J].IEEE transactions on image processing, 2004, 13:527-533. [22]ANFINSEN S N, 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. [23]ANFINSEN S N, DOULGERIS A P, ELTOFT T.Estimation of the equivalent number of looks in polarimetric synthetic aperture radar imagery[J].IEEE transactions on geoscience remote sensing, 2009, 47(11):3795-3809. [24]时公涛, 赵凌君, 桂琳, 等.基于Mellin变换的K 分布参数估计新方法[J].电子学报, 2010, 38(9):2083-2089. SHI G T, ZHAO L J, GUI L, et al.Novel parameter estimation of K-distribution based on Mellin Transform[J].Acta electronica sinica, 2010, 38(9):2083-2089.(in Chinese) [25]SAMORODNITSKY G, TAQQU M S.Stable non-Gaussian random process:Stochastic models with infinity variance[M].New York:Chapman and Hall, 1994. [26]JAMALIZADEH A, BALAKRISHNAN N.On a generalization of bivariable Cauchy distribution[J].Communications in statistics-theory and methods, 2008,37(4):469-474. [27]APELBLAT A.Table of definite and infinite integrals[M].Amsterdam:Elsevier Science Publishing Company, 1983. [28]GOODMAN N R.Statistical analysis based on a certain complex Gaussian distribution (an introduction)[J].Annals of mathematical statistics, 1963, 34:152-177. [29]BATEMAN H.Tables of integral transforms, volume I[M].McGraw-Hill Book Company, Inc.1954 [30]ABRAMOWITZ M, STEGUN I A.Handbook of mathematical functions with formulas, graphs, and mathematical tables, national bureau of standards applied mathematics series[M].Washington D C:U.S.Government Printing Office, 1964. [31]DEAN~O A, TEMME N M.On modified asymptotic series involving confluent hypergeometric functions[J].Electronic transactions on numerical analysis, 2009, 35:88-103. [32]数学手册编写组.数学手册[M].北京:高等教育出版社, 2001. 毛滔 (1979-),男,江西人,海军联合海情中心高级工程师,博士,研究方向为地波雷达、雷达信号处理和数据处理. 刘涛 (1978-),男,山东人,海军工程大学电子工程学院教授,博士生导师,研究方向为雷达极化信息处理、雷达电子战建模与仿真. Modeling polarimetric SAR image based on Fisher distribution and its parameter estimation MAO Tao1LIU Tao2 (1.TheFourthDepartmentofHeadquartersofNavy,Beijing100841,China;2.NavalUniversityofEngineering,Wuhan430033,China) Based on the Fisher multiple variable product model, the Fisher statistical model of polarimetric synthetic aperture radar(SAR) images is put forward to describe the properties of heavy tail and peak, and the statistics and the estimation methods are obtained.Firstly, the probability density function (PDF) of developed t distribution is given based on Cauchy distribution, so are the fractal moments.Secondly, the Fisher distribution model is obtained through the multiple variable product model mixed with the Fisher variable, and the PDF and the fractal moments of the covariance matrix moments are derived.At last, two estimation methods via Mellin transform and Matrix moments are presented.The performances of the novel data model and novel parameter estimation methods are verified by the simulated data and real data.The research provides a simple analytical method to describe different distributed clutter, which is useful to target detection and recognition. Fisher distribution;Kummer-U function;polarimetric SAR;Mellin transform;multi-look;modeling;estimation 毛滔, 刘涛.基于Fisher分布的极化合成孔径雷达统计建模及其参数估计方法[J].电波科学学报,2016,31(5):948-956. 10.13443/j.cjors.2016030102 MAO T, LIU T.Modeling polarimetric SAR image based on Fisher distribution and its parameter estimation[J].Chinese journal of radio science,2016,31(5):948-956.(in Chinese).DOI:10.13443/j.cjors.2016030102 2016-03-01 国家自然科学基金(61372165) TN957.52 A 1005-0388(2016)05-0948-09 联系人:刘涛 E-mail:liutao1018@sina3 极化Fisher分布乘积模型的参数估计方法
4 极化SAR图像数据分析
5 结 论