CN116400356B - 一种基于同质区域联合的层析sar三维成像方法 - Google Patents
一种基于同质区域联合的层析sar三维成像方法 Download PDFInfo
- Publication number
- CN116400356B CN116400356B CN202310665210.1A CN202310665210A CN116400356B CN 116400356 B CN116400356 B CN 116400356B CN 202310665210 A CN202310665210 A CN 202310665210A CN 116400356 B CN116400356 B CN 116400356B
- Authority
- CN
- China
- Prior art keywords
- homogeneous region
- pixels
- pixel
- sar
- dimensional imaging
- 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
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 59
- 238000000034 method Methods 0.000 claims abstract description 48
- 238000011084 recovery Methods 0.000 claims abstract description 13
- 230000009466 transformation Effects 0.000 claims abstract description 4
- 239000013598 vector Substances 0.000 claims description 35
- 239000011159 matrix material Substances 0.000 claims description 23
- 230000006870 function Effects 0.000 claims description 13
- 230000011218 segmentation Effects 0.000 claims description 11
- 230000001427 coherent effect Effects 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000004590 computer program Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 5
- 238000005315 distribution function Methods 0.000 claims 1
- 230000000007 visual effect Effects 0.000 claims 1
- 238000012545 processing Methods 0.000 abstract description 14
- LWYJUZBXGAFFLP-OCNCTQISSA-N menogaril Chemical compound O1[C@@]2(C)[C@H](O)[C@@H](N(C)C)[C@H](O)[C@@H]1OC1=C3C(=O)C(C=C4C[C@@](C)(O)C[C@H](C4=C4O)OC)=C4C(=O)C3=C(O)C=C12 LWYJUZBXGAFFLP-OCNCTQISSA-N 0.000 abstract description 11
- 238000009825 accumulation Methods 0.000 abstract description 4
- 238000013461 design Methods 0.000 abstract description 3
- 238000004422 calculation algorithm Methods 0.000 description 30
- 238000004088 simulation Methods 0.000 description 9
- 238000005457 optimization Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000012163 sequencing technique Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 3
- 238000011524 similarity measure Methods 0.000 description 3
- 239000000654 additive Substances 0.000 description 2
- 230000000996 additive effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000003708 edge detection Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000017525 heat dissipation Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
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
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)
- Image Processing (AREA)
Abstract
本发明公开了一种基于同质区域联合的层析SAR三维成像方法，属于雷达成像处理技术领域。所述方法设计了一种结合欧氏距离、像素强度和边缘信息的测度函数，将SAR图像中三维位置相近、电磁散射强度相似的像素点分割成局部连通的同质区域。然后，将TomoSAR单像素观测信号模型扩展到同质区域联合观测信号模型，通过同质区域积累得到高程约束中心，根据同质区域大小和成像参数自动确定高程搜索范围，稀疏逼近估计同质区域高程。最后，对斜高进行几何变换和地理编码，得到三维成像结果。该方法能减少成像结果中由固有相干斑噪声带来的噪点，有效提升了信噪比，提高了对弱小目标的恢复能力。
Description
技术领域
本发明涉及雷达成像处理技术领域，具体涉及一种基于同质区域联合的层析SAR三维成像方法。
背景技术
层析合成孔径雷达（Tomographic Synthetic Aperture Radar, TomoSAR）技术可以获得城市建筑、陡峭地形等复杂对象的三维结构和电磁散射信息，显著提升场景解译和目标识别能力，已经成为SAR技术发展的重要前沿方向。
TomoSAR 3维成像在距离向和方位向与传统2维SAR成像相同，距离向通过发射大带宽信号结合脉冲压缩技术获得距离分辨率，方位向通过合成孔径获取分辨能力。为了获得高程向分辨能力，利用航迹精确控制、满足奈奎斯特定律的多次观测构建高程向的等效阵列，实现高度维的合成孔径，获得高程向分辨能力。
图1展示了TomoSAR的成像几何，每个灰度扇形区域表示波束覆盖范围，对应生成一幅SAR图像，不同航迹获得的较小视角差异的多幅SAR图像，经图像配准后，逐像素解算目标的高度信息。
假设对某个区域有个相干的观测，相当于在垂直斜距的高程向（也即斜高向）上
形成了等效合成孔径，单个像素的连续信号模型可以表示为：
，
其中，表示第个观测采集到的复数观测值，是斜高向散射体的散射系数
分布，是由等效基线、雷达波长和斜距决定的等效空间频率，为噪
声。
对连续信号模型进行离散化后，可以表达为矩阵形式：
其中,是大小的观测矩阵，为高程向上划分的网格数量，矩阵中第行第列的元素为，为第次观测的等效空间频率，为第个预先划分的
网格位置，是噪声向量。将离散化的信号模型展开得到如下形式：
该方程中的未知数为，每个的值代表了不同位置网格的信号
强度，对进行估计可以得到目标的三维位置和散射强度。
目前，压缩感知是城市场景TomoSAR三维重建最常用的方法，它较好的解决了稀疏观测条件下超分辨率重构问题，主要包括凸优化类算法、贪婪类算法和贝叶斯类算法等。凸优化类算法的代表是Xiaoxiang Zhu等提出了基于L1范数正则化的压缩感知TomoSAR方法（Scale-down by L1 norm Minimization, Model selection, and EstimationReconstruction，SL1MMER），该方法通过L1范数的凸优化重构稀疏信号，同时采用一个正则化因子调整噪声水平。研究验证了压缩感知方法具备高程向的超分辨率，在多个散射体的识别中性能良好。Yilei Shi等将Nonlocal滤波方法与SL1MMER相结合，用于建筑物的精确三维成像，在提高高度向分辨率和重构精度的同时，保持了重构目标距离-方位向的分辨率。凸优化方法的不足是计算量较大，较难适用与大范围三维成像和应急测绘。贪婪类算法的代表是匹配追踪算法（Matching Pursuit, MP），通过在观测信号与观测残差之间进行匹配来恢复信号，包括正交匹配追踪算法（Orthogonal MP, OMP）、正则化正交匹配追踪算法（Regularized OMP，ROMP）、广义正交匹配追踪算法（Generalized OMP, GOMP）等。此类算法不直接求解信号的稀疏解，而是通过稀疏逼近的方式渐进的逼近测量值，极大的减少了运算量，但此类方法需要将信号的稀疏度作为先验信息，在复杂区域信号的稀疏度（即叠掩数）往往是未知的，迭代停止的条件只能通过预先估计叠掩数或者设置阈值的方式来确定。贝叶斯类算法的代表是稀疏贝叶斯算法，通过对待恢复的信号进行先验概率假设，结合观测数据的分布概率以及贝叶斯准则，对信号后验概率的最大值进行学习，以此来获得信号功率的空间分布。Zekun Jiao等提出了一种离格稀疏贝叶斯学习方法，用于阵列SAR系统的三维成像，通过局部高斯马尔可夫随机场(Local Markov Random Field ，LGMRF) 来探索和利用SAR图像中矩形邻域像素的之间的结构联系，在一定程度上克服了传统逐像素三维成像数据处理过程中噪声或其他因素导致的异常值和伪影影响。
综上所述，现有压缩感知算法大多逐像素对进行估计。贪婪类算法中的OMP算法
根据目标可由有限个观测矩阵中的原子线性表示的思想，通过所选择的观测矩阵中的原子
与观测信号的内积比较，每次从观测矩阵中选出与观测信号最相关的原子，将其
纳入重构稀疏基。计算该稀疏基下得到的估计值与之间的残差，此时将之前已经选出的
原子从观测矩阵中剔除，以保证每次选出的原子之间是正交的，并且每次迭代的残差与稀
疏基中每个原子也是正交的。再在观测矩阵中寻找与迭代残差最匹配的原子。由于正交化
处理，使得算法能够快速收敛，得到满足迭代误差精度的重构稀疏基，从而重构出原始层析
向的散射向量。
逐像素进行估计的方式无法利用邻近像素的信息，易受到SAR图像中固有的相干斑噪声的影响，产生大量离群值，且当观测矩阵中原子相关性较高时，算法的重构精度会变差。
TomoSAR三维重建性能主要由多角度观测次数与成像处理方法决定。通过增加观测次数来实现重建性能提升的代价成本极高。因此，亟需创新发展高性能三维成像处理方法，降低TomoSAR技术应用成本。
发明内容
本发明要解决的技术问题是，提供一种能够利用邻域像素的信息，减少固有相干斑噪声影响的层析SAR三维成像方法、装置和存储介质。
为解决上述技术问题，本发明提供一种基于同质区域联合的层析SAR 三维成像方法，包括以下步骤：
S1，对SAR图像进行同质区域分割；
S2，进行同质区域的联合稀疏恢复，计算同质区域所有像素斜高估计值；
S3，进行几何坐标变换，将斜高向转换为实际空间的高度向，得到三维成像结果。
如上所述步骤S1中对SAR图像进行同质区域分割的步骤为：
S101，在SAR图像上均匀初始化尺寸为像素的个网格聚类中心，将聚类中
心以元素的形式存入优先级队列，其中、为大于1的整数，
为空间位置，为切片的强度向量，为边缘强度信息，元素的第四列为，表
示第个待处理像素点到第个聚类中心的相似度，初始化聚类中心的为0，初始化与
SAR图像尺寸相同的标签图；
S102，计算聚类中心四邻域像素到所述聚类中心自身的相似度，并以元素的形式
推入优先级队列，每一次向优先级队列中推入一个元素，更新一次排序，使优先级队列中元素保持从大到小的顺序，推出优先级队列的顶端元素，所述顶端元素为与所述
聚类中心相似度最大的像素，将所述聚类中心的标签赋给顶端元素对应的像素，用顶端
元素的空间位置和切片强度向量更新所述聚类中心的空间位置和切片强度；
S103，依次计算顶端元素的四邻域像素点到聚类中心的相似度，创建相应的元
素，若所述顶端元素的四邻域像素点已被标记，则已被标记的像素点所创建的元素不推
入优先级队列，只将未标记的像素点创建的元素推入优先级队列，每一次向优先级队
列中推入一个元素，更新一次排序，使优先级队列中元素保持从大到小的顺序，推出优
先级队列的顶端元素，将所述聚类中心的标签赋给顶端元素对应的像素，用顶端元素的空间位置和切片强度向量更新所述聚类中心的空间位置和切片强度；
S104，若所有的像素点均被标记，并且优先级队列为空，输出同质区域聚类标签
图，否则返回步骤S103。
如上所述步骤S101所述相似度按照以下公式计算：
，
其中，为第个待处理像素点到第个聚类中心的相似度，为为第个
待处理像素点到第个聚类中心的空间相似度，为为第个待处理像素点到第个聚
类中心的像素强度相似度，为第个待处理像素点的边缘强度，、为加权系数；
所述空间相似度按照以下公式计算：，
其中，欧氏距离，、为第像素点的坐标，、
为第个聚类中心的坐标；
所述像素强度相似度按照以下公式计算：，
其中，是像素强度比值，，为待处理像素切片的强度，为聚类
中心切片的强度，是的概率密度函数，表示切片所包含的像素点个数，
为标准高斯核函数；所述概率密度函数按照以下公式计算：，其中，为gamma分布函数，为成像视数。
优选的，上述相似度公式中，加权系数的取值范围为[2, 6]，加权系数的取值
范围为[0.3, 0.7]。
优选的，步骤S101中所述初始化网格尺寸的取值为5或7。
优选的，步骤S2中进行同质区域的联合稀疏恢复，计算同质区域所有像素斜高估计值的步骤为：
S201，初始化斜高估计值，，稀疏集，残差，其中，
为观测信号向量，，为同质区域内部的第个像素的观测信号，，为同质区域内部的像素个数，，N为相干观测
的个数；
S202，构建同质区域内的观测矩阵，计算观测矩阵与观测信号向
量的相关性向量，其中，为同质区域内的第个像素的观测矩阵，，为第个相干观测的等效空间频率，，为第个
预先划分的高程位置，为高程向上划分的网格数量；
S203，根据相关性向量主峰数确定同质区域的叠掩数；设置k=1；
S204，计算观测矩阵与残差的相关性向量；
S205，由相关性向量的第一主峰确定约束中心；
S206，在约束中心的范围内搜索每个像素的极大值；
S207，将极大值对应的观测矩阵中的原子纳入同质区域重构稀疏基，
，其中，；
S208，计算该稀疏基下斜高估计值，，其中， ，为共轭转置；
S209，计算该稀疏基下的残差，，其中，；
S210，；
S211，若则结束并输出，否则返回步骤S204。
优选的，步骤S202中所述计算观测矩阵与观测信号向量的相关性向量的
步骤为：
S2021，计算同质区域内部每个高程位置与观测信号的相关性，，其中，为同质区域内部第个像素在高程位置
与观测信号的相关性；
S2022，将不同像素同一高度位置的相关性进行非相干叠加，拟合整个同质区域的
相关性向量， ；
优选的，步骤S204中所述计算观测矩阵与残差的相关性向量的步骤为：
S2041，计算同质区域内部每个高程位置与残差的相关性，，其中，为同质区域内部第个像素在高程位置
与残差的相关性；
S2042，将不同像素同一高度位置的相关性进行非相干叠加，拟合整个同质区域的
相关性向量， 。
优选的，步骤S206中所述约束中心的范围计算方法为：
，其中，为同质区域网格长度，为下视角，为距离向像素
尺寸。
本发明还提供一种基于同质区域联合的层析SAR 三维成像装置，包括非易失性存储器和一个或多个处理器，所述非易失性存储器中存储有可执行代码，其特征在于，所述处理器执行所述可执行代码时，用于实现前面所述的基于同质区域联合的层析SAR 三维成像方法。
本发明还提供一种计算机可读存储介质，所述计算机可读存储介质包括存储的计算机程序，其中，在所述计算机程序被处理器运行时控制所述存储介质所在设备执行前面所述基于同质区域联合的层析SAR 三维成像方法。
本发明相对于现有技术，具有以下有益效果：
针对现有TomoSAR 3维成像算法逐像素求解高程信息，无法利用邻域像素蕴含的丰富信息的问题，设计了一种结合欧式距离，像素强度和边缘信息的相似度测度，将三维位置相近、电磁散射强度相似的像素点聚集成一个局部同质连通区域，通过对同质区域联合进行稀疏恢复，设置约束中心对对高程范围进行有效约束。对于同质区域内部，高程上具有一定的相关性：在平地区域，同质区域内部的高程值完全一致；在斜坡区域，高程呈现上升趋势，但同质区域内部应当具有连续的特性，高程不会发生较大突变。充分利用邻域像素蕴含的结构、纹理等信息，减少了成像结果中由固有相干斑噪声带来的噪点，有效提升了信噪比，进而提高了对弱小目标的恢复能力，将微波和视觉做了一个较好的结合，有效提升了三维成像的效果。
附图说明
图1是层析SAR成像几何示意图。
图2是以同质区域为单位在高程向划分网格的示意图。
图3是仿真实验中，对双散射体逐像素估计的相关性和非相干累加的相关性的对比图。
图4是本发明实施例的流程图。
图5是本发明实施例中的仿真数据集。
图6是本发明实施例中，仿真实验同质区域分割结果图。
图7是采用4种算法分别对仿真数据集进行三维成像的结果图。
图8是本发明实施例提供的基于同质区域联合的层析SAR 三维成像装置结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”和 “包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当进一步理解，在本发明说明书和所附权利要求书中使用的术语“和/ 或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
具体实施例1
一种基于同质区域联合的层析SAR 三维成像方法，如图4所示，包括以下步骤：
S401，初始化C个聚类中心，具体步骤为：
在SAR图像上均匀初始化尺寸为像素的个网格聚类中心，将聚类中心以元
素的形式存入优先级队列，其中，为大于1的整数，为空间
位置，为切片的强度向量，为边缘强度信息，元素的第四列为，表示第
个待处理像素点到第个聚类中心的相似度，初始化聚类中心的为0，初始化与SAR图
像尺寸相同的标签图。
S402，根据位置和散射强度进行同质区域聚类，具体步骤为：
S4021，计算聚类中心四邻域像素到所述聚类中心自身的相似度，并以元素的形
式推入优先级队列，每一次向优先级队列中推入一个元素，更新一次排序，使优先级队
列中元素保持从大到小的顺序，推出优先级队列的顶端元素，所述顶端元素为与所
述聚类中心相似度最大的像素，将所述聚类中心的标签赋给顶端元素对应的像素，用顶
端元素的空间位置和切片强度向量更新所述聚类中心的空间位置和切片强度；
S4022，依次计算顶端元素的四邻域像素点到聚类中心的相似度，创建相应的元
素，若所述顶端元素的四邻域像素点已被标记，则已被标记的像素点所创建的元素不推
入优先级队列，只将未标记的像素点创建的元素推入优先级队列，每一次向优先级队
列中推入一个元素，更新一次排序，使优先级队列中元素保持从大到小的顺序，推出优
先级队列的顶端元素，将所述聚类中心的标签赋给顶端元素对应的像素，用顶端元素的空间位置和切片强度向量更新所述聚类中心的空间位置和切片强度；
S4023，若所有的像素点均被标记，并且优先级队列为空，则执行步骤S403，否则
返回步骤S4022。
S403，输出同质区域分割标签图。
S404，确定同质区域的K个约束中心，具体步骤为：
如图2所示，初始化斜高估计值，，稀疏集，残差，其
中，为观测信号向量，，为同质区域内部的第个像素的观测
信号，，为同质区域内部的像素个数，，N为相干
观测的个数；
构建同质区域内的观测矩阵，计算观测矩阵与观测信号向量的
相关性向量，其中，为同质区域内的第个像素的观测矩阵，，为第个相干观测的等效空间频率，，为第个
预先划分的高程位置，为高程向上划分的网格数量；
根据相关性向量的主峰确定同质区域的K个约束中心。
S405，若，在第个约束中心的范围内搜索每个像素的极大值，，其中，为同质区域网格长度，为下视角，为距离向像素尺寸。
S406，将极大值对应的观测矩阵中的原子纳入同质区域重构稀疏基，，其中，，，为同质区域内部的
像素个数，N为相干观测的个数，为第次观测的等效空间频率。
S407，，其中， ，
为观测信号向量，为共轭转置；
计算该稀疏基下的残差，，其中，；
。
S408，重复步骤S405至S407，直到达到叠掩数量K，得到第c个同质区域的高程估计
值，存储第c个同质区域的点云信息（包括距离、方位、斜高程信息）。
S409，若完成所有C个同质区域的联合稀疏恢复，将所有像素点的斜高变换为实际空间的高度向，得到三维成像结果。
具体实施例2
在实施例1的基础上，增加对相似度计算的优化，采用以下结合欧氏距离、像素强度相似度以及像素边缘强度相结合的相似度测度来计算像素点到聚类中心的相似度。相似度的计算公式为：
，
其中，为第个待处理像素点到第个聚类中心的相似度，为为第个
待处理像素点到第个聚类中心的空间相似度，为为第个待处理像素点到第个聚
类中心的像素强度相似度，为第个待处理像素点的边缘强度，、为加权系数；
所述空间相似度按照以下公式计算：，
其中，欧氏距离，、为第像素点的坐标，、
为第个聚类中心的坐标；
所述像素强度相似度按照以下公式计算：，
其中，是像素强度比值，，为待处理像素切片的强度，为聚类
中心切片的强度，是的概率密度函数，表示切片所包含的像素点个数，
为标准高斯核函数；所述概率密度函数按照以下公式计算：，其中，为gamma分布函数，为成像视数。
采用上述相似度计算方法，能够有效减少SAR图像中加性和乘性噪声的干扰，并能有效提高同质区域边缘与图像真实边缘的贴合率，分析如下：
在SAR图像上，包含同一目标像素点通常位置相近，采用欧氏距离来衡量两个切片
的空间距离，再使用标准高斯核函数将该距离映射成空间相似度，如下所示：
欧氏距离对于加性噪声具有较好的鲁棒性，但对于SAR图像中的乘性噪声则不具有。为了解决这个问题，Feng等人2011年提出一种基于像素强度比值距离的方法来衡量两个受噪声影响的切片之间的相似度，对SAR图像具有较强的鲁棒性。定义像素强度比值为：
其中，为待处理像素切片的强度向量，为聚类中心切片的强度向量，表示切
片所包含的像素点个数(通常取或)，表示标准高斯核函数，表示
和的比值距离，该比值距离对SAR图像的乘性噪声有较强的鲁棒性。定义像素强度相似度
为：
其中，是像素强度比值的商。的概率密度函数(Probability Density
Function，PDF)定义为：
其中，为gamma分布函数，为成像视数。
为了区分来自不同目标的像素，提取出图像中的真实边缘作为像素点间相似度测
度的一部分，可以有效提高同质区域边缘与图像真实边缘的贴合率。采用边缘检测算法检
测SAR图像中的边缘，得到边缘强度图。边缘强度图的数值越接近1表示边缘越强，数值越接
近0表示边缘越弱。定义边缘检测结果二值图为，当待处理像素点为边缘时，降低待处理
像素点到聚类中心的相似度。
具体实施例3
在实施例2的基础上，进一步优化相似度计算公式中，加权系数的取值范围为
[2, 6]，加权系数的取值范围为[0.3, 0.7]。
加权系数、采用以上取值范围的权重值，使得空间距离，像素强度以及边缘信
息的占比合理，能够尽可能的将三维位置接近且特征相似的像素点进行聚类，生成的同质
区域边缘与地物边缘贴合，聚类效果更好。
具体实施例4
在实施例1的基础上，将初始化网格尺寸的取值设定为5或7。
当同质区域过小时，内部包含的散射强度较大的噪声，容易对结果造成影响；当同
质区域的过大时，不准确的估计将影响较多像素的恢复。因此，初始化网格尺寸设定为5
或7是最佳方案。具体试验结果数据参见仿真试验说明。
具体实施例5
在实施例1的基础上，进一步优化相关性向量的计算方法，具体步骤为：
计算同质区域内部每个高程位置与观测信号的相关性，，其中，为同质区域内部第个像素在高程位置
与观测信号的相关性；
将不同像素同一高度位置的相关性进行非相干叠加，拟合整个同质区域的相关性
向量， 。
采用这种对同一高度位置不同像素进行非相干叠加的方式，能够有效消除逐像素估计恢复存在大量受相干斑噪声影响产生的错误估计。
如图3所示，在仿真实验中，设置同质区域包含的双散射体的高度分别为0米和22.5米，这个间隔为1.5倍瑞利分辨率。从图3可以直观的看出，逐像素估计的相关性存在大量受相干斑噪声影响产生的错误估计峰值，而采用非相干累加得到的相关性向量的峰值表征了目标所在高度为0.64米和22.04米，与设置的散射体位置基本一致，且由噪声影响产生的虚假目标被很好的剔除。
为了评估发明实施例和其它TomoSAR 3维成像方法的性能，并与其它算法进行对比分析，设计了仿真SAR图像如图5所示。
该仿真SAR图像基于典型单通道SAR图像乘积模型将SAR图像测量值分解为纹理
(地物的后向散射)和相干斑分量的乘积。其中，纹理分量服从逆Gamma分布，相干斑分量服
从单位均值Gamma分布，且设置相干斑分量的视数为4。在这些参数设置条件下，仿真得到带
有乘性相干斑噪声的SAR杂波图像，图像尺寸为像素。根据公式，计算得到仿真数据，即单个像素的一次观测获
得的复数值，设置基线长度为2m，波长为0.03m，最近斜距为2000m，
叠掩数，模拟高程差变化范围为[10m,50m]，再加入SNR=5dB的高斯噪声。场景的瑞利
分辨率为15米，不模糊高度为135米。二维SAR图像如图5中（a）所示，三维模型如图5中（b）所
示，灰度表示高程信息，模拟了平地，平缓斜坡和陡峭斜坡三种情况。
如图6所示，采用本发明实施例中的方法对仿真SAR图像进行同质区域分割，设置
初始化的网格尺寸，,和，分析不同初始化网格尺寸对重建结果的影
响。图6中（a）为同质区域分割结果，图6中（b）为同质区域分割结果，图6中（c）为同质区域分割结果，图6中（d）为同质区域分割结果。可以看出同质区域分割的
结果具有贴合地物边缘，内部像素同质的特点。
表1给出了在模拟了平地，平缓斜坡和陡峭斜坡三种情况下的仿真数据集，采用
MUSIC算法、OMP算法、GOMP算法和本发明实施例中的方法进行三维成像所得图像的重构成
功率RP、曼哈顿距离MD、准确性AC、完备性CO指标的评估结果。从评价指标的结果可以看出，
同质区域的大小对成像结果有一定影响，在时，各项指标效果相对最佳。当同质区域
过小时，内部包含的散射强度较大的噪声，容易对结果造成影响；当同质区域的过大时，不
准确的估计将影响较多像素的恢复。因此建议设置或。
从表1可以看出，本发明实施例中的方法具有较高的重构成功率，较小的曼哈顿距离，表明在叠掩数量的估计和散射体高程位置估计上较其他算法更为准确；OMP算法和GOMP算法的性能几乎一致，而MUSIC算法的效果较差。
表1 采用4种算法分别对仿真数据集进行三维成像的各项指标评估结果
| 方法/指标 | RP | MD | AC | CO |
| MUSIC | 0.1579 | 46.1899 | 5.8264 | 1.6841 |
| OMP | 0.3633 | 17.4431 | 4.9466 | 1.1192 |
| GOMP | 0.3635 | 17.2847 | 4.9047 | 1.1178 |
| 本发明实施例(ω=3) | 0.5146 | 3.2984 | 2.4684 | 1.1281 |
| 本发明实施例(ω=5) | 0.5627 | 3.0201 | 1.7727 | 1.1341 |
| 本发明实施例(ω=7) | 0.5690 | 3.0059 | 1.7138 | 1.1461 |
| 本发明实施例(ω=9) | 0.5685 | 3.0265 | 1.7329 | 1.1564 |
采用MUSIC算法、OMP算法、GOMP算法和本发明实施例中的方法（）时的进行成
像的结果如图7所示。从成像结果可以直观的看出，本发明实施例中的方法，恢复的结果只
有少数噪声点，在平地，平缓斜坡和陡峭斜坡出恢复效果均较好，其他算的则在斜坡变化区
域则出现了不同程度的噪点。
具体实施例6
如图8所示，本发明实施例6提供的一种基于同质区域联合的层析SAR 三维成像装置，包括非易失性存储器和一个或多个处理器，所述非易失性存储器中存储有可执行代码，所述处理器执行所述可执行代码时，用于实现上述实施例中的基于同质区域联合的层析SAR 三维成像方法。
本发明基于同质区域联合的层析SAR 三维成像装置的实施例可以应用在任意具备数据处理能力的设备上，该任意具备数据处理能力的设备可以为诸如计算机等设备或装置。装置实施例可以通过软件实现，也可以通过硬件或者软硬件结合的方式实现。以软件实现为例，作为一个逻辑意义上的装置，是通过其所在任意具备数据处理能力的设备的处理器将非易失性存储器中对应的计算机程序指令读取到内存中运行形成的。从硬件层面而言，如图8所示，为本发明基于同质区域联合的层析SAR 三维成像装置所在任意具备数据处理能力的设备的一种硬件结构图，除了图8所示的处理器、内存、网络接口、以及非易失性存储器之外，实施例中装置所在的任意具备数据处理能力的设备通常根据该任意具备数据处理能力的设备的实际功能，还可以包括其他硬件，对此不再赘述。
上述装置中各个单元的功能和作用的实现过程具体详见上述方法中对应步骤的实现过程，在此不再赘述。
对于装置实施例而言，由于其基本对应于方法实施例，所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本发明方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
具体实施例7
本发明实施例还提供一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时，实现上述实施例中的基于同质区域联合的层析SAR 三维成像方法。
所述计算机可读存储介质可以是前述任一实施例所述的任意具备数据处理能力的设备的内部存储单元，例如硬盘或内存。所述计算机可读存储介质也可以是任意具备数据处理能力的设备的外部存储设备，例如所述设备上配备的插接式硬盘、智能存储卡（Smart Media Card，SMC）、SD卡、闪存卡（Flash Card）等。进一步的，所述计算机可读存储介质还可以既包括任意具备数据处理能力的设备的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述任意具备数据处理能力的设备所需的其他程序和数据，还可以用于暂时地存储已经输出或者将要输出的数据。
综上，本发明方法设计了一种结合欧氏距离、像素强度和边缘信息的测度函数，将SAR图像中三维位置相近、电磁散射强度相似的像素点分割成局部连通的同质区域。然后，将TomoSAR单像素观测信号模型扩展到同质区域联合观测信号模型，通过同质区域积累得到高程约束中心，根据同质区域大小和成像参数自动确定高程搜索范围，稀疏逼近估计同质区域高程。最后，对斜高进行几何变换和地理编码，得到三维成像结果。该方法能减少成像结果中由固有相干斑噪声带来的噪点，有效提升了信噪比，提高了对弱小目标的恢复能力。
Claims (8)
1. 一种基于同质区域联合的层析SAR 三维成像方法，其特征在于，包括以下步骤：
S1，对SAR图像进行同质区域分割；
S2，进行同质区域的联合稀疏恢复，计算同质区域所有像素斜高估计值；
S3，进行几何坐标变换，将斜高向转换为实际空间的高度向，得到三维成像结果；
所述步骤S1中对SAR图像进行同质区域分割的步骤为：
S101，在SAR图像上均匀初始化尺寸为像素的/>个网格聚类中心，将聚类中心以元素/>的形式存入优先级队列/>，其中/>、/>为大于1的整数，/>为空间位置，/>为切片的强度向量，/>为边缘强度信息，元素/>的第四列为/>，/>表示第/>个待处理像素点到第/>个聚类中心的相似度，初始化聚类中心的/>为0，初始化与SAR图像尺寸相同的标签图；
S102，计算聚类中心四邻域像素到所述聚类中心自身的相似度，并以元素的形式推入优先级队列/>，每一次向优先级队列/>中推入一个元素，更新一次排序，使优先级队列/>中元素保持从大到小的顺序，推出优先级队列/>的顶端元素/>，所述顶端元素/>为与所述聚类中心相似度最大的像素，将所述聚类中心的标签赋给顶端元素/>对应的像素，用顶端元素/>的空间位置/>和切片强度向量/>更新所述聚类中心的空间位置和切片强度；
S103，依次计算顶端元素的四邻域像素点到聚类中心的相似度，创建相应的元素/>，若所述顶端元素/>的四邻域像素点已被标记，则已被标记的像素点所创建的元素不推入优先级队列/>，只将未标记的像素点创建的元素推入优先级队列/>，每一次向优先级队列/>中推入一个元素，更新一次排序，使优先级队列/>中元素保持从大到小的顺序，推出优先级队列/>的顶端元素/>，将所述聚类中心的标签赋给顶端元素/>对应的像素，用顶端元素/>的空间位置/>和切片强度向量/>更新所述聚类中心的空间位置和切片强度；
S104，若所有的像素点均被标记，并且优先级队列为空，输出同质区域聚类标签图，否则返回步骤S103；
所述相似度按照以下公式计算：
，
其中，为第/>个待处理像素点到第/>个聚类中心的相似度，/>为第/>个待处理像素点到第/>个聚类中心的空间相似度，/>为第/>个待处理像素点到第/>个聚类中心的像素强度相似度，/>为第/>个待处理像素点的边缘强度，/>、/>为加权系数；
所述空间相似度按照以下公式计算：/>，
其中，欧氏距离，/>、/>为第/>像素点的坐标，/>、/>为第个聚类中心的坐标；
所述像素强度相似度按照以下公式计算：/>，
其中，是像素强度比值，/>，/>为待处理像素切片的强度，/>为聚类中心切片的强度，/>是/>的概率密度函数，/>表示切片所包含的像素点个数，/>为标准高斯核函数，m表示切片所包含的像素点的序号，/>；所述概率密度函数/>按照以下公式计算：/>，其中，/>为gamma分布函数，/>为成像视数。
2. 根据权利要求1所述的基于同质区域联合的层析SAR 三维成像方法，其特征在于，所述加权系数的取值范围为[2, 6]，加权系数/>的取值范围为[0.3, 0.7]。
3. 根据权利要求1所述的基于同质区域联合的层析SAR 三维成像方法，其特征在于，所述的取值为5或7。
4. 根据权利要求1所述的基于同质区域联合的层析SAR 三维成像方法，其特征在于，所述步骤S2中进行同质区域的联合稀疏恢复，计算同质区域所有像素斜高估计值的步骤为：
S201，初始化斜高估计值，/>，稀疏集/>，残差/>，其中，/>为观测信号向量，/>，/>为同质区域内部的第/>个像素的观测信号，，/>为同质区域内部的像素个数，/>，N为相干观测的个数；
S202，构建同质区域内的观测矩阵，计算观测矩阵/>与观测信号向量/>的相关性向量/>，其中，/>为同质区域内的第/>个像素的观测矩阵，，/>为第/>个相干观测的等效空间频率，/>，/>为第/>个预先划分的高程位置，/>为高程向上划分的网格数量；
S203，根据相关性向量主峰数确定同质区域的叠掩数/>；设置k =1；
S204，计算观测矩阵与残差/>的相关性向量/>；
S205，由相关性向量的第一主峰确定约束中心；
S206，在约束中心的范围内搜索每个像素的极大值/>；
S207，将极大值对应的观测矩阵/>中的原子纳入同质区域重构稀疏基/>，
，其中，/>；
S208，计算该稀疏基下的斜高估计值，/>，其中， ， />为共轭转置；
S209，计算该稀疏基下的残差，/>，其中，/>；
S210，；
S211，若则结束并输出/>，否则返回步骤S204。
5. 根据权利要求4所述的基于同质区域联合的层析SAR 三维成像方法，其特征在于，所述步骤S202中所述计算观测矩阵与观测信号向量/>的相关性向量/>的步骤为：
S2021，计算同质区域内部每个高程位置与观测信号的相关性，，其中，/>为同质区域内部第/>个像素在/>高程位置与观测信号/>的相关性；
S2022，将不同像素同一高度位置的相关性进行非相干叠加，拟合整个同质区域的相关性向量，/> ；
所述步骤S204中所述计算观测矩阵与残差/>的相关性向量/>的步骤为：
S2041，计算同质区域内部每个高程位置与残差的相关性/>，，其中，/>为同质区域内部第/>个像素在/>高程位置与残差/>的相关性；
S2042，将不同像素同一高度位置的相关性进行非相干叠加，拟合整个同质区域的相关性向量，/> 。
6. 根据权利要求4所述的基于同质区域联合的层析SAR 三维成像方法，其特征在于，步骤S206中所述约束中心的范围计算方法为：
，其中，/>为同质区域网格长度，/>为下视角，/>为距离向像素尺寸。
7. 一种基于同质区域联合的层析SAR 三维成像装置，包括非易失性存储器和一个或多个处理器，所述非易失性存储器中存储有可执行代码，所述处理器执行所述可执行代码时，用于实现如权利要求1-6中任一项所述的基于同质区域联合的层析SAR 三维成像方法。
8. 一种计算机可读存储介质，其特征在于，所述计算机可读存储介质包括存储的计算机程序，其中，在所述计算机程序被处理器运行时控制所述存储介质所在设备执行如权利要求1-6中任一项所述的基于同质区域联合的层析SAR 三维成像方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310665210.1A CN116400356B (zh) | 2023-06-07 | 2023-06-07 | 一种基于同质区域联合的层析sar三维成像方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310665210.1A CN116400356B (zh) | 2023-06-07 | 2023-06-07 | 一种基于同质区域联合的层析sar三维成像方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116400356A CN116400356A (zh) | 2023-07-07 |
| CN116400356B true CN116400356B (zh) | 2023-08-18 |
Family
ID=87012696
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310665210.1A Active CN116400356B (zh) | 2023-06-07 | 2023-06-07 | 一种基于同质区域联合的层析sar三维成像方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116400356B (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011070538A2 (en) * | 2009-12-11 | 2011-06-16 | Eads Singapore Pte. Ltd | Method for despeckling of single-look dual-polarization synthetic aperture radar (sar) data |
| CN102662171A (zh) * | 2012-04-23 | 2012-09-12 | 电子科技大学 | 一种sar层析三维成像方法 |
| CN103969645A (zh) * | 2014-05-14 | 2014-08-06 | 中国科学院电子学研究所 | 基于压缩多信号分类的层析合成孔径雷达测量树高的方法 |
| CN105388476A (zh) * | 2015-12-28 | 2016-03-09 | 河南工业大学 | 一种基于联合稀疏模型的层析sar成像方法 |
| CN110109101A (zh) * | 2019-04-04 | 2019-08-09 | 电子科技大学 | 一种基于自适应阈值的压缩感知三维sar成像方法 |
| CN111948654A (zh) * | 2020-08-12 | 2020-11-17 | 中国科学院空天信息创新研究院 | 机载层析sar三维点云生成方法 |
| CN114004998A (zh) * | 2021-11-03 | 2022-02-01 | 中国人民解放军国防科技大学 | 基于多视张量积扩散的非监督极化sar图像地物分类方法 |
| CN114359744A (zh) * | 2021-12-07 | 2022-04-15 | 中山大学 | 一种基于激光雷达和事件相机融合的深度估计方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102496153B (zh) * | 2011-11-04 | 2013-08-14 | 西安电子科技大学 | 基于小波域中字典学习的sar图像相干斑抑制方法 |
| US10198858B2 (en) * | 2017-03-27 | 2019-02-05 | 3Dflow Srl | Method for 3D modelling based on structure from motion processing of sparse 2D images |
| CN108682017B (zh) * | 2018-04-11 | 2021-06-18 | 浙江工业大学 | 基于Node2Vec算法的超像素图像边缘检测方法 |
-
2023
- 2023-06-07 CN CN202310665210.1A patent/CN116400356B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011070538A2 (en) * | 2009-12-11 | 2011-06-16 | Eads Singapore Pte. Ltd | Method for despeckling of single-look dual-polarization synthetic aperture radar (sar) data |
| CN102662171A (zh) * | 2012-04-23 | 2012-09-12 | 电子科技大学 | 一种sar层析三维成像方法 |
| CN103969645A (zh) * | 2014-05-14 | 2014-08-06 | 中国科学院电子学研究所 | 基于压缩多信号分类的层析合成孔径雷达测量树高的方法 |
| CN105388476A (zh) * | 2015-12-28 | 2016-03-09 | 河南工业大学 | 一种基于联合稀疏模型的层析sar成像方法 |
| CN110109101A (zh) * | 2019-04-04 | 2019-08-09 | 电子科技大学 | 一种基于自适应阈值的压缩感知三维sar成像方法 |
| CN111948654A (zh) * | 2020-08-12 | 2020-11-17 | 中国科学院空天信息创新研究院 | 机载层析sar三维点云生成方法 |
| CN114004998A (zh) * | 2021-11-03 | 2022-02-01 | 中国人民解放军国防科技大学 | 基于多视张量积扩散的非监督极化sar图像地物分类方法 |
| CN114359744A (zh) * | 2021-12-07 | 2022-04-15 | 中山大学 | 一种基于激光雷达和事件相机融合的深度估计方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于广义似然比不规则建筑物层析SAR成像方法;周鹏 等;智能处理与应用(第6期);84-92 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116400356A (zh) | 2023-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| de Conto et al. | Performance of stem denoising and stem modelling algorithms on single tree point clouds from terrestrial laser scanning | |
| Çetin et al. | Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization | |
| Zheng et al. | Unsupervised saliency-guided SAR image change detection | |
| CN102662171B (zh) | 一种sar层析三维成像方法 | |
| Li et al. | Refraction corrected transmission ultrasound computed tomography for application in breast imaging | |
| Malmgren-Hansen et al. | Convolutional neural networks for SAR image segmentation | |
| Wu et al. | 3D Tree reconstruction from simulated small footprint waveform lidar | |
| CN104251991B (zh) | 一种基于稀疏度估计的分维度阈值迭代稀疏微波成像方法 | |
| CN116012364B (zh) | Sar图像变化检测方法和装置 | |
| Atkinson et al. | Multi-resolution soil-landscape characterisation in KwaZulu Natal: Using geomorphons to classify local soilscapes for improved digital geomorphological modelling | |
| CN112415515B (zh) | 一种机载圆迹sar对不同高度目标分离的方法 | |
| KR102262397B1 (ko) | 다중 시기에 획득된 위성 sar 영상 간의 정합 자동화 장치 및 방법 | |
| Rambour et al. | Urban surface reconstruction in SAR tomography by graph-cuts | |
| CN108509835B (zh) | 基于DFIC超像素的PolSAR图像地物分类方法 | |
| CN114509758A (zh) | 一种基于分形维数和极化分解的月表平坦区域选取方法 | |
| Huang et al. | Change detection method based on fractal model and wavelet transform for multitemporal SAR images | |
| Stuart et al. | Iceberg size and orientation estimation using SeaWinds | |
| D'Hondt et al. | Geometric primitive extraction for 3D reconstruction of urban areas from tomographic SAR data | |
| CN113608218A (zh) | 一种基于后向投影原理的频域干涉相位稀疏重构方法 | |
| CN116400356B (zh) | 一种基于同质区域联合的层析sar三维成像方法 | |
| CN109959933B (zh) | 一种基于压缩感知的多基线圆迹合成孔径雷达成像方法 | |
| Wu et al. | ISAR Image Registration Based on Normalized Correlation Coefficient | |
| Zhu et al. | Global LoD-1 building model from TanDEM-X data | |
| US20130332111A1 (en) | Generating information conditional on mapped measurements | |
| Yang et al. | AAE-Dpeak-SC: A novel unsupervised clustering method for space target ISAR images based on adversarial autoencoder and density peak-spectral clustering |
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 |