CN109061640B - 一种用于顺轨干涉sar海流反演的方位模糊抑制方法 - Google Patents

一种用于顺轨干涉sar海流反演的方位模糊抑制方法 Download PDF

Info

Publication number
CN109061640B
CN109061640B CN201810707683.2A CN201810707683A CN109061640B CN 109061640 B CN109061640 B CN 109061640B CN 201810707683 A CN201810707683 A CN 201810707683A CN 109061640 B CN109061640 B CN 109061640B
Authority
CN
China
Prior art keywords
doppler
orbit
insar
interference
curve
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810707683.2A
Other languages
English (en)
Other versions
CN109061640A (zh
Inventor
刘保昌
何宜军
李秀仲
端和阳
宋欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN201810707683.2A priority Critical patent/CN109061640B/zh
Publication of CN109061640A publication Critical patent/CN109061640A/zh
Application granted granted Critical
Publication of CN109061640B publication Critical patent/CN109061640B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904—SAR modes

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种用于顺轨干涉SAR海流反演的方位模糊抑制方法，属于海洋遥感技术领域，该方法将待测量海流的海面某区域对应的顺轨InSAR两通道SAR复图像进行配准后变换到距离频率‑多普勒频率域；然后利用距离频率维的独立同分布样本计算顺轨InSAR的协方差矩阵；再根据计算得到的协方差矩阵计算顺轨干涉特征谱熵随多普勒频率变化的曲线；根据上述曲线，分别找到曲线的左半部分和右半部分的极大值点；然后去除曲线左边极大值点以左的所有多普勒单元以及曲线右边极大值点以右的所有多普勒单元；再利用自适应迭代法确定特征谱熵的阈值；去除特征谱熵大于阈值的多普勒单元；最后基于得到的InSAR数据利用干涉法计算海流的距离向速度。

Description

一种用于顺轨干涉SAR海流反演的方位模糊抑制方法
技术领域
本发明属于海洋遥感技术领域，涉及一种针对近岸海域海流反演、适用于顺轨干涉SAR的方位模糊抑制方法。
背景技术
海流是海洋和气象领域中诸多过程中最基本、最重要的要素之一，对海洋中多种生物过程、化学过程和物理过程都有制约作用。海流支配着海水在全球范围内的运动，在整个海洋中输送和混合营养物质、盐、气体、生物和热量。这就意味着了解海流的存在、海流方向和海流速度对商业、社会服务和科学研究均具有重要的意义。例如，石油和天然气勘探等行业需要可靠的海流数据，以确保安全的工作条件；船舶需要依靠海流信息进行航线规划，以尽量减少所使用的燃料或在比赛中获得相对于竞争对手的战术优势。其他应用包括海上搜索和救援、水污染地图绘制和封堵、和全球热量输运。此外，海流信息对于沿海地区的生态环境具有深远的影响，获取高精度、时间连续、高分辨率的海流观测数据可以很好地服务于沿海地区的生态保护和经济发展。比如，近岸水深测量、潮汐图绘制、潮汐发电、近岸平台设施风险评估等。
观测海流比较直接的一种方法是“现场观测”。这种流场测量方法的优点是其测量精度较高；然而，这种方法的缺点是空间覆盖范围有限，且观测成本较高，使得“现场观测”难以满足实际应用的需求。高频地波雷达是一种专门用来观测沿岸海流的雷达设备，可以观测近岸200公里以内的海域。然而高频地波雷达测量的空间范围有限，难以进行全球观测，并且其站点建设和维护在环境较为恶劣的偏远地区非常困难。卫星高度计是另外一种观测海流的手段，利用高度计测量海面的高度以及高度变化可以获得海流的相关信息。然而，利用高度计测流也存在某些缺点。比如，高度计通常仅适用于测量大范围的地转流；存在重复观测周期较长的问题；并且，高度计接收到的雷达回波信号容易受到来自陆地回波信号的“污染”以及大气误差的影响。星载顺轨干涉合成孔径雷达(InSAR)，作为一种先进的主动式微波遥感系统已经广泛应用于海洋环境探测，尤其是海表面流场的测量。InSAR这种方法利用两个沿着轨道放置的天线，通过计算两天线的干涉相位即可获取径向海流信息。相比较而言，InSAR方法具有高分辨、宽覆盖(可以观测全球范围的海流信息)、不受气象条件限制等优势，弥补了传统海流测量方法的不足。
利用星载InSAR观测海流在带来优势的同时也产生了新的问题。其中一个问题就是“方位模糊”问题。该问题产生的原因是雷达的脉冲重复频率(PRF)有限所导致的SAR信号方位向“欠采样”。方位模糊问题在SAR图像域的表现形式为“重影”(Ghost Image)现象，即：某个位置的场景特征会在方位向上移动并“叠加”到另外的场景特征上。通常情况下，方位模糊问题对于“机载”InSAR系统并不构成太大的影响。这是因为对于机载SAR，PRF很容易超过6dB多普勒带宽，从而使得多普勒模糊分量几乎可以忽略。与机载SAR不同的是，星载SAR系统的PRF较难设计得超过6dB多普勒带宽。这是因为：卫星的速度远远大于飞机的速度，使得星载SAR的多普勒带宽较大；雷达脉冲的“占空比”也限制了使用较高的PRF；此外，为了获得较大的距离向刈幅宽度，PRF也不宜取得过大。上述原因造成了星载InSAR系统的方位模糊度较高。
尽管星载SAR的方位模糊度较高，对于均匀散射的“开阔洋”区域，方位模糊对流场反演的影响并不大。然而，对于“非均匀”的“近岸”区域，方位模糊对星载InSAR测流性能的影响非常大。原因可说明如下：其一，海岸(陆地)的雷达后向散射系数通常远大于海面(海水)的散射系数，因此，来自陆地的方位模糊分量(即其重影)会叠加到海面上，从而造成InSAR测流精度的严重下降；其二，陆地重影的速度与海表速度是不相同的(陆地表面的速度通常为零)，这也会改变海表流场的测量值。到目前为止，抑制方位模糊的方法主要分为两类：SAR图像域方法与多普勒频率域方法。对于第一类方法，现有文献给出了一种方位模糊抑制方法，该方法工作在SAR图像域，根据理论计算“重影”偏离其真实位置的距离，然后据此丢弃掉某些像素。然而，此方法存在一个明显的缺点：由于在SAR图像域丢弃了某些像素，因此，某些区域的流场信息会丢失。另外，该方法的精度也比较低，无法获得高精度的海流速度。现有文献给出了另外一种抑制方位模糊的方法。该方法属于第二类方法，工作在多普勒频率域。此方法是一种通道均衡算法，主要应用于“陆地”场景的地面运动目标检测。然而，当它应用于海面场景时，会产生额外的“虚假”干涉相位，使得该方法并不适合应用于海流反演。对于多普勒频率域方法(第二类方法)，现有文献提出了一种方法，该方法能够在一定程度上抑制方位模糊，不过该方法主要是用来估计顺轨InSAR系统的有效基线，且该方法仅适用于静止的均匀陆地场景。因此，该方法并不适合于海面场景，无法应用于流场反演。综上所述，到目前为止，还没有一种针对“近岸”海域海流反演的、有效的顺轨InSAR方位模糊抑制算法。因此，发明一种能够业务化运行的方位模糊抑制算法对于海面动力参数反演，尤其是近岸海流反演将具有重要的意义，能够有效地服务于沿岸生态保护和经济发展需求。
发明内容
为了解决现有技术中存在的问题，本发明的目的是提供一种针对“近岸”海域海流反演的、有效且能够实现业务化运行的顺轨InSAR方位模糊抑制算法。
为了达到上述目的，本发明提出的技术方案为：一种用于顺轨干涉SAR海流反演的方位模糊抑制方法，包括以下步骤：
步骤一、将待测量海流的海面某区域对应的顺轨InSAR两通道SAR复图像进行配准后变换到距离频率-多普勒频率域；然后利用距离频率维的独立同分布样本计算顺轨InSAR的协方差矩阵；再根据计算得到的协方差矩阵计算顺轨干涉特征谱熵随多普勒频率变化的曲线；
步骤二、根据计算得到的特征谱熵随多普勒频率变化的曲线，分别找到曲线的左半部分和右半部分的极大值点；然后去除曲线左边极大值点以左的所有多普勒单元以及曲线右边极大值点以右的所有多普勒单元；
步骤三、利用自适应迭代法确定特征谱熵的阈值；去除特征谱熵大于阈值的多普勒单元；
步骤四、基于步骤二和步骤三处理后的InSAR数据利用干涉法计算海流的距离向速度。
对上述技术方案的进一步设计为：所述步骤一的具体步骤如下：
步骤1、对InSAR两通道SAR复图像进行配准后利用一个矩形窗提取待反演海流的区域，然后把该区域对应的InSAR数据变换到距离频率-多普勒频率域，分别得到两通道数据信号：
Figure BDA0001715885930000031
与
Figure BDA0001715885930000032
其中，
Figure BDA0001715885930000033
(k＝1,2,…,K)表示第k个距离频率单元对应的距离频率，K为距离频率带宽内总的距离频率单元数；
Figure BDA0001715885930000034
(m＝1,2,…,M)表示第m个多普勒单元对应的多普勒频率，M为多普勒基带内总的多普勒单元数；
步骤2、根据二维频率域两通道InSAR数据
Figure BDA0001715885930000035
与
Figure BDA0001715885930000036
计算顺轨干涉特征谱熵
Figure BDA0001715885930000037
步骤2-1、基于
Figure BDA0001715885930000038
和
Figure BDA0001715885930000039
由下式计算“顺轨干涉协方差矩阵”
Figure BDA00017158859300000310
Figure BDA00017158859300000311
上式中，(·)*表示取复共轭，距离频率维用来提供估计协方差矩阵的独立同分布样本，协方差矩阵
Figure BDA00017158859300000312
为多普勒频率
Figure BDA00017158859300000313
的函数；
步骤2-2、计算协方差矩阵
Figure BDA00017158859300000314
的两个特征值：
Figure BDA00017158859300000315
与
Figure BDA00017158859300000316
特征值
Figure BDA00017158859300000317
与
Figure BDA00017158859300000318
满足如下关系：
Figure BDA00017158859300000319
上式中，P为
Figure BDA00017158859300000320
的特征矩阵，P-1表示P的逆矩阵。
步骤2-3、基于
Figure BDA00017158859300000321
与
Figure BDA00017158859300000322
由以下三式计算顺轨干涉“特征谱熵”
Figure BDA00017158859300000323
Figure BDA00017158859300000324
Figure BDA0001715885930000041
Figure BDA0001715885930000042
所述步骤二的具体步骤如下：
步骤3、搜索顺轨干涉特征谱熵曲线
Figure BDA0001715885930000043
分别找到该曲线左半部分和右半部分的极大值点：
Figure BDA0001715885930000044
与
Figure BDA0001715885930000045
Figure BDA0001715885930000046
与
Figure BDA0001715885930000047
满足如下两式：
Figure BDA0001715885930000048
Figure BDA0001715885930000049
上式中，
Figure BDA00017158859300000410
表示在多普勒范围
Figure BDA00017158859300000411
内取
Figure BDA00017158859300000412
的极大值；
Figure BDA00017158859300000413
表示在多普勒范围
Figure BDA00017158859300000414
内取
Figure BDA00017158859300000415
的极大值；
Figure BDA00017158859300000416
表示取函数
Figure BDA00017158859300000417
的自变量；fPRF表示雷达的脉冲重复频率(PRF)；
步骤4、丢弃掉特征谱熵曲线
Figure BDA00017158859300000418
左边极大值点
Figure BDA00017158859300000419
以左的所有多普勒单元以及特征谱熵曲线右边极大值点
Figure BDA00017158859300000420
以右的所有多普勒单元，得到多普勒子带
Figure BDA00017158859300000421
如下式所示：
Figure BDA00017158859300000422
上式中，
Figure BDA00017158859300000423
为多普勒子带
Figure BDA00017158859300000424
中第n个多普勒单元对应的多普勒频率，N为
Figure BDA00017158859300000425
中多普勒单元的总的数目。
所述步骤三的具体步骤如下：
步骤5、挑选位于多普勒子带
Figure BDA00017158859300000426
内的InSAR两通道数据信号，得到：
Figure BDA00017158859300000427
与
Figure BDA00017158859300000428
具体表达式如下：
Figure BDA00017158859300000429
Figure BDA00017158859300000430
步骤6、利用自适应迭代法确定顺轨干涉“特征谱熵”的阈值Hc；
步骤6-1、根据下式计算定义于多普勒子带
Figure BDA0001715885930000051
上的“特征谱熵”
Figure BDA0001715885930000052
Figure BDA0001715885930000053
步骤6-2、基于
Figure BDA0001715885930000054
与
Figure BDA0001715885930000055
利用下式计算定义于多普勒子带
Figure BDA0001715885930000056
上的干涉相位
Figure BDA0001715885930000057
Figure BDA0001715885930000058
上式中，∠{·}表示取一个复数的相位；
步骤6-3、基于
Figure BDA0001715885930000059
与
Figure BDA00017158859300000510
利用下式计算二维频率域的平均相干系数
Figure BDA00017158859300000511
Figure BDA00017158859300000512
步骤6-4、基于上一步计算得到的平均相干系数
Figure BDA00017158859300000513
根据下式计算得到一个“相位随机波动值”φCPF：
Figure BDA00017158859300000514
步骤6-5、设ε为一个变量(0≤ε≤1)，根据此变量以及
Figure BDA00017158859300000515
确定如下一个多普勒子带FD(ε)：
Figure BDA00017158859300000516
上式的含义为：挑选出“特征谱熵”的值小于ε的所有多普勒单元，然后将这些多普勒单元组成多普勒子带FD，该多普勒子带为变量ε的函数，设ε的初始值为1；
步骤6-6、根据下式，从集合
Figure BDA00017158859300000517
中挑选出定义于多普勒子带FD(ε)上的相位集合，得到相位集合
Figure BDA00017158859300000518
Figure BDA00017158859300000519
上式中，I为集合
Figure BDA00017158859300000520
中元素的数目；
步骤6-7、根据下式计算相位集合
Figure BDA00017158859300000521
的“干涉相位的平均变化值”φIPV(ε)：
Figure BDA0001715885930000061
上式中，干涉相位的平均变化值φIPV为参数ε的函数；
步骤6-8、判断条件φIPV(ε)＜φCPF是否成立。若该条件不成立，则按照一定的步长ε0(ε0＞0)减小ε的值，使其变为：
ε＝ε-ε0
利用更新后的ε值重复步骤6-5至步骤6-7，直到条件φIPV(ε)＜φCPF成立；若该条件成立，则执行下一步骤；
步骤6-9、当条件φIPV(ε)＜φCPF首次成立时，将此时的ε值确定为顺轨干涉特征谱熵的阈值Hc，即令Hc＝ε。
步骤7、丢弃掉特征谱熵的值大于阈值Hc的所有多普勒单元，得到多普勒子带
Figure BDA0001715885930000062
如下式所示：
Figure BDA0001715885930000063
上式中，
Figure BDA0001715885930000064
表示多普勒“子带”
Figure BDA0001715885930000065
中第l个(l＝1,2,…,L)多普勒单元对应的多普勒频率(其中，L为
Figure BDA0001715885930000066
中多普勒单元总的数目)。
所述步骤四的具体步骤如下：
步骤8、挑选位于多普勒子带
Figure BDA0001715885930000067
内的InSAR两通道数据信号，得到
Figure BDA0001715885930000068
与
Figure BDA0001715885930000069
具体表达式如下：
Figure BDA00017158859300000610
Figure BDA00017158859300000611
步骤9、基于
Figure BDA00017158859300000612
与
Figure BDA00017158859300000613
利用传统的顺轨干涉法计算距离向海流速度
Figure BDA00017158859300000614
Figure BDA00017158859300000615
上式中，λ为雷达波长，vs为雷达平台的有效速度，deff为InSAR有效基线长度，θinc为入射角。
本发明的有益效果为：
1、相比于传统的不考虑“方位模糊抑制”的顺轨InSAR测流方法，本发明提出“方位模糊抑制”算法与顺轨InSAR测流技术结合后，可以明显地提高“近岸”海域的流场反演精度。
2、相比于SAR图像域的“方位模糊抑制”算法，本发明提出“方位模糊抑制”算法不存在海面场景中某些区域海流信息丢失的缺点。
3、本发明提出的“方位模糊抑制”算法具有“自适应性”，即：能够“自动”地应对不同的雷达参数(比如：信噪比、方位模糊度等)以及不同的场景参数(比如海流速度、海岸陆地散射系数与海水散射系数之比、海浪波长等)。
4、除了两通道InSAR原始数据外，本发明提出“方位模糊抑制”算法需要的附加输入参数较少，因此便于“业务化”运行。
5、本发明方法属于海洋动力环境参数遥感技术领域，可直接服务于沿岸地区的科学研究、经济发展、以及生态保护等。
附图说明
图1为本发明实施例的算法流程图；
图2为近岸场景InSAR干涉对幅度图像
图3为近岸场景InSAR干涉相位图像；
图4为顺轨干涉协方差矩阵的两个特征值随多普勒频率变化的曲线；
图5为顺轨干涉特征谱熵随多普勒频率变化的曲线；
图6为自适应迭代法求解顺轨干涉特征谱熵的阈值的流程图；
图7为干涉相位平均变化值φIPV随参数ε变化的曲线、相位随机波动值φCPF、以及特征谱熵的阈值；
图8为方位模糊抑制后的InSAR干涉对幅度图像；
图9为方位模糊抑制后的InSAR干涉相位图像。
具体实施方式
下面结合附图以及具体实施例对本发明进行详细说明。
实施例
下面结合附图和实例进一步说明本发明的具体实施方式。
本发明提出的一种针对近岸海域海流反演、适用于顺轨干涉SAR的方位模糊抑制方法的总体流程图如图1所示，具体步骤如下：
1)利用常用的SAR成像算法(比如：距离-多普勒算法、Omega-K算法等)对顺轨InSAR两通道数据进行“成像”(即距离压缩和方位压缩)。图2给出了InSAR原始数据成像后的干涉对幅度图像。该SAR图像描述了一个“近岸”场景，其中既包括海岸陆地场景，也包括海面场景。从图2中可以看到海岸上的三个“目标”的“重影”叠加到了海面区域。图3给出了InSAR干涉相位图，从中可以看出在近岸海域，由于受到陆地“重影”的影响，海面的干涉相位受到了严重的影响。
2)对InSAR两通道SAR复图像进行配准，消除由顺轨基线所造成的两幅SAR复图像在“方位向”的偏移。
3)利用一个矩形窗提取待反演海流的区域(见图2中的矩形窗)，然后把该区域对应的InSAR数据变换到二维频率域(即距离频率/多普勒频率域)，分别得到两通道数据信号：
Figure BDA0001715885930000081
与
Figure BDA0001715885930000082
其中，
Figure BDA0001715885930000083
(k＝1,2,…,K)表示第k个距离频率单元对应的距离频率(K为距离频率带宽内总的距离频率单元数)；
Figure BDA0001715885930000084
(m＝1,2,…,M)表示第m个多普勒单元对应的多普勒频率(M为多普勒基带内总的多普勒单元数)。
4)根据二维频率域两通道InSAR数据
Figure BDA0001715885930000085
与
Figure BDA0001715885930000086
计算顺轨干涉“特征谱熵”
Figure BDA0001715885930000087
具体步骤如下：
Step 1基于
Figure BDA0001715885930000088
和
Figure BDA0001715885930000089
由下式计算“顺轨干涉协方差矩阵”
Figure BDA00017158859300000810
Figure BDA00017158859300000811
上式中，(·)*表示取复共轭，距离频率维用来提供估计协方差矩阵的独立同分布样本，协方差矩阵
Figure BDA00017158859300000812
为多普勒频率
Figure BDA00017158859300000813
的函数。
Step 2计算协方差矩阵
Figure BDA00017158859300000814
的两个特征值：
Figure BDA00017158859300000815
与
Figure BDA00017158859300000816
特征值
Figure BDA00017158859300000817
与
Figure BDA00017158859300000818
满足如下关系：
Figure BDA00017158859300000819
上式中，P为
Figure BDA00017158859300000820
的特征矩阵，P-1表示P的逆矩阵。图4给出了
Figure BDA00017158859300000821
与
Figure BDA00017158859300000822
随多普勒频率
Figure BDA00017158859300000823
变化的曲线。
Step 3基于
Figure BDA00017158859300000824
与
Figure BDA00017158859300000825
由以下三式计算顺轨干涉“特征谱熵”
Figure BDA00017158859300000826
Figure BDA00017158859300000827
Figure BDA0001715885930000091
Figure BDA0001715885930000092
图5给出了顺轨干涉“特征谱熵”
Figure BDA0001715885930000093
随多普勒频率
Figure BDA0001715885930000094
变化的曲线。顺轨干涉“特征谱熵”是本发明“方位模糊抑制算法”的核心参数，其物理意义为：该参数表征在某个多普勒单元内，方位模糊分量与非模糊分量之间“混合”的程度；该参数的值越大，前述两个分量“混合”的程度越高。基于“特征谱熵”及其物理意义，可在多普勒域实现“方位模糊”分量(对应于SAR图像中的“重影”)与非模糊信号分量之间的“分离”。
5)搜索顺轨干涉“特征谱熵”曲线
Figure BDA0001715885930000095
分别找到该曲线“左半部分”和“右半部分”的两个“极大值点”：
Figure BDA0001715885930000096
与
Figure BDA0001715885930000097
Figure BDA0001715885930000098
与
Figure BDA0001715885930000099
满足如下两式：
Figure BDA00017158859300000910
Figure BDA00017158859300000911
上式中，
Figure BDA00017158859300000912
表示在多普勒范围
Figure BDA00017158859300000913
内取
Figure BDA00017158859300000914
的极大值；
Figure BDA00017158859300000915
表示在多普勒范围
Figure BDA00017158859300000916
内取
Figure BDA00017158859300000917
的极大值；
Figure BDA00017158859300000918
表示取函数
Figure BDA00017158859300000919
的自变量；fPRF表示雷达的脉冲重复频率(PRF)。图5给出了
Figure BDA00017158859300000920
的两个极大值点
Figure BDA00017158859300000921
与
Figure BDA00017158859300000922
6)丢弃掉“特征谱熵”曲线
Figure BDA00017158859300000923
左边“极大值点”
Figure BDA00017158859300000924
以左的所有多普勒单元以及“特征谱熵”曲线右边“极大值点”
Figure BDA00017158859300000925
以右的所有多普勒单元，得到多普勒“子带”
Figure BDA00017158859300000926
(见图5)，如下式所示：
Figure BDA00017158859300000927
上式中，
Figure BDA00017158859300000928
为多普勒子带
Figure BDA00017158859300000929
中第n个多普勒单元对应的多普勒频率，N为
Figure BDA00017158859300000930
中多普勒单元的总的数目。
7)挑选位于多普勒子带
Figure BDA00017158859300000931
(见图5)内的InSAR两通道数据信号，得到：
Figure BDA00017158859300000932
与
Figure BDA0001715885930000101
具体表达式如下：
Figure BDA0001715885930000102
Figure BDA0001715885930000103
此步骤的作用是剔除掉方位模糊分量占支配地位的那些多普勒单元。
8)利用“自适应迭代法”确定顺轨干涉“特征谱熵”的阈值Hc。图6给出了“自适应迭代法”的流程图。结合图6，具体步骤描述如下：
Step 1基于上一步得到的
Figure BDA0001715885930000104
与
Figure BDA0001715885930000105
利用式(1)-(5)计算定义于多普勒子带
Figure BDA0001715885930000106
上的“特征谱熵”
Figure BDA0001715885930000107
Figure BDA0001715885930000108
与
Figure BDA0001715885930000109
满足如下关系：
Figure BDA00017158859300001010
Step 2基于
Figure BDA00017158859300001011
与
Figure BDA00017158859300001012
利用下式计算定义于多普勒子带
Figure BDA00017158859300001013
上的干涉相位
Figure BDA00017158859300001014
Figure BDA00017158859300001015
上式中，∠{·}表示取一个复数的相位。
Step 3基于
Figure BDA00017158859300001016
与
Figure BDA00017158859300001017
利用下式计算二维频率域的平均相干系数
Figure BDA00017158859300001018
Figure BDA00017158859300001019
Step 4基于上一步计算得到的平均相干系数
Figure BDA00017158859300001020
根据下式计算得到一个“相位随机波动值”φCPF：
Figure BDA00017158859300001021
上式中的“相位随机波动值”φCPF在物理意义上为干涉相位
Figure BDA00017158859300001022
的“克拉美罗界”。
Step 5设ε为一个变量(0≤ε≤1)，根据此变量以及
Figure BDA00017158859300001023
确定如下一个多普勒“子带”FD(ε)：
Figure BDA0001715885930000111
上式的含义为：挑选出“特征谱熵”的值小于ε的所有多普勒单元，然后将这些多普勒单元组成多普勒子带FD，该多普勒子带为变量ε的函数。设ε的初始值为1。
Step 6根据下式，从集合
Figure BDA0001715885930000112
中挑选出定义于多普勒子带FD(ε)上的相位集合，得到相位集合
Figure BDA0001715885930000113
Figure BDA0001715885930000114
上式中，I为集合
Figure BDA0001715885930000115
中元素的数目。
Step 7根据下式计算相位集合
Figure BDA0001715885930000116
的“干涉相位的平均变化值”φIPV(ε)：
Figure BDA0001715885930000117
上式中，“干涉相位的平均变化值”φIPV为参数ε的函数。
Step 8判断条件φIPV(ε)＜φCPF是否成立。若该条件不成立，则按照一定的步长ε0(ε0＞0)减小ε的值，使其变为：
ε＝ε-ε0 (18)
利用更新后的ε值重复步骤Step 5—Step 7，直到条件φIPV(ε)＜φCPF成立。若该条件成立，则执行下一步骤。在本实施例中，步长ε0的值取为0.01。
Step 9当条件φIPV(ε)＜φCPF“首次”成立时，将此时的ε值确定为顺轨干涉“特征谱熵”的阈值Hc(即令Hc＝ε)。图7给出了“干涉相位的平均变化值”φIPV随参数ε变化的曲线、“相位随机波动值”φCPF、以及本实施例中所确定的“特征谱熵”的“阈值”Hc。
9)丢弃掉“特征谱熵”的值大于阈值Hc的所有多普勒单元，得到多普勒“子带”
Figure BDA0001715885930000118
(见图5)，如下式所示：
Figure BDA0001715885930000119
上式中，
Figure BDA00017158859300001110
表示多普勒“子带”
Figure BDA00017158859300001112
中第l个(l＝1,2,…,L)多普勒单元对应的多普勒频率(其中，L为
Figure BDA00017158859300001111
中多普勒单元总的数目)。
10)挑选位于多普勒子带
Figure BDA0001715885930000121
内的InSAR两通道数据信号，得到
Figure BDA0001715885930000122
与
Figure BDA0001715885930000123
具体表达式如下：
Figure BDA0001715885930000124
Figure BDA0001715885930000125
图8给出了利用本发明提出的算法进行方位模糊抑制后的InSAR干涉对幅度图，而图9给出了利用本发明提出的算法进行方位模糊抑制后的InSAR干涉相位图。对比图2和图8，可以发现：利用本发明提出的算法进行方位模糊抑制后，海岸上的三个目标的“重影”在海面上消失了。此外，对比图3和图9，可以发现：利用本发明提出的算法进行方位模糊抑制后，近岸海面的干涉相位发生了明显的变化，变得与其它海域的相位一致了，从而表明了本发明“方位模糊抑制”算法的有效性。
11)基于
Figure BDA0001715885930000126
与
Figure BDA0001715885930000127
利用传统的顺轨干涉法计算距离向海流速度
Figure BDA0001715885930000128
Figure BDA0001715885930000129
上式中，λ为雷达波长，vs为雷达平台的有效速度，deff为InSAR有效基线长度，θinc为入射角。为了验证本发明提出的“方位模糊抑制”算法对流场反演精度的提高效果，对比了传统顺轨干涉法(无方位模糊抑制步骤)的海流反演结果和加入本发明算法后的干涉法海流反演结果。对于本实施例，以图8中“矩形窗”内的海面区域为例，传统的顺轨干涉法海流反演的误差为3.1m/s，而加入本发明提出的“方位模糊抑制”算法后，海流的反演误差减小为0.02m/s。上述结果进一步验证了本发明提出的“方位模糊抑制”算法的有效性。
本发明的方法不仅适用于近岸海域海表流场反演，也适用于河流表面流速的反演。
本发明的技术方案不局限于上述实施例，凡采用等同替换方式得到的技术方案均落在本发明要求保护的范围内。

