一种基于分布式压缩感知的星载差分层析SAR成像方法

一种基于分布式压缩感知的星载差分层析sar成像方法
技术领域
1.本发明属于高维合成孔径雷达成像处理技术领域，具体涉及一种基于分布式 压缩感知的星载差分层析sar成像方法。


背景技术：

2.合成孔径雷达(synthetic aperture radar,sar)是一种能够全天时、全天候工 作的微波成像雷达，对地表具有一定的穿透能力。基于sar合成孔径实现方位 向高分辨率，不改变现有sar成像体制，层析sar(tomographic sar,tomosar) 利用观测时间相同的多幅sar图像沿垂直于雷达视线与飞行轨迹的高程向合成 孔径，可对复杂地形高分辨率三维成像，因而也称为3d sar。
3.近年来，层析sar被广泛应用于建设领域。然而，3d sar对雷达天线数 目与空间分布要求苛刻，这大大增加了3d sar系统实现的难度，尤其对于星载 平台。当前国内外普遍采用观测时间接近的多幅星载重访sar图像对高信噪比 稳定强点目标进行三维成像。差分层析sar(differential tomographic sar, d-tomosar)技术利用长时间序列sar图像在高程向与时间维合成孔径来四维 成像，其中四维信息包括三维数字地形与高程随时间变化(简称形变)。借助卫 星平台高宽幅、高时空分辨率等优势，星载差分层析sar技术可对运动目标四 维成像。因此，差分层析sar技术综合了时序干涉sar高精度测量地表形变以 及层析sar三维高分辨率成像的能力，还能分辨出叠加在一个分辨单元内不同 运动速率的多个目标。星载差分层析sar已被广泛应用于国民经济建设以及自 然科学研究等众多领域，如：地下水、矿产资源等资源调查领域，地面沉降、同 位震波偏移、山体滑坡等环境监测领域，城市规划、铁路选址与监测等工业工程 领域，火山、动土、冰川等自然科学领域，建筑、桥梁、高压电塔等人工目标三 维可视化及稳定性监测，地形测绘、三维战场模拟、武器精确打击等军事侦察领 域，森林生物量与结构调查、森林植被与农业庄稼高度测量、冰层厚度测量等。
4.当前，基于压缩感知(compressive sensing,cs)的差分层析sar成像方法， 由于其超高分辨能力及对旁瓣噪声的抑制能力，因而被广泛应用于差分层析sar 成像参数谱估计中，以德宇航zhu等人在其论文“tomographic sar inversion byl1 norm regularization—the compressive sensing approach”(《ieeetransaction on geoscience and remote sensing》，2010)提出的 sl1mmer(scale-down by l1 norm minimization,model selection,and estimationreconstruction)方法最为著名。随后，结合gis数据寻找相同稀疏支撑的堆栈， zhu等人在其论文“joint sparsity in sar tomography for urban mapping”(《ieeejournal of selected topics in signal processing》，2015)应用分布式压缩感知 (distributed compressive sensing,dcs)模型实现了更高精度的差分层析sar 成像，但该方法的预处理流程较为复杂，对地理信息系统(geographic informationsystem，gis)辅助数据的精度要求过高。此外，两种方法均未能实现对均匀场 景分布式散射体目标的高程与形变速率信息获取，因而限制了差分层析sar
的 实用化。


技术实现要素：

