专利名称：多视全极化合成孔径雷达数据压缩方法
技术领域：
本发明涉及极化SAR(Synthetic Aperture Radar，合成孔径雷达)技术领域，尤其涉及一种基于极化相干矩阵的多视全极化SAR数据压缩方法。
背景技术：
SAR已经成为地学遥感领域中信息获得的一个重要来源，得到了广泛的应用，现在已有许多SAR系统正在实际使用运行中，如日本的ALOS卫星、德国的Terra-SAR、加拿大的Radarsat-2等等。
普通SAR只有一个极化通道，电磁波水平发射水平接收或垂直发射垂直接收，分别用HH和VV表示。单极化通道测得的数据为复数据，包含实部和虚部两部分。
全极化是SAR发展的一个重要趋势，全极化SAR(简称为POLSAR)拥有四个极化通道，分别为电磁波水平发射后水平接收和垂直接收，垂直发射后垂直接收和水平接收，分别用HH、VH、VV、HV表示，每个通道测得的也均为复数据。
通道数的增多在提供了更多观测目标信息的同时，也带来了数据量的显著增加。普通单通道SAR的数据量已十分巨大，四通道全极化SAR的数据量则是普通单通道SAR的4倍。为使数据的存储和传输更加便利，有必要对成像后的POLSAR数据进行压缩。全极化SAR数据压缩方法要求在减少数据量的同时，应尽可能减少引入的误差，并保持极化信息不变。其压缩对象为POLSAR一次观测成像后所得到的一幅全极化SAR图像数据。
为了降低SAR图像固有的斑点噪声，SAR系统在成像时都会采用多视处理。通常SAR系统的方位分辨率高于距离分辨率，可以在方位向进行多视处理。多视处理就是在成像方法最后进行方位向逆傅里叶变换成像之前，将数据在方位向频谱上分成几个带宽相等的独立部分，对每部分分别进行逆傅里叶变换成像，然后对各部份数据进行非相干叠加生成多视全极化SAR图像。
多视全极化SAR图像中每个象素点都对应一组全极化参数。多视全极化参数有多种表示形式，常见的有极化相关矩阵(简称为C)、极化相干矩阵(简称为T)、Kennaugh矩阵(简称为K)。C、T、K与HH、VH、VV、HV四通道数据具体关系如下 极化相干矩阵C 其中，上标t和*分别表示转置与共轭，N为多视处理的视数，即多视处理时将数据在方位向频谱上分成了N份。
极化相干矩阵T T＝Q1·C·Q1H/2， 其中，上标H表示共轭转置。
Kennaugh矩阵K 其中，
为Kronecker积。
目前，多数POLSAR系统电磁波的发射和接收都采用同一部天线，理论上VH与HV通道的数据应相等，因此POLSAR通常进行如下操作
其中φ为HV与VH通道的平均相位差，通常为0。(4)式也可理解为使得VH＝HV，这也是(1)和(2)式中计算C、T时并没有涉及VH通道数据的原因。在VH＝HV的情况下C、T均为3×3的非负定Hermitian矩阵，K为4×4实对称矩阵。
由上述可知C、T、K矩阵都是由POLSAR四个通道数据计算出来的，因此它们是等价的，相互之间可以转化，知道其中一种表示形式就可计算出另两种。
1993年1月，G.De Grandi，C.Lavalle和A.J.Sieber发表了论文《A coding algorithm for the covariance matrix representation ofpolarimetric radar data》，提出了一种基于极化相关矩阵的极化SAR数据压缩方法。该方法假设HH、VV通道与HV通道的共轭相乘为0，对于对称目标具有很好的压缩效果，但绝大部分实际数据不能满足这个假设，所以该方法应用十分有限，没有得到推广使用。
1989年6月，P.Dubois，L.Norikane，J.J.van Zyl和H.Zebker发表了论文《Data volume reduction for imaging radar polarimeter》，提出了一种对于Kennaugh矩阵的10字节压缩方法(简称为K10BC)。K10BC根据每个象素点自身极化数据间的相互关系进行压缩，对象素点的极化特征保持较好，压缩精度较高。压缩所引入的误差仅为量化误差，并不影响后续的滤波、分类、目标检测等操作。该方法已经成功的应用于美国NASA/JPL的机载AIRSAR系统。K10BC方法的缺点是压缩比固定且较低，同时压缩精度仍有待提高。K10BC方法主要对K矩阵表示的极化数据进行压缩，K矩阵与雷达的接收功率直接相关，早期对雷达极化的研究比较关注雷达的接收功率，主要希望利用改变雷达的极化方式达到最大的接收功率。而目前绝大部分的POLSAR数据处理过程，如滤波、分类、目标检测等，都涉及极化数据的概率统计分布。而C与T矩阵直接与极化数据的统计分布相关，其中尤以T矩阵使用最为广泛。K10BC方法需先解压出K矩阵再转换为T矩阵，此过程必然会造成一定程度的误差扩散，使得T精度较差。
发明内容
本发明的目的是提供一种全极化SAR数据压缩方法，以解决现有技术中极化SAR数据压缩方法存在的上述缺陷。
为了达到上述目的，本发明的技术方案提出一种多视全极化合成孔径雷达数据压缩方法，该方法包括以下步骤 读入多视全极化合成孔径雷达图像每一象素点对应的极化相干矩阵； 根据所述极化相干矩阵对极化总功率进行动态范围可变的压缩； 利用所述动态范围可变压缩的解压值对所述极化相干矩阵的各元素进行归一化处理。
上述多视全极化合成孔径雷达数据压缩方法中，所述极化相干矩阵T为3×3非负定厄米特矩阵，如式(5)所示 其中，*代表复共轭，T11、T22、T33均为非负实数，T12、T13、T23为三个复数；则所述极化总功率与极化相干矩阵的关系如式(6)所示 Aspan＝T11+T22+T33 (6) 其中，Aspan即为极化总功率。
上述多视全极化合成孔径雷达数据压缩方法中，所述对极化总功率进行动态范围可变的压缩具体包括 统计所述极化总功率的最大值max(Aspan)与最小值min(Aspan)； 利用式(7)计算压缩所需的参数α和β 使用求得的β值根据式(8)对Aspan值进行调整 根据式(9)使用两字节对调整后的Aspan值进行存储
(9) 其中，
表示不大于logαAspan的最大整数。
上述多视全极化合成孔径雷达数据压缩方法中，所述对极化相干矩阵各元素进行归一化处理具体包括 利用式(10)计算Aspan的解压值x 利用x值根据式(11)对极化相干矩阵的对角线元素T22和T33进行归一化处理 (11) 利用x值根据式(12)对极化相干矩阵的右上角复元素T12、T13和T23进行实部、虚部分别归一化处理 (12) 其中，real()为取复数实部，imag()为取复数虚部。
上述多视全极化合成孔径雷达数据压缩方法中，该方法还包括压缩之后的数据存储步骤 将byte1先加上128，再采用无符号8比特整形数据对byte1～byte10进行存储； 采用double型数据对参数α和β进行存储。
上述多视全极化合成孔径雷达数据压缩方法中，该数据压缩方法对应的数据解压缩方法包括以下步骤 读入压缩后的10字节无符号8比特整形数据，将每字节赋值给byte1～byte10，并对byte1减去128； 读入α和β的值，根据式(10)计算解压值x 根据式(13)解压极化相干矩阵的对角线元素T22和T33 T22＝(byte3+0.5)2/65536·x(13) T33＝(byte4+0.5)2/65536·x 根据式(14)解压极化相干矩阵的对角线元素T11 T11＝x-T22-T33 (14) 根据式(15)解压极化相干矩阵的右上角复元素T12、T13和T23 本发明的技术方案采用极化相干矩阵表示全极化参数，并对极化总功率进行动态范围可变的压缩；不使用极化总功率，而使用它的解压值对极化相干矩阵的元素进行归一化；同时对复数采用虚部、实部分别量化的方法；使该方案在压缩精度、信噪比、保证极化相干矩阵特征值非负和极化特征保持等方面都要优于现有方法，并且运算效率高、速度快，具有良好的实时性和工程应用价值。
图1为本发明多视全极化SAR数据压缩方法实施例流程图； 图2为图1所示实施例对应的解压方法流程图； 图3为采用现有技术的K10BC方法对Aspan进行压缩的相对误差直方图； 图4为采用本发明方法对Aspan进行压缩的相对误差直方图。
具体实施例方式 以下实施例用于说明本发明，但不用来限制本发明的范围。
一幅多视POLSAR图像每个象素点对应的极化相干矩阵T为3×3非负定Hermitian矩阵，如下式所示 其中，*代表复共轭，T11、T22、T33均为非负实数，T12、T13、T23为三个复数，T矩阵左下角元素为右上角元素的复共轭。由公式(5)可知T共包含9个独立的实变量，本发明提供的压缩方法将利用10个字节存储每个象素点的T矩阵，以下简称T10BC方法。
本发明T10BC压缩方法的实施例如图1所示，包括以下步骤 S101、读入多视极化SAR图像每一象素点对应的极化相干矩阵数据。
S102、根据待压缩数据的动态范围求取压缩所需参数； 对于极化参数存在如下的关系 Aspan＝T11+T22+T33(6) 其中，Aspan为极化总功率，是极化雷达中非常重要的一个参数。首先统计待压缩数据的动态范围，即Aspan的最大值与最小值，然后计算α和β两个参数如下 其中max()为选取待压缩POLSAR图像中括号内参数的最大值，min()为选取最小值。
S103、根据α和β压缩Aspan值； 首先，利用β参数调整Aspan的动态范围，方法如下 经上式调整后的Aspan值，指数部分正、负范围相等。本实施例的方法利用2字节存储Aspan值，1字节存储指数，1字节存储尾数，方法如下
(9) 其中，
表示不大于logα Aspan的最大整数。
对数底使用α，即(7)式中根据数据实际动态范围计算出的一个合适的对数底，一般α值很接近于1。
S104、计算Aspan的解压值； 计算Aspan的解压值用x表示，如下 x将用来对T矩阵中的其它元素进行归一化。由于解压时x可精确算得，所以用x进行归一化将比直接用Aspan的误差更小，进而提高压缩精度。
S105、压缩极化相干矩阵的对角线元素T22和T33； T22和T33为T矩阵对角线上的元素，均为非负实数，压缩方法如下 (11) T22和T33经x归一化后，值域范围为
，分布主要集中于0附近，所以采用开方操作调整分布后再进行压缩，以提高压缩精度。
S106、压缩极化相干矩阵的右上角复元素T12、T13和T23； T12、T13和T23为T矩阵右上角元素，均为复数，压缩方法如下 (12) 其中real( )为取复数实部，imag( )为取复数虚部。由于T12、T13和T23的实部、虚部相对Aspan来说都较小，所以先将它们的动态范围都调整到
再进行开方操作调整分布。经试验发现，对于复数T12、T13和T23采用实部、虚部分别量化的压缩精度要高于采用模值、相位的压缩方法，因此T10BC选用实部、虚部分别量化的方法。
S107、数据存储； 数据存储时byte1先加上128，然后byte1～10用uint8型，即无符号8比特整形数据存储。参数α和β采用double型数据存储。
由式(11)和(12)可知，对T矩阵中各元素采用的是8比特均匀量化，因此T10BC方法引入的误差为量化误差。8比特量化精度较高，压缩比适中，8比特为1字节，各种编程语言如C语言、Matlab等都可以很方便的直接对字节操作，编程实现容易、运行效率高。
但量化比特数并非必须固定为8比特不变。如果需要高的数据精度可以采用较多比特的量化，如10比特、12比特均匀量化，这将提高压缩精度，但降低压缩比；若需要较高的压缩比可采用较少比特的量化，如6比特、4比特均匀量化，这将降低压缩精度，提高压缩比。采用非8比特量化实际编程实现时，需要对数据进行按位操作，编程实现时将会较复杂，同时降低一定的运行效率。量化比特数的具体选择应根据实际应用需求确定。
T矩阵包含9个独立实变量，未经压缩时如果用float型数据存储，则每个T矩阵需要9×4＝36字节。若T10BC采用8比特量化，则压缩后每个T矩阵用10字节表示，则压缩比为36/10＝3.6。
上面介绍的是T10BC压缩方法，相应的解压方法描述如下 S201、数据读入； 读入压缩后的10字节uint8型数据，并把每字节赋值给byte1～10，并对byte1减去128。
S202、计算极化总功率Aspan的解压值x； 读入参数α和β，并根据上式(10)计算出x x值将被用来解压T矩阵的其它元素。
S203、解压极化相干矩阵的对角线元素T11、T22和T33； 先解压出T22和T33，方法如下 T22＝(byte3+0.5)2/65536·x(13) T33＝(byte4+0.5)2/65536·x 然后，根据极化参数之间的关系解压出T11，方法如下 T11＝x-T22-T33(14) S204、解压极化相干矩阵的右上角复元素T12、T13和T23； 解压T12、T13和T23的具体方法如下 以下结合另一具体实施例对上述本发明的T10BC压缩方法进一步加以阐述。
本实施例的压缩数据，选用德国机载全极化E-SAR系统对德国Oberpfaffenhofen地区进行观测的多视极化数据，图像大小为1300×1200。
量化比特数选用8比特，即T10BC方法，压缩过程如下 1)读入原float型极化相干矩阵数据； 2)利用式(7)计算出α和β两个参数； 3)根据α和β利用式(8)和式(9)压缩Aspan值； 4)利用式(10)计算出Aspan的解压值x； 5)利用式(11)压缩T22、T33； 6)利用式(12)压缩T12、T13和T23； 7)将byte1加上128，然后byte1～10用uint8型数据进行存储，α和β采用double型数据存储。
相应的解压缩过程如下 1)读入uint8型压缩后数据分别赋值给byte1～10，并对byte1减去128； 2)利用式(10)计算出每个象素点对应的x； 3)利用式(13)和式(14)解压出T矩阵的对角线元素T22、T33和T11； 4)利用式(15)解压出T矩阵右上角的复元素T12、T13和T23。
Oberpfaffenhofen地区全极化数据原大小为56160000字节。T10BC与K10BC方法压缩后数据大小均为15600000字节，压缩比3.6。其中T10BC采用的是8比特量化方法。
为了更好的说明本发明方法在压缩精度及极化特征保持上的性能，下面将对T10BC方法与AIRSAR系统使用的K10BC方法进行4个方面的比较。
1).Aspan压缩精度比较 Aspan是极化中十分重要的一个参数，因此对其压缩精度进行单独比较。
K10BC方法中对于Aspan压缩采用的是固定以2为底的指数量化。而本发明提供的T10BC方法引入α和β两个参数对Aspan采用动态范围可变的压缩。用两种方法分别对数据进行压缩，比较x值，即Aspan解压后的值，误差为(Aspan-x)。由于Oberpfaffenhofen地区全极化数据中包含城镇、森林、农田等各种目标，Aspan本身的动态范围就较大，所以统计误差的绝对大小并没有实际意义。下面首先给出相对误差如式(16)的统计直方图。
K10BC方法ER直方图如图3所示，本发明的T10BC方法ER直方图如图4所示。由图3和4可看出T10BC对Aspan压缩的相对误差要远小于K10BC方法。
分别计算两种方法对Aspan压缩的信噪比。信噪比如下式所示 根据上式计算得到两种方法对Aspan压缩的SNR分别为 SNRK10BC＝56.081dB，SNRT10BC＝86.197dB。
可见通过引入α和β两个参数使得对Aspan的压缩精度提高了30个dB。
2).整体压缩精度比较 计算所有数据的信噪比，以比较两种方法整体的压缩精度。信噪比是一个比较方法压缩精度的标准参数，如下式所示 其中，mean()代表对图中所有象素对应的该参数求平均，data表示原数据，datal表示解压后的数据。对于K矩阵，data包括K11、K12、K13、K14、K23、K24、K33、K34和K44九个K矩阵中的元素，这些是在K10BC方法中要用到的独立参量。对于T矩阵表示的数据，data包括T11、T22、T33、T12、T13和T23。SNR计算结果如表1所示 表1 两种压缩方法的信噪比
注下标K、T分别代表用K或T矩阵表示的极化数据 由表1中数据可知，T10BC方法整体压缩精度要稍优于K10BC方法。
3).保证极化相干矩阵特征值非负 相干矩阵特T的征值分解在极化数据处理中占有重要地位。T为3×3非负定的Hermitian矩阵，三个特征值应全为非负实数。压缩引入量化误差后，虽然误差很小仍会使一些较小的特征值变为负数。仍以德国Oberpfaffenhofen地区的数据为例，原数据经过特征值分解后特征值全为正数。经过K10BC和T10BC压缩后，再进行特征值分解，特征值变为负的个数如下表2所示 表2 对Oberpfaffenhofen地区的E-SAR数据 进行K10BC与T10BC压缩后特征值为负的个数 其中，|λ1|≥|λ2|≥|λ3|
由表2可知在保持特征值非负这一特性上T10BC要优于K10BC。
4)基于Signature的极化特征保持分析 全极化SAR数据是有实际物理意义的，在1)和2)中将极化SAR数据只作为数字单纯分析数据压缩后引入的信噪比并不全面。下面将结合极化的Signature理论研究压缩方法对极化特征的保持性能。
极化Signature能很好的描述各种目标的散射和极化特征，对于一个给定的目标，设其Kennaugh矩阵为K，目标的Signature定义如下 g(ε，τ)＝[1，cos 2τ cos 2ε，sin 2τ cos 2ε，sin 2ε]t 其中，上标t表示转置，椭圆倾角[-π/2，π/2]、椭圆率角[-π/4，π/4]，在共极化(简称为C-POL)方式下h(ε，τ)＝g(ε，τ)，在交叉极化(简称为X-POL)方式下 h(ε，τ)＝[1，-cos 2τ cos 2ε，-sin 2τ cos 2ε，-sin 2ε]t(20) 以共极化为例，令P1(ε，τ)与P2(ε，τ)分别为原数据与经过压缩后数据所对应的Signature。则计算Mean Relative Error of Signatures(简称为MRES)作为比较指标如下式所示，此类指标常用来分析压缩方法对全极化数据的极化特征保持能力。
其中q＝[K11 K12 K13 K14 K22 K23 K24 K33 K34]t，Kij为Kenn-augh矩阵元素，Wc为一加权矩阵。(21)式中分子为两个Signature之间的距离，然后用分母进行归一化，即为Signature的相对误差。对于交叉极化情况下计算MRES时，只需把(21)式中加权矩阵Wc替换为Wx。Wc、Wx如下 分别对K10BC，T10BC计算MRES结果如下表3所示 表3 两种方法的MRES参数比较 由表3可知T10BC对于极化特征的保持要优于K10BC方法。
综上所述，本发明提供的基于极化相干矩阵的多视全极化SAR数据压缩方法，由于采用了如下几项操作，使其在压缩精度、信噪比、保证极化相干矩阵特征值非负和极化特征保持等方面都要优于现有方法 1)采用极化相干矩阵表示全极化参数； 2)利用α和β两个参数对极化总功率进行动态范围可变的压缩； 3)不使用极化总功率，而使用它的解压值对极化相干矩阵的元素进行归一化； 4)对复数采用虚部、实部分别量化的方法； 5)对归一化后的T矩阵元素通过开方操作调整分布之后再进行均匀量化； 6)量化比特数可根据实际需求动态调整。
而且本发明提供的压缩方法程序编写容易，运算效率高、速度快，具有良好的实时性和工程应用价值。
以上为本发明的最佳实施方式，依据本发明公开的内容，本领域的普通技术人员能够显而易见地想到一些雷同、替代方案，均应落入本发明保护的范围。
权利要求
1、一种多视全极化合成孔径雷达数据压缩方法，其特征在于，该方法包括以下步骤
读入多视全极化合成孔径雷达图像每一象素点对应的极化相干矩阵；
根据所述极化相干矩阵对极化总功率进行动态范围可变的压缩；
利用所述动态范围可变压缩的解压值对所述极化相干矩阵的各元素进行归一化处理并压缩。
2、如权利要求1所述多视全极化合成孔径雷达数据压缩方法，其特征在于，所述极化相干矩阵T为3×3非负定厄米特矩阵，如式(5)所示
其中，*代表复共轭，T11、T22、T33均为非负实数，T12、T13、T23为三个复数；则所述极化总功率与极化相干矩阵的关系如式(6)所示
Aspan＝T11+T22+T33(6)
其中，Aspan即为极化总功率。
3、如权利要求2所述多视全极化合成孔径雷达数据压缩方法，其特征在于，所述对极化总功率进行动态范围可变的压缩具体包括
统计所述极化总功率的最大值max(Aspan)与最小值min(Aspan)；
利用式(7)计算压缩所需的参数α和β
使用求得的β值根据式(8)对Aspan值进行调整
根据式(9)使用两字节对调整后的Aspan值进行存储
(9)
其中，
表示不大于logαAspan的最大整数。
4、如权利要求3所述多视全极化合成孔径雷达数据压缩方法，其特征在于，所述对极化相干矩阵各元素进行归一化处理具体包括
利用式(10)计算Aspan的解压值x
利用x值根据式(11)对极化相干矩阵的对角线元素T22和T33进行归一化处理
(11)
利用x值根据式(12)对极化相干矩阵的右上角复元素T12、T13和T23进行实部、虚部分别归一化处理
(12)
其中，real()为取复数实部，imag()为取复数虚部。
5、如权利要求4所述多视全极化合成孔径雷达数据压缩方法，其特征在于，该方法还包括压缩之后的数据存储步骤
将byte1先加上128，再采用无符号8比特整形数据对byte1～byte10进行存储；
采用double型数据对参数α和β进行存储。
6、如权利要求5所述多视全极化合成孔径雷达数据压缩方法，其特征在于，该数据压缩方法对应的数据解压缩方法包括以下步骤
读入压缩后的10字节无符号8比特整形数据，将每字节赋值给byte1～byte10，并对byte1减去128；
读入α和β的值，根据式(10)计算解压值x
根据式(13)解压极化相干矩阵的对角线元素T22和T33
T22＝(byte3+0.5)2/65536·x(13)
T33＝(byte4+0.5)2/65536·x
根据式(14)解压极化相干矩阵的对角线元素T11
T11＝x-T22-T33 (14)
根据式(15)解压极化相干矩阵的右上角复元素T12、T13和T2全文摘要
本发明涉及一种多视全极化合成孔径雷达数据压缩方法，包括以下步骤首先读入多视全极化合成孔径雷达图像每一象素点对应的极化相干矩阵；然后根据极化相干矩阵对极化总功率进行动态范围可变的压缩；再利用动态范围可变压缩的解压值对极化相干矩阵的各元素进行归一化处理并压缩。本发明的技术方案采用极化相干矩阵表示全极化参数，并对极化总功率进行动态范围可变的压缩；不使用极化总功率，而使用它的解压值对极化相干矩阵的元素进行归一化；使该方案在压缩精度、信噪比、保证极化相干矩阵特征值非负和极化特征保持等方面都要优于现有方法，并且运算效率高、速度快，具有良好的实时性和工程应用价值。
文档编号G01S17/00GK101363911SQ20081022294
公开日2009年2月11日 申请日期2008年9月23日 优先权日2008年9月23日
发明者安文韬, 张卫杰, 健 杨 申请人:清华大学