Claims (3)

1.一种用于顺轨干涉SAR海流反演的方位模糊抑制方法，其特征在于，包括以下步骤：
步骤一、将待测量海流的海面某区域对应的顺轨InSAR两通道SAR复图像进行配准后变换到距离频率-多普勒频率域；然后利用距离频率维的独立同分布样本计算顺轨InSAR的协方差矩阵；再根据计算得到的协方差矩阵计算顺轨干涉特征谱熵随多普勒频率变化的曲线；
步骤二、根据计算得到的特征谱熵随多普勒频率变化的曲线，分别找到曲线的左半部分和右半部分的极大值点；然后去除曲线左边极大值点以左的所有多普勒单元以及曲线右边极大值点以右的所有多普勒单元；
步骤三、利用自适应迭代法确定特征谱熵的阈值；去除特征谱熵大于阈值的多普勒单元；具体方法如下：
步骤5、挑选位于多普勒子带
Figure FDA0003634385850000011
内的InSAR两通道数据信号，得到：
Figure FDA0003634385850000012
与
Figure FDA0003634385850000013
具体表达式如下：
Figure FDA0003634385850000014
Figure FDA0003634385850000015
其中，
Figure FDA0003634385850000016
表示第k个距离频率单元对应的距离频率，K为距离频率带宽内总的距离频率单元数；
Figure FDA0003634385850000017
表示第m个多普勒单元对应的多普勒频率，M为多普勒基带内总的多普勒单元数；
Figure FDA0003634385850000018
为多普勒子带
Figure FDA0003634385850000019
中第n个多普勒单元对应的多普勒频率，N为
Figure FDA00036343858500000110
中多普勒单元的总的数目；
步骤6、利用自适应迭代法确定顺轨干涉“特征谱熵”的阈值Hc；
步骤6-1、根据下式计算定义于多普勒子带
Figure FDA00036343858500000111
上的“特征谱熵”
Figure FDA00036343858500000112
Figure FDA00036343858500000113
步骤6-2、基于
Figure FDA00036343858500000114
与
Figure FDA00036343858500000115
利用下式计算定义于多普勒子带
Figure FDA00036343858500000116
上的干涉相位
Figure FDA00036343858500000117
Figure FDA00036343858500000118
上式中，∠{·}表示取一个复数的相位；
步骤6-3、基于
Figure FDA0003634385850000021
与
Figure FDA0003634385850000022
利用下式计算二维频率域的平均相干系数
Figure FDA0003634385850000023
Figure FDA0003634385850000024
步骤6-4、基于上一步计算得到的平均相干系数
Figure FDA0003634385850000025
根据下式计算得到一个“相位随机波动值”φCPF：
Figure FDA0003634385850000026
步骤6-5、设ε为一个变量(0≤ε≤1)，根据此变量以及
Figure FDA0003634385850000027
确定如下一个多普勒子带FD(ε)：
Figure FDA0003634385850000028
上式的含义为：挑选出“特征谱熵”的值小于ε的所有多普勒单元，然后将这些多普勒单元组成多普勒子带FD，该多普勒子带为变量ε的函数，设ε的初始值为1；
步骤6-6、根据下式，从集合
Figure FDA0003634385850000029
中挑选出定义于多普勒子带FD(ε)上的相位集合，得到相位集合
Figure FDA00036343858500000210
Figure FDA00036343858500000211
上式中，I为集合
Figure FDA00036343858500000212
中元素的数目；
步骤6-7、根据下式计算相位集合
Figure FDA00036343858500000213
的“干涉相位的平均变化值”φIPV(ε)：
Figure FDA00036343858500000214
上式中，干涉相位的平均变化值φIPV为参数ε的函数；
步骤6-8、判断条件φIPV(ε)＜φCPF是否成立，若该条件不成立，则按照一定的步长ε0(ε0＞0)减小ε的值，使其变为：
ε＝ε-ε0
利用更新后的ε值重复步骤6-5至步骤6-7，直到条件φIPV(ε)＜φCPF成立；若该条件成立，则执行下一步骤；
步骤6-9、当条件φIPV(ε)＜φCPF首次成立时，将此时的ε值确定为顺轨干涉特征谱熵的阈值Hc，即令Hc＝ε；
步骤7、丢弃掉特征谱熵的值大于阈值Hc的所有多普勒单元，得到多普勒子带
Figure FDA0003634385850000031
如下式所示：
Figure FDA0003634385850000032
上式中，
Figure FDA0003634385850000033
表示多普勒“子带”
Figure FDA0003634385850000034
中第l个(l＝1,2,…,L)多普勒单元对应的多普勒频率，其中，L为
Figure FDA0003634385850000035
中多普勒单元总的数目；
步骤四、基于步骤二和步骤三处理后的InSAR数据利用干涉法计算海流的距离向速度；
具体方法如下：
步骤8、挑选位于多普勒子带
Figure FDA0003634385850000036
内的InSAR两通道数据信号，得到
Figure FDA0003634385850000037
与
Figure FDA0003634385850000038
具体表达式如下：
Figure FDA0003634385850000039
Figure FDA00036343858500000310
步骤9、基于
Figure FDA00036343858500000311
与
Figure FDA00036343858500000312
利用传统的顺轨干涉法计算距离向海流速度
Figure FDA00036343858500000313
Figure FDA00036343858500000314
上式中，λ为雷达波长，vs为雷达平台的有效速度，deff为InSAR有效基线长度，θinc为入射角。
2.根据权利要求1所述的用于顺轨干涉SAR海流反演的方位模糊抑制方法，其特征在于：所述步骤一的具体步骤如下：
步骤1、对InSAR两通道SAR复图像进行配准后利用一个矩形窗提取待反演海流的区域，然后把该区域对应的InSAR数据变换到距离频率-多普勒频率域，分别得到两通道数据信号：
Figure FDA00036343858500000315
与
Figure FDA00036343858500000316
步骤2、根据二维频率域两通道InSAR数据
Figure FDA0003634385850000041
与
Figure FDA0003634385850000042
计算顺轨干涉特征谱熵
Figure FDA0003634385850000043
步骤2-1、基于
Figure FDA0003634385850000044
和
Figure FDA0003634385850000045
由下式计算“顺轨干涉协方差矩阵”
Figure FDA0003634385850000046
Figure FDA0003634385850000047
上式中，(·)*表示取复共轭，距离频率维用来提供估计协方差矩阵的独立同分布样本，协方差矩阵
Figure FDA0003634385850000048
为多普勒频率
Figure FDA0003634385850000049
的函数；
步骤2-2、计算协方差矩阵
Figure FDA00036343858500000410
的两个特征值：
Figure FDA00036343858500000411
与
Figure FDA00036343858500000412
特征值
Figure FDA00036343858500000413
与
Figure FDA00036343858500000414
满足如下关系：
Figure FDA00036343858500000415
上式中，P为
Figure FDA00036343858500000416
的特征矩阵，P-1表示P的逆矩阵；
步骤2-3、基于
Figure FDA00036343858500000417
与
Figure FDA00036343858500000418
由以下三式计算顺轨干涉“特征谱熵”
Figure FDA00036343858500000419
Figure FDA00036343858500000420
Figure FDA00036343858500000421
Figure FDA00036343858500000422
3.根据权利要求2所述的用于顺轨干涉SAR海流反演的方位模糊抑制方法，其特征在于：所述步骤二的具体步骤如下：
步骤3、搜索顺轨干涉特征谱熵曲线
Figure FDA00036343858500000423
分别找到该曲线左半部分和右半部分的极大值点：
Figure FDA00036343858500000424
与
Figure FDA00036343858500000425
与
Figure FDA00036343858500000426
满足如下两式：
Figure FDA0003634385850000051
Figure FDA0003634385850000052
上式中，
Figure FDA0003634385850000053
表示在多普勒范围
Figure FDA0003634385850000054
内取
Figure FDA0003634385850000055
的极大值；
Figure FDA0003634385850000056
表示在多普勒范围
Figure FDA0003634385850000057
内取
Figure FDA0003634385850000058
的极大值；
Figure FDA0003634385850000059
表示取函数
Figure FDA00036343858500000510
的自变量；fPRF表示雷达的脉冲重复频率(PRF)；
步骤4、丢弃掉特征谱熵曲线
Figure FDA00036343858500000511
左边极大值点
Figure FDA00036343858500000512
以左的所有多普勒单元以及特征谱熵曲线右边极大值点
Figure FDA00036343858500000513
以右的所有多普勒单元，得到多普勒子带
Figure FDA00036343858500000514
如下式所示：
Figure FDA00036343858500000515
CN201810707683.2A 2018-07-02 2018-07-02 一种用于顺轨干涉sar海流反演的方位模糊抑制方法 Active CN109061640B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810707683.2A CN109061640B (zh) 2018-07-02 2018-07-02 一种用于顺轨干涉sar海流反演的方位模糊抑制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810707683.2A CN109061640B (zh) 2018-07-02 2018-07-02 一种用于顺轨干涉sar海流反演的方位模糊抑制方法