5.本发明针对低采样率条件下强散身体目标高程与形变速率难以高精度估计、 分布式散射体目标的高程与形变速率难以获取的问题，提供一种基于分布式压缩 感知的星载差分层析sar成像方法，无需辅助数据且能够对分布式散射体参数 高精度估计，突破复杂场景下多类目标差分层析sar 3d/4d信息重构这一关键 问题，解决了基于dcs方法对gis辅助数据的精度要求过高且处理流程复杂以 及分布式散射体难以成像等问题。
6.本发明的上述目的通过以下技术方案实现：
7.一种基于分布式压缩感知的星载差分层析sar成像方法，包括以下步骤：
8.步骤1、预处理sar图像，得到差分层析sar数据堆栈；
9.步骤2、基于ks检测理论，识别差分层析sar数据堆栈中相邻同质像元， 判定sar图像像元内是否可能存在分布式散射体和强散射体；
10.步骤3、假设模型散射体数目为1时，估计sar图像像元内能量最强散射体 的残余高程与形变速率；
11.步骤4、构建有无强散射体假设检验问题模型，逐像元判别sar像元内有 无强散射体；
12.步骤5、对于步骤4中判别为无强散射体的sar图像像元且步骤2中判别 为可能存在分布式散射体的sar图像像元，对sar图像像元内分布式散射体的 高程与分布式散射体的形变速率信息进行估计；
13.步骤6、对于步骤4中判别为有强散射体的sar图像像元且步骤2中判别为 可能存在强散射体的sar图像像元，判别sar图像像元内有单个还是双个散射体；
14.步骤7、对于步骤6中有强散射体的sar图像像元，估计sar图像像元内 强散射体的高程与形变速率；
15.步骤8、根据步骤5估计的sar图像像元内分布式散射体的高程、分布式散 射体的形变速率以及散射系数，还根据步骤7中估计的sar图像像元内强散射 体的高程、形变速率以及散射系数，逐个像元判定sar图像像元内分布式散射 体和强散射体的有效数目，并获得有效的分布式散射体和强散射体的高程与形变 速率信息，实现星载差分层析sar成像。
16.如上所述的步骤1包括以下步骤：对所有sar图像进行幅度及相位校正， 然后基于地理编码原理，对每幅sar图像进行去斜操作，得到差分层析sar数 据堆栈。
17.如上所述的步骤2包括以下步骤；
18.步骤2.1、采用两个像素统计所得的概率分布函数差的绝对值dm来度量任 意sar图像像元p与相邻sar图像像元p'的相似性；
19.步骤2.2、设定阈值c与置信度α，概率分布函数差的绝对值dm的概率分 布函数为h(t)，若dm≤c且α≥1-h(t)时，则两个sar图像像元识别为相邻同 质像元；否则，两个sar图像像元不为相邻同质像元；
20.步骤2.3、设置第一邻域窗口和第二邻域窗口，第二邻域窗口尺寸小于第一 邻域窗口尺寸；若在第一邻域窗口内与sar图像像元相邻的相似同质像元数目大 于第一设定阈值，则判定该sar图像像元内可能存在分布式散射体；若在第二邻 域窗口与sar图像像元相
邻的相似同质像元数目大于第二设定阈值，则判定该 sar图像像元内可能存在强散射体，且该强散射体像元具有多个相似同质像元。
21.如上所述的步骤4包括以下步骤；
22.构建有无强散射体假设检验问题模型：
[0023][0024]
其中与分别表示sar像元内有强散射体的假设、sar像元内无强散射 体的假设，和分别表示在有强散射体的假设和无散射体假设条 件下的概率密度函数，为差分层析sar数据堆栈中任意同名像元的平均亮度， sar时序观测值构成的矢量x＝[x1,
…
,xm]
t
，t为转置符，ε表示热噪声标准差，为sar图像像元内强散射体的高程，为sar图像像元内强散射体的形变速 率；
[0025]
若成立，则sar图像像元被判定为无强散射体；否则， sar图像像元被判定为有强散射体，k为sar图像像元内强散射体的个数，p0与 pk表示sar图像像元内没有强散射体的先验概率与有k个强散射体的先验概率。
[0026]
如上所述的步骤5包括以下步骤：对于步骤4中无强散射体的sar图像像 元且步骤2中判别可能存在分布式散射体的sar图像像元，采用多视压缩感知模 型对sar图像像元内分布式散射体的高程与分布式散射体的形变速率信息 进行估计，
[0027]
多视压缩感知模型如下式所示：
[0028][0029]
其中，||
·
||f表示frobenius范数，||||1为矢量的1范数，g表示多个差分层析 sar数据堆栈按列依次联合而成的观测数据矩阵，γ1表示待估计的多个sar图 像像元内的相同分布式散射体散射系数按列依次联合而成的矩阵，即 γ1＝[γ,γ,
…
,γ]，γ表示待估计的sar图像像元的相同分布式散射体散射系数， 而g表示所有相邻同质像元的sar时序观测值构成的矢量联合而成的矩阵，观测 矩阵φ的m行n列元素为其中δhn与vn分 别为分布式散射体的高程与形变速率的离散值。
[0030]
如上所述的步骤7包括以下步骤：
[0031]
对于步骤6中判别为有强散射体且步骤2中判别为无相邻同质像元的sar 图像像元，采用基于压缩感知的sl1mmer方法估计像元内强散射体的高程与形变 速率；
[0032]
对于步骤6中判别为有强散射体且步骤2中判别为有相邻同质像元的sar 图像像元，采用基于分布式压缩感知的方法估计具有相邻同质像元的sar图像像 元内强散射体的高程与形变速率。
[0033]
如上所述的步骤8包括以下步骤：
[0034]
采用基于贝叶斯信息准则的惩罚约束判别准则逐个sar图像像元判定分布 式散射体的有效个数和强散射体的有效个数，为有效个数，基于贝叶斯信息准 则的惩罚约束判别准则基于以下公式：
[0035][0036]
其中，n为高程与形变速率向离散化间隔的数目，k表示sar图像像元内 可能的散射体的数目，γk表示步骤5和7估计的sar像元内散射体(包括分布 式散射体和强散射体)的散射系数γ能量由强到弱排列后取其前k个数值构成的 向量，φk表示观测矩阵φ中与γk在散射系数γ相同位置的k列向量构成的矩 阵，
[0037]
记sar图像像元内分布式散射体的有效个数为记录能量前名的有 效分布式散射体的位置，对应截取步骤5估计的sar像元内对应位置的分布式 散射体的高程、形变速率以及散射系数，即获得有效的分布式散射体的高程与形 变速率信息，
[0038]
记sar图像像元内强散射体的有效数目记录能量前名的有效强散 射体的位置，对应截取步骤7估计的sar图像像元内对应位置的强散射体的高 程、形变速率以及散射系数，即获得有效的强散射体的高程与形变速率信息。
[0039]
本发明相对于现有技术，具有以下有益效果：
[0040]
(1)由于联合分布式散射体相邻同类像元进行多视压缩感知重构，本发 明提出的一种基于分布式压缩感知的星载差分层析sar成像方法，能够有效降 低旁瓣噪声被错误成像的风险且提高分布式散射体目标参数的估计精度，解决了 由分布式散射体像元导致的成像产品信息空白问题。
[0041]
(2)由于联合强散射体相邻同类像元进行联合压缩感知重构，本发明提 出的一种基于分布式压缩感知的星载差分层析sar成像方法，能够有效提高相 邻同类强散射体的成像处理精度以及高维分辨能力。
附图说明
[0042]
图1为本发明的流程图
[0043]
图2为基于分布式压缩感知的星载差分层析sar成像结果地理编码后附在 google earth上的高程与形变速率图，其中(a)(b)分别为高程与形变速率图；
[0044]
图3为德宇航sl1mmer成像方法与本发明基于分布式压缩感知成像方法 的结果对比，其中(a)(c)(e)分别为sl1mmer方法获得的高度图、平均形变速率 图与散射系数图，而(b)(d)(f)分别为本发明方法获得的高度图、平均形变速率图 与散射系数图。
具体实施方式
[0045]
为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明 作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发 明，并不用于限定本发明。
[0046]
一种基于分布式压缩感知的星载差分层析sar成像方法，如图1所示，包 括以下步骤：
[0047]
步骤1：预处理sar图像，得到差分层析sar数据堆栈；
[0048]
步骤1具体包括以下步骤：
[0049]
对所有sar图像进行幅度及相位校正，然后基于地理编码原理，对每幅sar 图像进行去斜操作，得到差分层析sar数据堆栈。
[0050]
步骤2：基于kolmogorov-smirnov(ks)检测理论，识别差分层析sar数据 堆栈中相邻同质像元；
[0051]
步骤2具体包括以下步骤：
[0052]
步骤2.1：采用两个像素统计所得的概率分布函数差的绝对值dm来度量任 意sar图像像元p与相邻sar图像像元p'的相似性，即
[0053][0054]
其中m为sar图像的总数，r为实数，与分别表示随机变量x 在任意sar图像像元p及其相邻sar图像像元p'位置的累积分布函数，而
[0055][0056][0057]
其中表示差分层析sar数据堆栈在同一sar图像像元p处的观测样 本值，表示差分层析sar数据堆栈在同一相邻sar图像像元p'处的观 测样本值，m为sar图像序号，m∈{1,2,...,m}，|
·
|表示复数据的幅度算子；
[0058]
步骤2.2：设定阈值c与置信度α，概率分布函数差的绝对值dm的分布逼 近于ks分布，概率分布函数差的绝对值dm的概率分布函数h(t)为
[0059][0060]
其中，p(
·
)表示条件概率密度函数，t为相似度阈值常数，e为自然对数的 底数，z为中间变量，h(t)不依赖于单个像素数据概率分布函数。若dm≤c且 α≥1-h(t)时，则两个sar图像像元具有相同的统计分布，则两个sar图像像元 可识别为相邻同质像元；否则，两个sar图像像元不具有相似性，即不为相邻同 质像元。
[0061]
步骤2.3：设置第一邻域窗口为9*9像元大小，若在第一邻域窗口内与sar 图像像元相邻的相似同质像元数目大于20(第一设定阈值)，则判定该sar图像 像元内可能存在分布式散射体；设置第二领域窗口大小为3*3像元大小，若在第 二邻域窗口与sar图像像元相邻的相似同质像元数目大于3(第二设定阈值)， 则判定该sar图像像元内可能存在强散射体，且该强散射体像元具有多个相似同 质像元。
[0062]
步骤3：假设模型散射体数目为1时，估计sar图像像元内能量最强散射体 的残余高程与形变速率；
[0063]
步骤3具体包括以下步骤：
[0064]
采用基于振幅归一化度量的周期图最大化方法，估计散射体数目为1的假设 下的任意sar图像像元内能量最强散射体的残余高程δh与sar图像像元内能量 最强散射体的形变速率v，即
[0065][0066]
其中，式(4)的b
⊥m与tm分别为第m次观测中相对主图像的有效基线与时间 差，xm表示时序sar第m次观测值，*表示共轭算子，j表示虚数单位，λ为雷 达波长，r为天线相位中心到参考地形的斜距，θ为入射角。
[0067]
步骤4：构建有无强散射体假设检验问题模型，逐像元判别sar像元内有 无强散射体；
[0068]
上述步骤4包括以下步骤：
[0069]
构建有无强散射体假设检验问题模型，采用广义似然比检测(generalized likelihood ratio test,glrt)理论贝叶斯决策判别sar图像像元内有无强散射 体，构建的有无强散射体假设检验问题模型，即：
[0070][0071]
其中与分别表示sar像元内有强散射体的假设、sar像元内无强散射 体的假设，和分别表示在有强散射体的假设和无散射体假设条 件下的概率密度函数，为差分层析sar数据堆栈中任意同名像元的平均亮度， sar时序观测值构成的矢量x＝[x1,
…
,xm]
t
，t为转置符，ε表示热噪声标准差，为sar图像像元内强散射体的高程，为sar图像像元内强散射体的形变速 率，采用glrt理论的贝叶斯决策为：
[0072][0073]
其中，k为sar图像像元内强散射体的个数，p0与pk表示sar图像像元内没 有强散射体的先验概率与有k个强散射体的先验概率，均设置为50％。若(6)式 成立，则sar图像像元被判定为无强散射体；反之，sar图像像元被判定为有强 散射体。
[0074]
步骤5：对于步骤4中判别为无强散射体的sar图像像元且步骤2中判别 为可能存在分布式散射体的sar图像像元，对sar图像像元内分布式散射体的 高程与分布式散射体的形变速率信息进行估计；
[0075]
上述步骤5包括以下具体步骤：
[0076]
对于步骤4中无强散射体的sar图像像元且步骤2中判别可能存在分布式 散射体
的sar图像像元，采用分布式压缩感知理论中的多视压缩感知(multilook compressive sensing,mcs)模型谱估计方法对sar图像像元内分布式散射体的 高程与分布式散射体的形变速率信息进行估计。所述多视压缩感知 (multilook compressive sensing,mcs)模型如下
[0077][0078]
其中||
·
||f表示frobenius范数，||||1为矢量的1范数，g表示多个差分层析 sar数据堆栈按列依次联合而成的观测数据矩阵，γ1表示待估计的多个sar图 像像元内的相同分布式散射体散射系数按列依次联合而成的矩阵，即 γ1＝[γ,γ,
…
,γ]，γ表示待估计的sar图像像元的相同分布式散射体散射系数， 而g表示所有相邻同质像元的sar时序观测值构成的矢量联合而成的矩阵，观测 矩阵φ的m行n列元素为其中δhn与vn分 别为分布式散射体的高程与形变速率的离散值。
[0079]
步骤6：对于步骤4中判别为有强散射体的sar图像像元且步骤2中判别为 可能存在强散射体的sar图像像元，判别sar图像像元内有单个还是双个散射体 情形；
[0080]
上述的步骤6包括以下步骤：
[0081]
若有强散射体的sar图像像元的幅度离差小于或等于所有sar图像平均幅度 的70％，则sar图像像元为单个散射体像元；若有散射体的sar图像像元的幅度 离差大于所有sar图像平均幅度的70％，则该sar图像像元内为有双个散射体叠 加。
[0082]
步骤7：对于步骤6中有强散射体的sar图像像元，估计sar图像像元内 强散射体的高程与形变速率
[0083]
步骤7具体包括以下步骤：
[0084]
对于步骤6中判别为有强散射体且步骤2中判别为无相邻同质像元的sar 图像像元，采用式(8)基于压缩感知(compressive sensing,cs)的sl1mmer 方法估计像元内强散射体的高程与形变速率
[0085][0086]
对于步骤6中判别为有强散射体且步骤2中判别为有相邻同质像元的sar 图像像元，采用式(9)基于分布式压缩感知(distributed compressive sensing, dcs)的方法估计具有相邻同质像元的sar图像像元内强散射体的高程与形 变速率
[0087][0088]
其中，||
·
||
2,1
表示矩阵行的l2范数和的最小化，γ2表示待估计的多个sar图像 像元内的不同强散射体散射系数按列依次联合而成的矩阵，即γ2＝[γ1,γ2,
…
,γ
l
]， [γ1,γ2,
…
,γ
l
]表示待估计的sar图像像元的不同强散射体散射系数。
[0089]
步骤8：根据步骤5估计的sar图像像元内分布式散射体的高程分布 式散射体的形变速率以及散射系数，还根据步骤7中估计的sar图像像元内 强散射体的高程形变速率以及散射系数，逐个像元判定sar图像像元 内分布式散射体和强散射体的有效数目，并获得有效的分布式散射体和强散射体 的高程与形变速率信息，实现星载差分层析
sar成像。
[0090]
具体包括以下步骤：
[0091]
采用基于贝叶斯信息准则(bayesian information criterion,bic)的惩罚 约束判别准则逐个sar图像像元判定分布式散射体的有效个数和强散射体的有 效个数，为有效个数，基于贝叶斯信息准则(bayesian information criterion, bic)的惩罚约束判别准则基于公式(10)：
[0092][0093]
其中n为高程与形变速率向离散化间隔的数目，k表示sar图像像元内可能的 散射体的数目，其中γk表示步骤5和7估计的sar像元内散射体(包括分布式 散射体和强散射体)的散射系数γ能量由强到弱排列后取其前k个数值构成的向 量，φk表示观测矩阵φ中与γk在散射系数γ相同位置的k列向量构成的矩阵。
[0094]
最后，记sar图像像元内分布式散射体的有效个数为记录能量前名 的有效分布式散射体的位置，对应截取步骤5估计的sar像元内对应位置的分 布式散射体的高程、形变速率以及散射系数，即获得有效的分布式散射体的高程 与形变速率信息。
[0095]
记sar图像像元内强散射体的有效数目记录能量前名的有效强散 射体的位置，对应截取步骤7估计的sar图像像元内对应位置的强散射体的高 程、形变速率以及散射系数，即获得有效的强散射体的高程与形变速率信息。
[0096]
对比例：
[0097]
通过terrasar-x/tandem-x获取的31幅首都国际机场t3 e航站楼区域的 真实数据，验证本发明一种基于分布式压缩感知的星载差分层析sar成像方法 的有效性。31幅terrasar-x/tandem-x图像获取时间与基线参数信息如表1所 示，表2列出了与本发明有关的terrasar-x/tandem-x差分层析sar成像系统 参数。将基于分布式压缩感知的差分层析sar处理结果与当前先进的同样不依 赖于外部辅助数据的sl1mmer方法的结果对比，说明本发明所提方法的优势。 处理中对应的垂直于水平面的高度范围设定为70m，采样点数为601；沿视线向 形变速率范围设定为20mm/year，采样点数为101。
[0098]
表1 31幅terrasar-x/tandem-x图像获取时间与基线参数
[0099]
[0100][0101]
表2terrasar-x/tandem-x差分层析sar成像系统参数
[0102]
[0103][0104]
将本发明基于分布式压缩感知方法估计的高程与形变速率信息经过地理编 码后附在google earth上，结果如图2所示，可以看出本发明基于分布式压缩 感知方法估计的t3 e三维空间位置与google光学图像匹配得较好，这表明了本 发明基于分布式压缩感知方法估计的三维空间位置信息的准确性，由于形变速率 信息与高程信息是同时估计的，进一步，表明了本发明基于分布式压缩感知方法 估计形变速率信息的准确性。从图3可以看出，德宇航sl1mmer成像方法与本发 明基于分布式压缩感知方法获取的高程([m])与形变速率([mm/year])信息以及 散射系数整体上相一致，表明了本发明基于分布式压缩感知方法的有效性。
[0105]
随后，统计了德宇航sl1mmer成像方法与本发明基于分布式压缩感知方法重 构的像元数目，表3列出了sl1mmer与本发明基于分布式压缩感知方法重构的像 元数目以及本发明基于分布式压缩感知方法三种模型各自重构的像元数目。由于 绝大多数分辨单元内的强散射体个数不超过2，所以将像元内散射体个数的最大 阈值设置为2。本发明基于分布式压缩感知方法比sl1mmer方法能够重构出多余 的134个分布式散射体像元，且本发明基于分布式压缩感知方法重构出的叠掩散 射体更多。
[0106]
表3德宇航sl1mmer方法与本发明方法重构的像元数
[0107]
差分层析sar方法sl1mmer本文方法csdcsmcs重构的全部像元数266542678826449205134一个散射体像元数251122518424941109134两个散射体像元数154216041508960
[0108]
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施 例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于 本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润 饰，应视为本发明的保护范围。