CN114723794B - 基于lsd直线检测的sar图像配准方法 - Google Patents

基于lsd直线检测的sar图像配准方法 Download PDF

Info

Publication number
CN114723794B
CN114723794B CN202210398024.1A CN202210398024A CN114723794B CN 114723794 B CN114723794 B CN 114723794B CN 202210398024 A CN202210398024 A CN 202210398024A CN 114723794 B CN114723794 B CN 114723794B
Authority
CN
China
Prior art keywords
image
sar
calculating
digital map
sar image
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
CN202210398024.1A
Other languages
English (en)
Other versions
CN114723794A (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 Thunderbolt Information Technology Co.,Ltd.
Original Assignee
Nanjing Leading Information Technology Co ltd
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 Leading Information Technology Co ltd filed Critical Nanjing Leading Information Technology Co ltd
Priority to CN202210398024.1A priority Critical patent/CN114723794B/zh
Publication of CN114723794A publication Critical patent/CN114723794A/zh
Application granted granted Critical
Publication of CN114723794B publication Critical patent/CN114723794B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了基于LSD（Line Segment Detector）直线检测的SAR图像配准方法，包括：读取SAR图像与卫星数字地图，构建待配准图像对；对待配准图像对进行预处理；采用LSD直线检测方法检测预处理图像结果中的直线特征，以点数据的形式进行存储；计算出角度参数、尺度参数、变换矩阵和平移参数，然后通过尺度变换、平移变换以及变换矩阵的处理，得到SAR图像中的目标点在卫星数字地图中对应的匹配点；载入卫星数字地图配套的经纬度数据库，根据匹配点坐标获取目标点经纬度信息。本发明通过使用基于线特征对SAR图像与卫星数字地图进行配准工作，提升了配准的准确度，对视角和光照变化具有很好的鲁棒性。

Description

基于LSD直线检测的SAR图像配准方法
技术领域
本发明涉及SAR图像配准技术领域，具体涉及基于LSD直线检测的SAR图像配准方法。
背景技术
合成孔径雷达(Synthetic Aperture Radar，SAR)成像具有全天时、全天候、可穿透性等优点，以及它在不同频段、不同极化下可得到目标的高分辨率雷达图像，为人们提供非常有用的目标信息。因此，一直以来针对SAR图像配准的研究热度经久不息，大量研究成果已经被广泛应用于经济和科技等众多领域，有着广泛的应用前景和发展潜力。
图像配准技术是将相同地区，在不同时刻、不同视角、不同传感器或不同光照条件下拍摄的图像进行空间对准的过程。由于采集位置、空间分辨率、采集时间等不同，使得获取的同一场影的光学与SAR影像间存在较大的辐射与几何差异。因此，要综合利用不同遥感影像必须使它们实现几何意义上的对准，即需要进行图像配准。
图像配准方法一般分为两种：一种是基于灰度的方法，另一种是基于特征的方法。基于灰度的方法使用图像或者图像块的灰度信息，通过相似性度量算法，测度两幅对应图像或对应图像块之间的相似性，搜索参数空间，寻找相似性最优的变换模型的参数值，从而实现图像配准，但是因为其全局优化的计算复杂度太高是不太适用于对时间要求比较严格的配准的。基于特征的方法用对图像中具有稳定性和重复性图像特征的分析，代替对整个图像的分析，降低计算量以及配准算法对传感器变化的敏感性。
对比基于灰度的方法，基于特征的方法因为具有更好的重复性和稳定性，所以它具有更强大的区分能力，对局部有差异的图像对有更好的鲁棒性。而且基于特征的方法使用少量图像特征，对相同图像来说，算法的运算性能普遍高于基于灰度的方法。同时，因为SAR图像分辨率高、数据量大、受噪声影响严重的问题，现有的配准算法处理速度很慢，所以为了高效性和高精度的合成孔径雷达图像的配准，对原有的配准算法改进使其具有高效性和高精度特性是一种比较普遍存在的方式。因为基于特征方法相较于基于灰度方法具有计算量低的特征，所以基于特征的方法是当前SAR图像使用最多的配准方法。但是，由于合成孔径雷达成像方式与可见光成像方式的不同，造成两种不同图像在点特征上的差异，即二者之间能检测出来的角点并不完全是一一对应的，因此基于点特征配准方法降低了匹配精度。
发明内容
发明目的：本发明目的在于针对现有技术的不足，提供一种基于LSD(LineSegment Detector)直线检测的SAR图像配准方法，与有重叠区域的可见光图像进行比对，减弱或者去除因为传感器或者天气变换等原因导致的图像变化的影响，获取更加精准的目标信息。
技术方案：本发明所述基于LSD直线检测的SAR图像配准方法，包括如下步骤：
S1：读取SAR图像与卫星数字地图，构建待配准图像对Msar、Msate；
S2：对待配准图像对Msar、Msate进行预处理，得到预处理后的图像Sarcanny、Satecanny；
S3：采用LSD直线检测方法检测S2预处理后的图像Sarcanny、Satecanny中的直线特征，保存检测到的直线信息，以点数据的形式进行存储，并表示为Vsar、Vsate；
S4：根据Vsar、Vsate计算出角度参数theta、尺度参数s、变换矩阵transMatric和平移参数，然后对图像Msar进行尺度变换、平移变换和变换矩阵处理，得到处理后图像中的目标点在卫星数字地图Msate中对应的匹配点，表示为MP(x，y)；
S5：载入卫星数字地图Msate配套的经纬度数据库，根据匹配点坐标MP(x，y)获取目标点经纬度信息，表示为Lon、Lat。
进一步完善上述技术方案，所述S1中SAR图像与卫星数字地图以多维矩阵的格式进行存储。
进一步地，所述S2中对Msate进行的预处理操作包括灰度化、滤波、二值化、形态学运算、高斯滤波、Canny边缘检测，对Msar进行的预处理操作包括灰度化、滤波、二值化、形态学运算、相干斑抑制、高斯滤波、Canny边缘检测。
进一步地，所述预处理操作过程如下：
S21：将待配准图像对进行灰度化处理
设r、g、b分别表示待处理图像的三个通道，采用平均值法进行灰度化处理，将三通道转换为单通道，得到图中(i，j)处的像素灰度值为：
grayVal_mean(i，j)＝(r(i，j)+g(i，j)+b(i，j))/3 (1)
r(i，j)、g(i，j)、b(i，j)表示三个通道r、g、b在图中(i，j)处的像素灰度值；
将待配准图像对Msar、Msate按公式(1)处理，得到SAR图像与卫星数字地图的灰度图；
S22：对灰度化的待配准图像对进行滤波处理，然后对滤波结果进行二值化处理，再根据二值化图像中非零像素个数的占比重复滤波及二值化操作，迭代至非零像素个数占比达到预设阈值，包括：
对S21中灰度图进行归一化处理，使像素灰度值在0～1范围内，设max[grayVal(i，j)]为图像中最大灰度值，min[grayVal(i，j)]为图像中最小灰度值，通过下式进行图像归一化：
Figure BDA0003592648550000031
设置卷积核kernel的大小，以核内中心位置对应的像素为目标像素，使用kernel对待配准图像对进行卷积操作，卷积结果为：
Figure BDA0003592648550000032
按从左到右、从上到下滑动卷积核，直到处理完图像最后一个像素，
对卷积结果进行二值化处理，得到二值化后的结果Sarthrehold、Satethrehold，凸显待配准图像中关键区域的轮廓边缘信息；
S23：对Sarthrehold、Satethrehold进行形态学处理，将断裂区域连接起来，得到处理结果Sardilate、Satedilate；
S24：根据SAR图像Msar的距离向分辨率Dre与方位向分辨率Dis构建滤波核kernelnoise，大小为(9/Dis)×(9/Dre)，以核中心点为原点，通过下式构建：
Figure BDA0003592648550000033
相干斑抑制步骤如下：
S241、采用构建的滤波核在二值SAR图像Sardilate区域进行滑动操作，从左到右，从上到下，
S242、判断滤波核中元素与对应二值SAR图像Sardilate中元素是否相等且等于255，记录符合条件的元素个数，记为Pixel255num，
S243、判断滤波核中元素值为255，且对应二值SAR图像Sardilate中元素值为0，记录符合条件的元素个数，记为Pixel0/255num，
S244、当Pixel255num与Pixel0/255num的值大于滤波核元素个数的1/3时，判断图像中每行左右邻域或每列上下邻域中像素灰度值为255的像素个数是否超过5，如果超过则设置目标像素的灰度值为0，
S245、重复S241～S244，直到滤波核中心位置与图像右下角像素重合，得到处理结果Sarcoh；
S25：对卫星数字地图Satedilate和相干斑抑制后的SAR图像Sarcoh进行高斯滤波去噪、Canny边缘检测，得到处理后的图像Sarcanny、Satecanny。
进一步地，所述S3采用LSD直线检测方法检测图像Sarcanny与Satecanny中的直线特征，包括：
S31：首先计算图像中所有像素点的梯度大小，记为grd，根据梯度大小grd计算出水平角度agl与梯度幅值|grd|，如下式所示：
Figure BDA0003592648550000041
agl＝arctan(grdi/-grdj) (6)
Figure BDA0003592648550000042
利用伪排序得到梯度幅值大的点作为种子点，以该点水平线角度agl作为区域的初始角度θr，在八邻域中寻找与θr的偏差小于预设阈值threshold的点，将该点加入到区域中并更新θr，得到线支持区域，最后通过矩形逼近得到直线特征，检测到的直线特征以点数据的形式进行存储，分别表示为Psar、Psate；
S32：筛选Psar、Psate中最具代表性的直线特征来进行配准工作，筛选步骤如下：
S321、计算Psar、Psate中所有直线段的长度，并按降序排列，表示为Lsar、Lsate，
S322、根据预设直线特征数量Lnum，在Lsar、Lsate中提取前Lnum条直线特征，
S323、提取距离间隔较大的直线特征，
S324、保存直线特征的端点坐标，表示为Vsar、Vsate。
进一步地，所述S4包括：
S41：将直线特征的端点坐标Vsar、Vsate转化为极坐标，表示为Pcsar(angle，solar)、Pcsate(angle，solar)，
S411、计算每条直线段的斜率k与截距b；
S412、计算直线段的中心点坐标(xmid，ymid)、kmid及bmid
Figure BDA0003592648550000051
kmid＝-1/k (9)
bmid＝ymid-kmid×xmid (10)
S413、计算极坐标的极角与极径：
当kmid＞0时，获得如下三个点：
O(x3，y3)＝O(xmid，ymid)
A(x4，y4)＝A[xmid-1，kmid×(xmid-1)+bmid]
B(x5，y5)＝B(xmid-1，ymid)
当kmid＜0时，获得如下三个点：
O(x3，y3)＝O(xmid，ymid)
A(x4，y4)＝A[xmid+1，kmid×(xmid+1)+bmid]
B(x5，y5)＝B(xmid-1，ymid)
得到极角angle：
Figure BDA0003592648550000052
Figure BDA0003592648550000053
dot＝(x3-x4)×(x3-x5)+(y3-y4)×(y3-y5)
angle＝[(a cos(dot/(m1×m2)))×180]/3.1415926
得到极径solar：
令
Figure BDA0003592648550000054
则有
Figure BDA0003592648550000055
S42：以卫星数字地图为基准，计算SAR图像中直线与卫星数字地图中直线的角度差theta.
计算Pcsar与Pcsate中angle的差值：angle_diff＝Pcsar(angle)-Pcsate(angle)，将相同的angle-diff值进行计数，出现次数最多的angle_diff值记为theta；
S43：计算尺度参数s，计算步骤如下：
S431、采用Pcsar中的angle统计出SAR图像中平行直线组；
S432、计算SAR图像与卫星数字地图中每组平行直线极坐标的极径差
solarsar＝abs(Pcsar(solarn)-Pcsar(solarn-1))
solarsate＝abs(Pcsate(solarn)-Pcsate(solarn-1))
scale＝solarsate/solarsar
将相同的scale值进行计数，将出现次数最多的scale值记为尺度参数s；
S44、计算变换矩阵transMatric与平移参数
S441、分别计算SAR图像与卫星数字地图中直线特征的交点，分别表示为nodesar、nodesate，
S442、分别计算交点集合nodesar与nodesate中各点间坐标的差分node_diff(xdiff，ydiff)，
S443、将node_diff中的横坐标差值与纵坐标差值分别与阈值比较，保存满足条件的node_diff，同时将相同的node_diff值进行计数；
S444、获取出现次数最多的node_diff中的xdiff、ydiff；
S445、定义与SAR图像尺寸相关参数imgSarx与imgSary，与卫星数字地图尺寸相关参数imgSatex与imgSatey
Figure BDA0003592648550000061
Figure BDA0003592648550000062
S446、计算变换矩阵
Figure BDA0003592648550000063
Figure BDA0003592648550000064
S447、计算平移参数
Figure BDA0003592648550000071
S45、结合尺度参数s、变换矩阵transMatric、平移参数计算出SAR图像在卫星数字地图中的匹配点MP(x，y)，设SAR图像中的目标点为C(x，y)，则得到匹配点MP(x，y)：
Figure BDA0003592648550000072
有益效果：与现有技术相比，本发明的优点在于：本发明通过使用基于线特征对SAR图像与卫星数字地图进行配准工作，提升了配准的准确度，相比于点特征配准方法，对视角和光照变化具有很好的鲁棒性，克服了因不同传感器成像方式差异带来的误差；在特征提取方面，本发明在效率上也高于现有的基于点特征匹配方法；与现有的基于Hough直线特征检测的配准技术相比，本发明具有结果准确、误检可控，同时不需要调节参数的优点。
附图说明
图1是本发明的方法流程图；
图2是本发明中滤波卷积核的示意图；
图3是本发明中待卷积图像区域；
图4是本发明中待配准图像对；
图5是本发明中配准后的效果图；
图6是本发明中SAR图像与卫星数字地图匹配结果；
图7是本发明SAR图像中已匹配目标的经纬度信息结果。
具体实施方式
下面通过附图对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
本发明采用基于LSD(Line Segment Detector)直线检测的SAR图像配准方法：首先SAR图像与卫星数字地图对进行灰度化、滤波、二值化等预处理，然后进行相干斑抑制、高斯滤波、以及边缘检测，获取可以进行直线检测的图像样本，接着使用由直线特征计算得到的尺度参数、角度参数以及变换矩形获取SAR图像中目标点在卫星数字地图中的匹配点，最后根据卫星数字地图中匹配点坐标与卫星数字地图配套经纬度数据库计算SAR图像中目标点的地理位置信息。
如图1所示的基于LSD直线检测的SAR图像配准方法，包括如下步骤：
S1：读取SAR图像与卫星数字地图，使用多通道矩阵存储，构建待配准图像对，并表示为Msar、Msate；
S2：对待配准图像对进行预处理，主要包括灰度化、滤波、二值化、形态学运算，相干斑抑制，高斯滤波、Canny边缘检测；
S3：对S2中的图像结果使用LSD直线检测方法检测图像结果中的直线特征，保存从SAR图像中与卫星数字地图中检测到的直线信息，以点数据的形式进行存储，分别表示为Vsar、Vsate；
S4：根据上述Vsar和Vsate计算出尺度参数、角度参数及变换矩阵，分别表示为s，theta，transMatric，然后通过尺度变换、平移变换以及变换矩阵的处理，得到SAR图像中的目标点在卫星数字地图中对应的匹配点，表示为Mp(x，y)；
S5：载入经纬度数据库，根据匹配点坐标MP(x，y)计算目标点经纬度信息，表示为Lon、Lat。
具体地，S1中使用的SAR图像与卫星数字地图，以多维矩阵的格式进行存储，读取图像关键像素信息。
S2中对图像进行预处理的操作，包括灰度化、滤波、二值化、形态学运算，相干斑抑制(只有SAR图像需要)，高斯滤波、Canny边缘检测，具体地，包括：
S21：设r、g、b分别是SAR图像与卫星数字地图的三个通道，采用平均值法进行灰度化处理，将三通道转换为单通道，即将SAR图像与卫星数字地图的三个通道r、g、b在(i，j)处的像素灰度值r(i，j)、g(i，j)、b(i，j)相加，求平均值grayVal_mean(i，j)，得到在图中(i，j)处的像素灰度值为：
grayVal_mean(i，j)＝(r(i，j)+g(i，j)+b(i，j))/3 (1)
最终将待配准图像对按公式1处理，得到SAR图像与卫星数字地图的灰度图。
S22：对待配准图像进行滤波处理，提高对比度，突显图像边缘等特征，然后对滤波结果进行二值化处理，再根据二值化图像中非零像素个数的占比重复滤波与二值化操作，迭代至非零像素个数占比达到预设阈值，包括
对S21中灰度图进行归一化处理，使像素灰度值在0～1范围内，设max[grayVal(i，j)]为图像中最大灰度值，min[grayVal(i，j)]为图像中最小灰度值，通过下式进行图像归一化：
Figure BDA0003592648550000081
设卷积核kernel的大小为3×3(kernel的大小根据实际情况可以设置)，如图2所示。以核内中心位置对应的像素为目标像素，使用kernel对待配准图像对进行卷积操作，对应待卷积的图像区域如图3所示，卷积结果为：
Figure BDA0003592648550000091
按从左到右，从上到下滑动卷积核，直到处理完图像最后一个像素。
接下来，对卷积结果进行二值化处理，得到二值化后的结果Sarthrehold，Satethrehold凸显待配准图像中关键区域的轮廓等边缘信息。
S23：对Sarthrehold与Satethrehold进行形态学处理，将断裂区域连接起来，得到处理结果Sardilate、Satedilate；
S24：对二值化后的SAR图像Sardilate进行滤波处理、，抑制相干斑噪声
根据SAR图像的距离向分辨率与方位向分辨率(分别表示为Dre，Dis)构建滤波核kernelnoise，大小为(9/Dis)×(9/Dre)，以核中心点为原点，通过下式构建：
Figure BDA0003592648550000092
相干斑噪声抑制步骤如下：
1)采用构建的滤波核在SAR图像Sardilate区域进行滑动操作，从左到右，从上到下；
2)判断滤波核中元素与对应二值SAR图像元素是否相等且等于255，记录符合条件的元素个数，记为Pixel255num；
3)判断滤波核中元素值为255，且对应二值SAR图像元素值为0，记录符合条件的元素个数，记为Pixel0/255num；
4)当Pixel255hum与Pixel0/255num的值大于滤波核元素个数的1/3时，判断图像中每行左右邻域或每列上下邻域中像素灰度值为255的像素个数是否超过5，如果超过则设置目标像素的灰度值为0；
5)重复上述操作，直到滤波核中心位置与图像右下角像素重合，得到结果Sarcoh。
S25：对卫星数字地图Satedilate与相干斑抑制后的SAR图像Sarcoh进行高斯滤波去噪、Canny边缘检测，得到处理后的图像Sarcanny、Satecanny。
S3采用LSD直线检测方法能在线性的时间内得出亚像素级精度的检测结果，且无需调节任何参数。
S31：首先计算图像Sarcanny、Satecanny中所有像素点的梯度大小和水平角度，分别记为grd、agl，根据梯度大小计算出agl与梯度幅值|grd|，如下式所示：
Figure BDA0003592648550000101
/>
agl＝arctan(grdi/-grdj) (6)
Figure BDA0003592648550000102
利用伪排序得到梯度幅值大的点作为种子点，以该点水平线角度agl作为区域的初始角度θr，在八邻域中寻找与θr的偏差小于预设阈值threshold的点，将该点加入到区域中并更新θr，得到线支持区域，最后通过矩形逼近得到直线特征，检测到的直线特征以点数据的形式进行存储，分别表示为Psar、Psate。
S32：筛选Psar、Psate中最具代表性的直线特征用来进行配准工作，Psar、Psate中存储了检测到的所有直线段，筛选步骤如下：
1)计算所有直线段的长度，并按降序排列，表示为Lsar、Lsate；
2)根据预设直线特征数量Lnum，在Lsar、Lsate中提取前Lnum条直线特征；
3)提取距离间隔较大的直线特征；
保存直线特征的端点坐标，表示为Vsar、Vsate。
S4中根据上述Vsar、Vsate计算出尺度参数、角度参数及变换矩阵，分别表示为s、theta、transMatric；然后通过尺度变换、平移以及变换矩阵的处理，得到SAR图像中的目标点在卫星数字地图中对应的匹配点，表示为MP(x，y)。
S41：将Vsar、Vsate中直线特征的端点坐标转化为极坐标，表示为PCsar(angle，solar)、Pcsate(angle，solar)
1)计算每条直线段的斜率k与截距b；
2)计算直线段的中心点坐标(xmid，ymid)、kmid及bmid
Figure BDA0003592648550000111
kmid＝-1/k (9)
bmid＝ymid-kmid×xmid (10)
计算极坐标的极角与极径：
当kmid＞0时，由2中可以获得如下三个点：
O(x3，y3)＝O(xmid，ymid)
A(x4，y4)＝A[xmid-1，kmid×(xmid-1)+bmid]
B(x5，y5)＝B(xmid-1，ymid)
当kmid＜0时，由2中可以获得如下三个点：
O(x3，y3)＝O(xmid，ymid)
A(x4，y4)＝A[xmid+1，kmid×(xmid+1)+bmid]
B(x5，y5)＝B(xmid-1，ymid)
可以得到极角angle：
Figure BDA0003592648550000112
Figure BDA0003592648550000113
dot＝(x3-x4)×(x3-x5)+(y3-y4)×(y3-y5)
angle＝[(a cos(dot/(m1×m2)))×180]/3.1415926
可以得到极径solar：
令
Figure BDA0003592648550000114
则有
Figure BDA0003592648550000115
S42：以卫星数字地图为基准，计算SAR图像中直线与卫星数字地图中直线的角度差theta。
根据在SAR图像与卫星数字地图中检测到的直线特征计算出SAR图像与卫星数字地图之间的旋转角度，S41中通过直线端点坐标将直线用极坐标，分别表示为Pcsar(angle，solar)与Pcsate(angle，solar)；
计算Pcsar与Pcsate中angle的差值：
angle_diff＝Pcsar(angle)-Pcsate(angle)
将相同的angle_diff值进行计数，出现次数最多的angle_diff值即为theta；
S43：计算尺度参数s
尺度参数计算步骤如下：
(1)采用Pcsar中的angle统计出SAR图像中平行直线组
(2)计算SAR图像与卫星数字地图中每组平行直线极坐标的极径差
solarsar＝abs(Pcsar(solarn)-Pcsar(solarn-1))
solarsate＝abs(Pcsate(solarn)-Pcsate(solarn-1))
scale＝solarsate/solarsar
将相同的scale值进行计数，出现次数最多的scale值即为尺度参数s；
使用s、theta值可以将SAR图像的尺度与朝向变换为与卫星数字地图一致，如图4、图5所示，其中图4为原图，图5为处理结果。
S44：计算变换矩阵transMatric与平移参数
(1)分别计算SAR图像与卫星数字地图中直线特征的交点，分别表示为nodesar与nodesate；
(2)分别计算交点集合nodesar与nodesate中各点间坐标的差分node_diff(xdiff，ydiff)；
(3)将node_diff中的横坐标差值与纵坐标差值分别与阈值比较，保存满足条件的node_diff，同时将相同的node_diff值进行计数；
(4)获取出现次数最多的node_diff中的xdiff与ydiff；
(5)定义与SAR图像尺寸相关参数imgSarx与imgSary，与卫星数字地图尺寸相关参数imgSatex与imgSatey；
Figure BDA0003592648550000121
Figure BDA0003592648550000122
(6)计算变换矩阵
Figure BDA0003592648550000123
则有：/>
Figure BDA0003592648550000124
(7)计算平移参数
Figure BDA0003592648550000131
S45：结合尺度参数s与变换矩阵transMatric计算出SAR图像在卫星数字地图中的匹配点MP(x，y)，设SAR图像中的目标点为C(x，y)，则可以得到匹配点MP(x，y)：
Figure BDA0003592648550000132
匹配结果如图6所不，其中A点为SAR图像中的目标点，B点为卫星数字地图中与A对应的匹配点。
S5：采用匹配点坐标计算目标点经纬度信息
采用载入的经纬度数据库，结合匹配点坐标MP(x，y)计算出目标点经纬度信息，如图7所示。
如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (5)

1.基于LSD直线检测的SAR图像配准方法，其特征在于，包括如下步骤：
S1：读取SAR图像与卫星数字地图，构建待配准图像对
Figure DEST_PATH_IMAGE001
、
Figure 675103DEST_PATH_IMAGE002
；
S2：对待配准图像对
Figure 890184DEST_PATH_IMAGE001
、
Figure 651466DEST_PATH_IMAGE002
进行预处理，得到预处理后的图像
Figure DEST_PATH_IMAGE003
、
Figure 711826DEST_PATH_IMAGE004
；
S3：采用LSD直线检测方法检测S2预处理后的图像
Figure 140534DEST_PATH_IMAGE003
、
Figure 260936DEST_PATH_IMAGE004
中的直线特征，保存检测到的直线信息，以点数据的形式进行存储，并表示为
Figure DEST_PATH_IMAGE005
、
Figure 243936DEST_PATH_IMAGE006
；
S4：根据
Figure 373566DEST_PATH_IMAGE005
、
Figure 391200DEST_PATH_IMAGE006
计算出角度参数
Figure DEST_PATH_IMAGE007
、尺度参数
Figure 682505DEST_PATH_IMAGE008
、变换矩阵
Figure DEST_PATH_IMAGE009
和平移参数，然后对图像
Figure 884291DEST_PATH_IMAGE001
进行尺度变换、平移变换和变换矩阵处理，得到处理后图像中的目标点在卫星数字地图
Figure 552033DEST_PATH_IMAGE002
中对应的匹配点，表示为
Figure 424174DEST_PATH_IMAGE010
；
所述S4包括：
S41：将直线特征的端点坐标
Figure 417538DEST_PATH_IMAGE005
、
Figure 640709DEST_PATH_IMAGE006
转化为极坐标，表示为
Figure DEST_PATH_IMAGE011
、
Figure 315404DEST_PATH_IMAGE012
；
S411、计算每条直线段的斜率
Figure DEST_PATH_IMAGE013
与截距
Figure 42051DEST_PATH_IMAGE014
；
S412、计算直线段的中心点坐标
Figure DEST_PATH_IMAGE015
、
Figure 675158DEST_PATH_IMAGE016
及
Figure DEST_PATH_IMAGE017
Figure 588887DEST_PATH_IMAGE018
(8)
Figure DEST_PATH_IMAGE019
(9)
Figure 332852DEST_PATH_IMAGE020
(10)
S413、计算极坐标的极角与极径：
当
Figure DEST_PATH_IMAGE021
时，获得如下三个点：
Figure DEST_PATH_IMAGE023
Figure DEST_PATH_IMAGE025
Figure DEST_PATH_IMAGE027
当
Figure 600759DEST_PATH_IMAGE028
时，获得如下三个点：
Figure 404767DEST_PATH_IMAGE023
Figure 868109DEST_PATH_IMAGE030
Figure 415765DEST_PATH_IMAGE027
得到极角
Figure DEST_PATH_IMAGE031
：
Figure DEST_PATH_IMAGE033
Figure DEST_PATH_IMAGE035
Figure DEST_PATH_IMAGE037
Figure DEST_PATH_IMAGE039
得到极径
Figure 992371DEST_PATH_IMAGE040
：
令
Figure DEST_PATH_IMAGE041
，
则有
Figure 701701DEST_PATH_IMAGE042
；
S42：以卫星数字地图为基准，计算SAR图像中直线与卫星数字地图中直线的角度差
Figure 386760DEST_PATH_IMAGE007
，
计算
Figure DEST_PATH_IMAGE043
与
Figure 472528DEST_PATH_IMAGE044
中
Figure DEST_PATH_IMAGE045
的差值：
Figure 28274DEST_PATH_IMAGE046
，
将相同的
Figure DEST_PATH_IMAGE047
值进行计数，出现次数最多的
Figure 905576DEST_PATH_IMAGE047
值记为
Figure 77931DEST_PATH_IMAGE007
；
S43：计算尺度参数
Figure 967390DEST_PATH_IMAGE048
，计算步骤如下：
S431、采用
Figure DEST_PATH_IMAGE049
中的
Figure 643222DEST_PATH_IMAGE031
统计出SAR图像中平行直线组；
S432、计算SAR图像与卫星数字地图中每组平行直线极坐标的极径差
Figure DEST_PATH_IMAGE051
Figure DEST_PATH_IMAGE053
Figure DEST_PATH_IMAGE055
将相同的
Figure 632038DEST_PATH_IMAGE056
值进行计数，将出现次数最多的
Figure 26110DEST_PATH_IMAGE056
值记为尺度参数
Figure 719259DEST_PATH_IMAGE008
；
S44、计算变换矩阵
Figure 780756DEST_PATH_IMAGE009
与平移参数
S441、分别计算SAR图像与卫星数字地图中直线特征的交点，分别表示为
Figure DEST_PATH_IMAGE057
、
Figure 2790DEST_PATH_IMAGE058
，
S442、分别计算交点集合
Figure 884158DEST_PATH_IMAGE057
与
Figure 380999DEST_PATH_IMAGE058
中各点间坐标的差分
Figure DEST_PATH_IMAGE059
，
S443、将
Figure 765844DEST_PATH_IMAGE060
中的横坐标差值与纵坐标差值分别与阈值比较，保存满足条件的
Figure 692867DEST_PATH_IMAGE060
，同时将相同的
Figure 61532DEST_PATH_IMAGE060
值进行计数；
S444、获取出现次数最多的
Figure 96484DEST_PATH_IMAGE060
中的
Figure DEST_PATH_IMAGE061
、
Figure 335835DEST_PATH_IMAGE062
；
S445、定义与SAR图像尺寸相关参数
Figure DEST_PATH_IMAGE063
、
Figure 899672DEST_PATH_IMAGE064
，与卫星数字地图尺寸相关参数
Figure DEST_PATH_IMAGE065
、
Figure 490053DEST_PATH_IMAGE066
Figure 63117DEST_PATH_IMAGE068
Figure 688133DEST_PATH_IMAGE070
S446、计算变换矩阵
Figure DEST_PATH_IMAGE071
，则有：
Figure DEST_PATH_IMAGE073
S447、计算平移参数
Figure 157292DEST_PATH_IMAGE074
；
S45、结合尺度参数
Figure 234969DEST_PATH_IMAGE008
、变换矩阵
Figure 142882DEST_PATH_IMAGE009
、平移参数计算出处理后SAR图像在卫星数字地图中的匹配点
Figure 622405DEST_PATH_IMAGE010
，设处理后SAR图像中的目标点为
Figure DEST_PATH_IMAGE075
，则得到匹配点
Figure 525115DEST_PATH_IMAGE010
：
Figure 355667DEST_PATH_IMAGE076
；
S5：载入卫星数字地图
Figure 270534DEST_PATH_IMAGE002
配套的经纬度数据库，根据匹配点坐标
Figure 870142DEST_PATH_IMAGE010
获取目标点经纬度信息，表示为
Figure 743420DEST_PATH_IMAGE078
、
Figure DEST_PATH_IMAGE079
。
2.根据权利要求1所述的基于LSD直线检测的SAR图像配准方法，其特征在于：所述S1中SAR图像与卫星数字地图以多维矩阵的格式进行存储。
3.根据权利要求1所述的基于LSD直线检测的SAR图像配准方法，其特征在于：所述S2中对
Figure 264532DEST_PATH_IMAGE002
进行的预处理操作包括灰度化、滤波、二值化、形态学运算、高斯滤波、Canny边缘检测，对
Figure 248668DEST_PATH_IMAGE001
进行的预处理操作包括灰度化、滤波、二值化、形态学运算、相干斑抑制、高斯滤波、Canny边缘检测。
4.根据权利要求3所述的基于LSD直线检测的SAR图像配准方法，其特征在于，所述预处理操作过程如下：
S21：将待配准图像对进行灰度化处理
设
Figure 437204DEST_PATH_IMAGE080
、
Figure DEST_PATH_IMAGE081
、
Figure 950225DEST_PATH_IMAGE082
分别表示待处理图像的三个通道，采用平均值法进行灰度化处理，将三通道转换为单通道，得到图中
Figure DEST_PATH_IMAGE083
处的像素灰度值为：
Figure 958632DEST_PATH_IMAGE084
(1)
Figure DEST_PATH_IMAGE085
、
Figure 215301DEST_PATH_IMAGE086
、
Figure DEST_PATH_IMAGE087
表示三个通道
Figure 258344DEST_PATH_IMAGE080
、
Figure 207845DEST_PATH_IMAGE081
、
Figure 234707DEST_PATH_IMAGE082
在图中
Figure 551857DEST_PATH_IMAGE083
处的像素灰度值；
将待配准图像对
Figure 980564DEST_PATH_IMAGE001
、
Figure 100967DEST_PATH_IMAGE002
按公式(1)处理，得到SAR图像与卫星数字地图的灰度图；
S22：对灰度化的待配准图像对进行滤波处理，然后对滤波结果进行二值化处理，再根据二值化图像中非零像素个数的占比重复滤波及二值化操作，迭代至非零像素个数占比达到预设阈值，包括：
对S21中灰度图进行归一化处理，使像素灰度值在0~1范围内，设
Figure 615125DEST_PATH_IMAGE088
为图像中最大灰度值，
Figure DEST_PATH_IMAGE089
为图像中最小灰度值，通过下式进行图像归一化：
Figure 948017DEST_PATH_IMAGE090
(2)
设置卷积核kernel的大小，以核内中心位置对应的像素为目标像素，使用kernel对待配准图像对进行卷积操作，卷积结果为：
Figure DEST_PATH_IMAGE091
按从左到右、从上到下滑动卷积核，直到处理完图像最后一个像素，
对卷积结果进行二值化处理，得到二值化后的结果
Figure DEST_PATH_IMAGE093
、
Figure 903335DEST_PATH_IMAGE094
，凸显待配准图像中关键区域的轮廓边缘信息；
S23：对
Figure DEST_PATH_IMAGE095
、
Figure 194639DEST_PATH_IMAGE096
进行形态学处理，将断裂区域连接起来，得到处理结果
Figure DEST_PATH_IMAGE097
、
Figure 399355DEST_PATH_IMAGE098
；
S24：根据SAR图像
Figure 67097DEST_PATH_IMAGE001
的距离向分辨率Dre与方位向分辨率Dis构建滤波核
Figure DEST_PATH_IMAGE099
，大小为
Figure 939238DEST_PATH_IMAGE100
，以核中心点为原点，通过下式构建：
Figure DEST_PATH_IMAGE101
相干斑抑制步骤如下：
S241、采用构建的滤波核在二值SAR图像
Figure 135864DEST_PATH_IMAGE097
区域进行滑动操作，从左到右、从上到下，
S242、判断滤波核中元素与对应二值SAR图像
Figure 90526DEST_PATH_IMAGE097
中元素是否相等且等于255，记录符合条件的元素个数，记为
Figure 561959DEST_PATH_IMAGE102
，
S243、判断滤波核中元素值为255，且对应二值SAR图像
Figure 819765DEST_PATH_IMAGE097
中元素值为0，记录符合条件的元素个数，记为
Figure DEST_PATH_IMAGE103
，
S244、当
Figure 921713DEST_PATH_IMAGE102
与
Figure 632180DEST_PATH_IMAGE103
的值大于滤波核元素个数的1/3时，判断图像中每行左右邻域或每列上下邻域中像素灰度值为255的像素个数是否超过5，如果超过则设置目标像素的灰度值为0，
S245、重复S241~S244，直到滤波核中心位置与图像右下角像素重合，得到处理结果
Figure 907304DEST_PATH_IMAGE104
；
S25：对卫星数字地图
Figure 754037DEST_PATH_IMAGE098
和相干斑抑制后的SAR图像
Figure DEST_PATH_IMAGE105
进行高斯滤波去噪、Canny边缘检测，得到处理后的图像
Figure 558045DEST_PATH_IMAGE003
、
Figure 490229DEST_PATH_IMAGE004
。
5.根据权利要求4所述的基于LSD直线检测的SAR图像配准方法，其特征在于，所述S3采用 LSD直线检测方法检测图像
Figure 569043DEST_PATH_IMAGE003
与
Figure 535862DEST_PATH_IMAGE004
中的直线特征，包括：
S31：首先计算图像中所有像素点的梯度大小，记为
Figure 510772DEST_PATH_IMAGE106
，根据梯度大小
Figure 195831DEST_PATH_IMAGE106
计算出水平角度
Figure DEST_PATH_IMAGE107
与梯度幅值
Figure 16019DEST_PATH_IMAGE108
，如下式所示：
Figure DEST_PATH_IMAGE109
(5)
Figure 309116DEST_PATH_IMAGE110
(6)
Figure DEST_PATH_IMAGE111
(7)
其中，
Figure 189347DEST_PATH_IMAGE112
为图像
Figure DEST_PATH_IMAGE113
在第
Figure 564965DEST_PATH_IMAGE114
行第
Figure DEST_PATH_IMAGE115
个像素上的像素值，
利用伪排序得到梯度幅值大的点作为种子点，以该点水平线角度
Figure 454424DEST_PATH_IMAGE107
作为区域的初始角度
Figure 130256DEST_PATH_IMAGE116
，在八邻域中寻找与
Figure 978126DEST_PATH_IMAGE116
的偏差小于预设阈值
Figure DEST_PATH_IMAGE117
的点，将该点加入到区域中并更新
Figure 841040DEST_PATH_IMAGE116
，得到线支持区域，最后通过矩形逼近得到直线特征，检测到的直线特征以点数据的形式进行存储，分别表示为
Figure 799768DEST_PATH_IMAGE118
、
Figure DEST_PATH_IMAGE119
；
S32：筛选
Figure 64528DEST_PATH_IMAGE118
、
Figure 817720DEST_PATH_IMAGE119
中最具代表性的直线特征来进行配准工作，筛选步骤如下：
S321、计算
Figure 699088DEST_PATH_IMAGE118
、
Figure 195929DEST_PATH_IMAGE119
中所有直线段的长度，并按降序排列，表示为
Figure 111932DEST_PATH_IMAGE120
、
Figure DEST_PATH_IMAGE121
，
S322、根据预设直线特征数量
Figure 501938DEST_PATH_IMAGE122
，在
Figure 136181DEST_PATH_IMAGE120
、
Figure 171134DEST_PATH_IMAGE121
中提取前
Figure 676064DEST_PATH_IMAGE122
条直线特征，
S323、提取距离间隔较大的直线特征，
S324、保存直线特征的端点坐标，表示为
Figure 36638DEST_PATH_IMAGE005
、
Figure 627020DEST_PATH_IMAGE006
。
CN202210398024.1A 2022-04-12 2022-04-12 基于lsd直线检测的sar图像配准方法 Active CN114723794B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210398024.1A CN114723794B (zh) 2022-04-12 2022-04-12 基于lsd直线检测的sar图像配准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210398024.1A CN114723794B (zh) 2022-04-12 2022-04-12 基于lsd直线检测的sar图像配准方法

Publications (2)

Publication Number Publication Date
CN114723794A CN114723794A (zh) 2022-07-08
CN114723794B true CN114723794B (zh) 2023-03-24

Family

ID=82243152

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210398024.1A Active CN114723794B (zh) 2022-04-12 2022-04-12 基于lsd直线检测的sar图像配准方法

Country Status (1)

Country Link
CN (1) CN114723794B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115457222B (zh) * 2022-09-14 2023-06-16 北京建筑大学 一种在地理信息系统中地理配准三维模型的方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103514606A (zh) * 2013-10-14 2014-01-15 武汉大学 一种异源遥感影像配准方法
CN107292922A (zh) * 2017-06-23 2017-10-24 电子科技大学 一种用于光学与合成孔径雷达图像配准的方法
CN113077514A (zh) * 2021-06-03 2021-07-06 浙江大学 一种sar图像船只弱尾迹增强与检测方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102663725B (zh) * 2012-03-05 2014-07-16 西北工业大学 基于线特征和控制点的可见光和sar图像配准方法
CN103345757B (zh) * 2013-07-19 2015-12-02 武汉大学 多层次多特征约束下的光学和sar影像自动配准方法
CN106886977B (zh) * 2017-02-08 2021-02-05 徐州工程学院 一种多图自动配准及融合拼接方法
CN109308715A (zh) * 2018-09-19 2019-02-05 电子科技大学 一种基于点特征和线特征结合的光学图像配准方法
CN111145228B (zh) * 2019-12-23 2023-05-26 西安电子科技大学 基于局部轮廓点与形状特征融合的异源图像配准方法
CN113870324A (zh) * 2020-06-29 2021-12-31 上海微创卜算子医疗科技有限公司 多模态图像的配准方法及其配准装置和计算机可读存储介质

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103514606A (zh) * 2013-10-14 2014-01-15 武汉大学 一种异源遥感影像配准方法
CN107292922A (zh) * 2017-06-23 2017-10-24 电子科技大学 一种用于光学与合成孔径雷达图像配准的方法
CN113077514A (zh) * 2021-06-03 2021-07-06 浙江大学 一种sar图像船只弱尾迹增强与检测方法

Also Published As

Publication number Publication date
CN114723794A (zh) 2022-07-08

Similar Documents

Publication Publication Date Title
US11403839B2 (en) Commodity detection terminal, commodity detection method, system, computer device, and computer readable medium
CN111626190B (zh) 基于聚类分区进行刻度识别的水位监测方法
CN111222474B (zh) 一种任意尺度的高分辨率图像小目标检测方法
Rao et al. Textural analysis of IRS-1D panchromatic data for land cover classification
Palenichka et al. Automatic extraction of control points for the registration of optical satellite and LiDAR images
CN112085772B (zh) 一种遥感图像配准方法及装置
Choi et al. Vehicle detection from aerial images using local shape information
Wang et al. A robust multisource image automatic registration system based on the SIFT descriptor
Mukherjee et al. Interest points for hyperspectral image data
CN108564532B (zh) 大尺度地距星载sar图像镶嵌方法
CN107862319B (zh) 一种基于邻域投票的异源高分光学影像匹配误差剔除方法
CN112734729B (zh) 适用于夜间补光条件的水尺水位线图像检测方法、装置及存储介质
CN114723794B (zh) 基于lsd直线检测的sar图像配准方法
CN110569861A (zh) 一种基于点特征和轮廓特征融合的图像匹配定位方法
CN111144250A (zh) 融合雷达和光学遥感数据的土地覆盖分类方法
Xing et al. Integrating change magnitude maps of spectrally enhanced multi-features for land cover change detection
CN114897705A (zh) 一种基于特征优化的无人机遥感图像拼接方法
CN110246165B (zh) 提高可见光图像与sar图像配准速度的方法及系统
CN113793309B (zh) 一种基于形态学特征的亚像素级椭圆检测方法
CN114972453B (zh) 改进的基于lsd与模板匹配的sar图像区域配准方法
Xia et al. A table method for coded target decoding with application to 3-D reconstruction of soil specimens during triaxial testing
CN107832793A (zh) 一种高光谱图像的分类方法及系统
WO2024040937A1 (zh) 体积测量容器读数方法、装置、存储介质和计算机设备
CN116206139A (zh) 一种基于局部自卷积的无人机图像升尺度匹配方法
CN115049708B (zh) 基于lsd直线检测与模板匹配的sar图像配准方法

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
CP01 Change in the name or title of a patent holder
CP01 Change in the name or title of a patent holder

Address after: 210000 Room 301, floor 3, building 75, zone B, entrepreneurship and innovation city, No. 15, Fengji Avenue, Yuhuatai District, Nanjing, Jiangsu Province

Patentee after: Nanjing Thunderbolt Information Technology Co.,Ltd.

Address before: 210000 Room 301, floor 3, building 75, zone B, entrepreneurship and innovation city, No. 15, Fengji Avenue, Yuhuatai District, Nanjing, Jiangsu Province

Patentee before: NANJING LEADING INFORMATION TECHNOLOGY Co.,Ltd.