Publications (2)

Publication Number Publication Date
CN109061640A CN109061640A (zh) 2018-12-21
CN109061640B true CN109061640B (zh) 2022-06-21

Family

ID=64818255

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810707683.2A Active CN109061640B (zh) 2018-07-02 2018-07-02 一种用于顺轨干涉sar海流反演的方位模糊抑制方法

Country Status (1)

Country Link
CN (1) CN109061640B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110501708B (zh) * 2019-08-29 2021-03-30 北京航空航天大学 一种多通道星载topsar方位模糊度分析方法
CN110554377B (zh) * 2019-09-05 2021-04-09 中国科学院电子学研究所 基于多普勒中心偏移的单通道sar二维流场反演方法及系统
CN111398911B (zh) * 2020-03-24 2022-03-08 中国人民解放军海军航空大学 Mimo雷达目标检测方法及装置
CN112862868B (zh) * 2021-01-31 2023-12-01 南京信息工程大学 一种基于线性变换和小波分析的运动海浪图像配准融合方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104360325A (zh) * 2014-11-26 2015-02-18 西安电子科技大学 机载前视阵雷达的空时自适应处理方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5745069A (en) * 1996-09-10 1998-04-28 Ball Corporation Reduction of radar antenna area
JP5419749B2 (ja) * 2009-06-03 2014-02-19 キヤノン株式会社 画像処理装置、画像処理方法およびプログラム
US8493262B2 (en) * 2011-02-11 2013-07-23 Mitsubishi Electric Research Laboratories, Inc. Synthetic aperture radar image formation system and method
CN102520404B (zh) * 2011-11-30 2013-07-03 北京理工大学 一种基于图像质量最优的sar多普勒模糊数估计方法
CN103413146B (zh) * 2013-08-23 2017-03-29 西安电子科技大学 基于Freeman熵和自学习的极化SAR图像精细分类方法
CN104698431B (zh) * 2015-03-17 2017-06-30 河海大学 基于模糊分量doa估计的多通道sar方位解模糊方法
CN105445711B (zh) * 2015-11-27 2017-08-01 南京信息工程大学 一种基于逆Omega‑K算法的海面要素SAR原始数据仿真方法
CN107610160B (zh) * 2017-09-21 2020-05-08 中南大学 一种极化sar船舶速度估计方法及系统
CN108152817A (zh) * 2017-11-20 2018-06-12 上海无线电设备研究所 流速测量中顺轨干涉sar的多普勒中心误差补偿方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104360325A (zh) * 2014-11-26 2015-02-18 西安电子科技大学 机载前视阵雷达的空时自适应处理方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于DDC模型的机载雷达多普勒参数估计研究;许稼等;《电子学报》;20040925(第09期);1421-1424 *

Also Published As

Publication number Publication date
CN109061640A (zh) 2018-12-21

Similar Documents

Publication Publication Date Title
CN109061640B (zh) 一种用于顺轨干涉sar海流反演的方位模糊抑制方法
CN109856635B (zh) 一种csar地面动目标重聚焦成像方法
Gogineni et al. Bed topography of Jakobshavn Isbræ, Greenland, and Byrd Glacier, Antarctica
CN110109102B (zh) 一种sar运动目标检测与速度估计的方法
CN110632594B (zh) 一种长波长星载sar成像方法
CN110568434B (zh) 一种多通道匀加速sar动目标二维速度估计方法
CN103529437A (zh) 系留气球载相控阵雷达在多目标下分辨空地目标的方法
CN108776342A (zh) 一种高速平台sar慢速动目标检测与速度估计方法
CN110261833A (zh) 高分辨星载sar成像误差估计与补偿方法
CN108107432B (zh) 基于时域扰动的高低轨双基sar保相成像方法
Sletten et al. Maritime signature correction with the NRL multichannel SAR
CN102565772B (zh) 基于sar子孔径序列图像的海洋动态信息提取方法
CN103091682B (zh) 基于时频分析InISAR多动目标成像和运动轨迹重建法
CN110133646B (zh) 基于nlcs成像的双基前视sar的多通道两脉冲杂波对消方法
Suchandt et al. X-band sea surface coherence time inferred from bistatic SAR interferometry
Kraus et al. Topography correction for distributed SAR imaging: A case study for TerraSAR-X and TanDEM-X
CN112162282A (zh) 一种基于合成孔径雷达海表流速反演方法
Sun et al. Remote sensing of ocean surface wind direction with shipborne high frequency surface wave radar
Liu et al. A new azimuth ambiguity suppression algorithm for surface current measurement in coastal waters and rivers with along-track InSAR
CN113204020B (zh) 一种基于频谱分割的海浪波谱仪斑点噪声谱估计方法
CN107422320A (zh) 一种消除降雨对x波段雷达观测海浪的影响的方法
Wu et al. Statistical analysis of monopulse‐synthetic aperture radar for constant false‐alarm rate detection of ground moving targets
Sun et al. Ocean waves inversion based on airborne radar images with small incident angle
Suchandt et al. First results of tandem-x along-track interferometry
Pang et al. A GMTI algorithm based on novel channel equalization method for multi-channel imaging radar

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant