CN113096163A - 一种无需先验升降轨信息的星载sar影像高精度匹配方法 - Google Patents

一种无需先验升降轨信息的星载sar影像高精度匹配方法 Download PDF

Info

Publication number
CN113096163A
CN113096163A CN202110408520.6A CN202110408520A CN113096163A CN 113096163 A CN113096163 A CN 113096163A CN 202110408520 A CN202110408520 A CN 202110408520A CN 113096163 A CN113096163 A CN 113096163A
Authority
CN
China
Prior art keywords
image
scene
borne sar
satellite
image plane
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.)
Granted
Application number
CN202110408520.6A
Other languages
English (en)
Other versions
CN113096163B (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.)
Xiangtan University
Original Assignee
Xiangtan University
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 Xiangtan University filed Critical Xiangtan University
Priority to CN202110408520.6A priority Critical patent/CN113096163B/zh
Publication of CN113096163A publication Critical patent/CN113096163A/zh
Application granted granted Critical
Publication of CN113096163B publication Critical patent/CN113096163B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/60Rotation of whole images or parts thereof
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/60Analysis of geometric attributes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/70Determining position or orientation of objects or cameras
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • G06T2207/10044Radar image

Landscapes

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

Abstract

本发明涉及一种无需先验升降轨信息的星载SAR影像高精度匹配方法,包括:计算一景星载SAR影像的四个顶点的地理坐标;根据步骤a所得地理坐标计算其在另一景星载SAR影像上的像平面坐标;根据步骤b所得像平面坐标形成的四边形,选取其最大内接矩形,作为两景星载SAR影像的重叠区域;判断两景星载SAR影像成像轨道与观测方向之间的关系,并将两景星载SAR影像翻转为同轨同侧;根据SAR‑SIFT进行两景星载SAR影像的匹配。本发明提供一种在未知星载SAR影像轨道和观测方向时,得到有重叠区域的两景星载SAR影像的同名点,整个算法具有准确度高、鲁棒性好的优点。

Description

一种无需先验升降轨信息的星载SAR影像高精度匹配方法
技术领域
本发明涉及通信领域,特别涉及一种无需先验升降轨信息的星载SAR影像高精度匹配方法。
背景技术
合成孔径雷达(SAR)技术的不断进步,使人们从SAR影像中获取更多地表信息成为可能。影像匹配是指通过一定的匹配算法在两幅或多幅影像之间识别同名点的过程。它是图像融合、目标识别、目标变化检测、计算机视觉等问题中的一个重要前期步骤,在遥感、数字摄影测量、计算机视觉、地图学及军事应用等多个领域都有着广泛的应用,对卫星影像的配准、拼接、3D重建等具有重要意义。
绝大多数针对星载影像的匹配都是同轨同侧SAR影像之间的匹配,针对如异轨异侧等其它情况的星载SAR影像匹配方法多是先将其中一景影像通过正射纠正的方式将其纠正到和另一景影像同一视角后再进行匹配,进行正射纠正后进行匹配,匹配到的同名点会增加正射纠正这一步所带来的累积误差。另外,星载影像较大,直接匹配结果较差,耗费的计算资源也较大。
发明内容
针对现有技术存在的上述问题,本发明提供一种通过先行判断两景SAR影像之间的成像轨道和观测方位之间的关系,然后直接进行匹配的方法,解决了非同轨同侧的两景影像需要正射纠正后匹配的问题,并通过分块匹配的方案对两景影像进行匹配。
为实现上述目的,本发明采取以下技术方案:一种无需先验升降轨信息的星载SAR影像高精度匹配方法,其具体包括以下步骤:
步骤a)计算一景星载SAR影像的四个顶点的地理坐标;
步骤b)根据步骤a)所得地理坐标计算其在另一景星载SAR影像上的像平面坐标;
步骤c)根据步骤b)所得像平面坐标形成的四边形,选取其最大内接矩形,作为两景星载SAR影像的重叠区域;
步骤d)判断两景星载SAR影像成像轨道与观测方向之间的关系,并将步骤c)所得两景星载SAR影像的重叠区域翻转为同轨同侧;
步骤e)进行两景星载SAR影像的匹配。
优选地,步骤a)的具体操作为:根据一景星载SAR影像的顶点坐标,设置初始高程值H0=100,代入公式(1)和公式(2)计算一景星载SAR影像该顶点对应的经纬度坐标,再根据全球高程模型获取该经纬度位置的实际高程值H1判断△H=|H1-H0|的大小,当△H≥0.001时,再将高程值H1代入公式(1)和公式(2)进行迭代,直至△H=|Hn-Hn-1|<0.001时求得其所对应的P、L后停止迭代,即得一景星载SAR影像的顶点坐标对应的地理坐标(P,L,Hn)。
更优选地,公式(1)和公式(2)具体为:
Figure BDA0003023245740000021
Figure BDA0003023245740000022
其中,X表示像平面的横坐标,Y表示像平面的纵坐标,H表示像平面坐标所在实际地理空间的高程,P表示像平面坐标所在实际地理空间的经度,L表示像平面坐标所在实际地理空间的纬度;且,
NumP(X,Y,H)=a1+a2Y+a3X+a4H+a5YX+a6YH+a7XH+a8Y2+a9X2+a10H2+a11XYH+a12Y3+a13YX2+a14YH2+a15Y2X+a16X3+a17XH2+a18Y2H+a19X2H+a20H3
DenP(X,Y,H)=b1+b2Y+b3X+b4H+b5YX+b6YH+b7XH+b8Y2+b9X2+b10H2+b11XYH+b12Y3+b13YX2+b14YH2+b15Y2X+b16X3+b17XH2+b18Y2H+b19X2H+b20H3
NumL(X,Y,H)=c1+c2Y+c3X+c4H+c5YX+c6YH+c7XH+c8Y2+c9X2+c10H2+c11XYH+c12Y3+c13YX2+c14YH2+c15Y2X+c16X3+c17XH2+c18Y2H+c19X2H+c20H3
DenL(X,Y,H)=d1+d2Y+d3X+d4H+d5YX+d6YH+d7XH+d8Y2+d9X2+d10H2+d11XYH+d12Y3+d13YX2+d14YH2+d15Y2X+d16X3+d17XH2+d18Y2H+d19X2H+d20H3
其中,ai,bi,ci,di;i=1,2,…,20为星载SAR影像的RPC参数。
优选地,步骤b)的具体操作为:根据步骤a)所得地理坐标,和另一景星载SAR影像的RPC参数,代入公式(3)和公式(4)求解一景星载SAR影像的顶点对应于另一景星载SAR影像上的像平面坐标。
更优选地,公式(3)和公式(4)具体为:
Figure BDA0003023245740000031
Figure BDA0003023245740000032
其中,
NumX(P,L,H)=a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2+a18L2H+a19P2H+a20H3
DemX(P,L,H)=b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2+b18L2H+b19P2H+b20H3
NumY(P,L,H)=c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2+c18L2H+c19P2H+c20H3
DemY(P,L,H)=d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2+。d18L2H+d19P2H+d20H3
优选地,步骤c)计算最大内接矩形的具体操作为:分别对解算的横纵坐标按从小到大的顺序排序,得从小到大的横坐标顺序和纵坐标顺序,再根据所得横坐标顺序和纵坐标顺序组成其最大内接矩形。
优选地,步骤d)的判断方法为:根据步骤a和步骤b计算一景星载SAR影像的四个顶点在其“同轨同侧”、“同轨异侧”、“异轨同侧”和“异轨异侧”的星载SAR影像的像平面坐标,对一景星载SAR影像的四个顶点进行标号并将其对应的像平面坐标同时对应进行标号,得到一景星载SAR影像对应其不同成像轨道和观测方向之间的对应关系,再将待测两景星载SAR影像的四个顶点和其对应像平面坐标的标号与所得对应关系进行比较,得两景星载SAR影像的实际成像轨道和观测方向之间的对应关系。
优选地,所述方法还可包括步骤d’)将两景星载SAR影像分割成若干不重叠的小块影像对,以便节省计算过程中所需计算机的算力;所述步骤d’)可添加在步骤d)前或后。
优选地,步骤e)包括如下步骤:
步骤e1)将已翻转的两景星载SAR影像(或分割后的若干组小块影像对)进行关键点检测,计算其特征描述子;
步骤e2)根据步骤e1)计算所得特征描述子,对检测出的关键点进行暴力匹配;
步骤e3)剔除步骤e2)匹配结果中的误匹配点。
更优选地,步骤e1)具体为:检测步骤d’)所得一两景小块影像对或步骤d)所得两景星载SAR影像的关键点,计算特征描述子。
进一步,步骤e1)中检测两景星载SAR影像或两景小块影像对的关键点的具体实现方式如下:
⑤计算影像水平方向和垂直方向的比值梯度:
Figure BDA0003023245740000051
Figure BDA0003023245740000052
⑥计算协方差矩阵:
Figure BDA0003023245740000053
⑦计算角点响应值:
RSH(x,y,α)=det(GSH(x,y,α))-d·(trace(GSH(x,y,α)))2 (8)
其中d的取值范围为0.04到0.08;
⑧通过角点响应局部极大值抑制,计算求得关键点坐标。
进一步,步骤e1)中计算关键点特征描述子的具体实现方式如下:
③计算影像每个像素点的梯度与幅值:
Figure BDA0003023245740000054
Figure BDA0003023245740000055
④以关键点为中心,统计大小为16*16的邻域范围内的梯度方向直方图,将16*16大小的区域分为16块4*4大小的小区域,统计每个小区域8个梯度方向的幅值直方图,合并所有小区域形成128维的特征向量即为该关键点的特征描述子。
更优选地,步骤e2)中,通过欧式距离(公式(11))的相似性度量,采用最近邻原则进行暴力匹配,得到初始匹配结果:
Figure BDA0003023245740000061
其中ml表示一景影像分割后的小块影像的特征描述子向量,mr表示另一景影像分割后的小块影像的特征描述子向量。
更优选地,步骤e3)中剔除步骤e2)匹配结果中的误匹配点的具体实现方式如下:
①随机从匹配结果的同名点集中随机抽出4个样本数据,计算出变换矩阵:
Figure BDA0003023245740000062
记为模型M,
②计算同名点集中所有数据与模型M的投影误差,若误差小于阈值,加入内点集I,
③如果当前内点集I元素个数大于最优内点集I_best,则更新I_best=I,同时更新迭代次数k,
④如果迭代次数大于k,则退出;否则迭代次数加1,并重复上述步骤。
进一步,若步骤e3)是匹配若干小块影像对,需合并匹配得到的所有同名点,在全局情况下再次剔除误匹配点,剔除误匹配点的具体实现方式如下:
①随机从匹配结果的同名点集中随机抽出4个样本数据,计算出变换矩阵:
Figure BDA0003023245740000063
记为模型M,
②计算同名点集中所有数据与模型M的投影误差,若误差小于阈值,加入内点集I,
③如果当前内点集I元素个数大于最优内点集I_best,则更新I_best=I,同时更新迭代次数K,
④如果迭代次数大于K=1000,则退出;否则迭代次数加1,并重复上述步骤。
与现有技术相比,本发明的优势如下:
1、本发明无需知晓两景星载SAR影像间的成像轨道与观测方向间的关系,通过星载影像附带的RPC参数信息计算可以得到两景影像间的关系。
2、采用分块匹配的策略,节省计算过程中所需计算机的算力。
3、通过自动判断两景影像间的关系,省略非同轨同侧影像匹配过程中的正射纠正步骤,省去了人为干预过程,大大节省人力成本。
附图说明
图1为本发明提供的无需先验升降轨信息的星载SAR影像高精度匹配方法的流程图;
图2为本发明提供的用于匹配的两景异轨同侧SAR影像;
图3为本发明提供的用于匹配的两景异轨同侧SAR影像的重叠区域;
图4为本发明提供的两景星载影像成像轨道和观测方向判定示例图;
图5为本发明提供的用于判断两景SAR影像为异轨同侧关系的实例图;
图6为本发明提供的用于匹配的两景异轨同侧SAR影像划分后的其中一对小块影像对;
图7为本发明提供的用于匹配的两景异轨同侧SAR影像划分后的其中一对小块影像对剔除误匹配点后的实际匹配效果图;
图8为本发明提供的用于匹配的两景异轨同侧SAR影像划分后的另一对小块影像对剔除误匹配点后的实际匹配效果图;
图9为本发明提供的用于匹配的两景异轨同侧SAR影像合并所有小块影像对同名点的匹配效果图;
图10为本发明提供的用于匹配的两景异轨同侧SAR影像最终剔除误匹配点后同名点的匹配效果图。
具体实施方式
下面结合附图和实施例对本发明做进一步具体的说明。
如图1所示,本发明提供的无需先验升降轨信息的星载SAR影像高精度匹配方法,具体包括以下步骤:
1)计算两景星载SAR影像重叠区域:
由于星载SAR影像较大,直接匹配比较耗时且对计算机硬件要求较高,故确定两景星载SAR影像的重叠区域有助于减小实际的匹配区域,如图2所示,为本实施例用于匹配的两景异轨同侧星载SAR影像。
现有技术中,常规的匹配过程需要搜索匹配区域,本发明省去了搜索匹配区域这一操作,通过公式(1)~公式(4)计算直接得到星载SAR影像中明确的匹配范围,即通过两景星载SAR影像包含的有理多项式系数参数(下称RPC参数)及全球高程模型图像文件,求解第一景星载SAR影像已知的四个顶点(0,0)、(wL,0)、(wL,hL)、(0,hL)在第二景星载SAR影像中对应的像平面坐标(X1′,Y1′)、(X2′,Y2′)、(X3′,Y3′)、(X4′,Y4′),其中hL、wL分别表示第一景星载SAR影像的高与宽,将第二景星载SAR影像上述像平面坐标所组成的四边形的最大内接矩形作为两景星载SAR影像的重叠区域。
具体步骤如下:
1.1向软件中输入第一景星载SAR影像的顶点坐标(0,0),设置初始高程值H0=100,代入公式(1)和公式(2)计算第一景星载SAR影像顶点对应的经纬度坐标,再根据全球高程模型中的数据获取所在经纬度位置的实际高程值H1,判断△H=|H1-H0|的大小,当△H≥0.001时,将新的高程值H1代入公式(1)和(2),求得新的经纬度坐标,再在全球高程模型中获取新的经纬度坐标对应的实际高程值H2,判断△H=|H2-H1|的大小,如此往复,直到ΔH<0.001,即得第一景星载SAR影像的顶点所对应的准确高程值Hn。
Figure BDA0003023245740000081
Figure BDA0003023245740000091
其中,X表示像平面的横坐标,Y表示像平面的纵坐标,H表示像平面坐标所在实际地理空间的高程,P表示像平面坐标所在实际地理空间的经度,L表示像平面坐标所在实际地理空间的纬度;且,
NumP(X,Y,H)=a1+a2Y+a3X+a4H+a5YX+a6YH+a7XH+a8Y2+a9X2+a10H2+a11XYH+a12Y3+a13YX2+a14YH2+a15Y2X+a16X3+a17XH2+a18Y2H+a19X2H+a20H3
DenP(X,Y,H)=b1+b2Y+b3X+b4H+b5YX+b6YH+b7XH+b8Y2+b9X2+b10H2+b11XYH+b12Y3+b13YX2+b14YH2+b15Y2X+b16X3+b17XH2+b18Y2H+b19X2H+b20H3
NumL(X,Y,H)=c1+c2Y+c3X+c4H+c5YX+c6YH+c7XH+c8Y2+c9X2+c10H2+c11XYH+c12Y3+c13YX2+c14YH2+c15Y2X+c16X3+c17XH2+c18Y2H+c19X2H+c20H3
DenL(X,Y,H)=d1+d2Y+d3X+d4H+d5YX+d6YH+d7XH+d8Y2+d9X2+d10H2+d11XYH+d12Y3+d13YX2+d14YH2+d15Y2X+d16X3+d17XH2+d18Y2H+d19X2H+d20H3
其中,ai,bi,ci,di;i=1,2,…,20为星载SAR影像的RPC参数,此处为随高分3号卫星发来的第一星载SAR影像所附带的RPC参数。
将顶点坐标(0,0)与求解得到的准确高程值Hn代入公式(1)和公式(2),即可得到顶点坐标(0,0)对应的实际经纬度坐标(P,L),与Hn结合后最终得顶点坐标(0,0)对应的地理坐标(P,L,Hn)。
1.2根据第二景星载SAR影像的RPC参数,和步骤1.1所得第一景星载SAR影像顶点(0,0)对应的地理坐标(P,L,Hn),根据公式(3)和公式(4)求解其对应于第二景星载SAR影像上的像平面坐标(X1′,Y1′):
Figure BDA0003023245740000101
Figure BDA0003023245740000102
其中,
NumX(P,L,H)=a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2+a18L2H+a19P2H+a20H3
DemX(P,L,H)=b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2+b18L2H+b19P2H+b20H3
NumY(P,L,H)=c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2+c18L2H+c19P2H+c20H3
DemY(P,L,H)=d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2+d18L2H+d19P2H+d20H3
其中,ai,bi,ci,di;i=1,2,…,20为星载SAR影像的RPC参数,此处为随高分3号卫星发来的第二星载SAR影像所附带的RPC参数。
1.3重复步骤1.1和1.2分别求得第一景星载SAR影像另外三个顶点(wL,0)、(wL,hL)、(0,hL)在第二景星载SAR影像中的像平面坐标(X2′,Y2′)、(X3′,Y3′)、(X4′,Y4′)。
1.4由于步骤1.1~1.3所得第一景星载SAR影像的四顶点在第二景星载SAR影像中的像平面坐标(X1′,Y1′)、(X2′,Y2′)、(X3′,Y3′)、(X4′,Y4′)所形成的四边形(两景星载SAR影像的实际重叠区域)图形不规则,故选取其最大内接矩形为两景星载SAR影像的重叠区域。具体选择方式为:
分别对解算的横纵坐标按从小到大的顺序排序,如对第一景星载SAR影像在第二景星载SAR影像中的像平面坐标(X1′,Y1′)、(X2′,Y2′)、(X3′,Y3′)、(X4′,Y4′)的横纵坐标分别排序,得到从小到大的横坐标顺序{X2′,X1′,X4′,X3′}、纵坐标顺序{Y3′,Y2′,Y4′,Y1′},则其对应的最大内接矩形的四顶点像平面坐标为(X1′,Y2′)、(X1′,Y4′)、(X4′,Y4′)、(X4′,Y2′)。
如图3中矩形框为选取的两景星载SAR影像通过上述步骤1所述方法得到的实际重叠区域,该重叠区域约为待匹配两景星载SAR影像的一半大小,实际匹配范围缩小,匹配效率明显提高。
类似地,本实施例也可以通过步骤1.1~步骤1.3求解第二景星载SAR影像四顶点(0,0)、(wR,0)、(wR,hR)、(0,hR)在第一景星载SAR影像中的像平面坐标(X1,Y1)、(X2,Y2)、(X3,Y3)、(X4,Y4),其中hR、wR分别表示第二景星载SAR影像的高与宽。再通过步骤1.4的方式得到第二景星载SAR影像在第一景星载SAR影像中最大内接矩形的四顶点坐标,在此不再赘述。
2)判断两景星载SAR影像成像轨道与观测方向之间的关系:
由于现有技术中多采用正射纠正的方式进行对同轨异侧、异轨同侧和异轨异侧的星载SAR影像进行视角纠正,本发明根据预先确定好的两景星载SAR影像成像轨道与观测方向之间的关系,即可将不是同轨同侧的两景星载SAR影像进行翻转操作,人为将两景星载SAR影像处理为“同轨同侧”,以避免正射纠正带来的累积误差。
发明人将一星载SAR影像(a)的四个顶点(0,0)、(w,0)、(w,h)、(0,h)分别编号为①、②、③、④,如图4(a)所示,为某一星载SAR影像的简略图,通过步骤1计算该星载SAR影像(a)的四个顶点在其同轨同侧的星载SAR影像(b)中的像平面坐标,分别将点①、②、③、④对应求得的像平面坐标在星载SAR影像(b)中按照其对应关系也标记为①、②、③、④,如图4(b)所示;类似地,将星载SAR影像(a)的四个顶点按照步骤1计算在其同轨异侧的星载SAR影像(c)、异轨同侧的星载SAR影像(d)和异轨异侧的星载SAR影像(e)中对应的像平面坐标并进行标号,如图4(c)、图4(d)、图4(e)所示。
通过上述方法,可得到一景星载SAR影像的四个顶点对应于另一景星载SAR影像中选取的重叠区域的四个顶点位置的四种可能情况。若如图4(b),即四顶点呈逆时针旋转且①号顶点在左下角,则两景影像为同轨同侧的关系;若如图4(c),即四顶点呈顺时针旋转且①号顶点在右下角,则两景影像为同轨异侧的关系;若如图4(d),即四顶点呈顺时针旋转且①号顶点在左上角,则两景影像为异轨同侧的关系;若如图4(e),即四顶点呈逆时针旋转且①号顶点在右上角,则两景影像为异轨异侧的关系。
判断影像成像轨道与观测方向之间的关系后,若为同轨同侧影像,在分割影像对时不进行操作;若为同轨异侧影像,在分割影像时对第二景影像进行左右翻转操作;若为异轨同侧影像,在分割影像时对第二景影像进行上下翻转操作;若为异轨异侧影像,在分割影像时对第二景影像进行旋转180度操作。
本实施例中,通过步骤1计算第二景星载SAR影像的顶点坐标在第一景星载SAR影像中的像平面坐标并对其标号,如图5所示,通过与图2进行对比可知,本实施例的两景影像为异轨同侧,应对第二景影像进行上下翻转操作,以便后续匹配。
3)将两景星载SAR影像分割成小块影像对:
通过步骤2)已经得到本实施例中两景影像为同轨异侧影像,并对第二景影像进行了上下翻转操作。从两景影像中分割出步骤1)选中的区域,记其重叠区域的大小分别为ML*NL、MR*NR,将第一景星载SAR影像的重叠区域分割成不重叠的若干600*600大小的小块影像,将翻转后的第二景星载SAR影像的重叠区域分割成不重叠的同等数量
Figure BDA0003023245740000121
大小的小块影像,分割后的两景影像的小块影像根据小块影像在重叠区域中的位置一一对应,得到若干分割后的小块影像对,如图6所示。
4)检测步骤3所得小块影像对关键点
关键点检测是图像匹配过程的核心步骤之一,通过角点检测的思路进行关键点检测,具体步骤为:
4.1计算小块影像对的各影像水平方向和垂直方向的比值梯度:
Figure BDA0003023245740000131
Figure BDA0003023245740000132
其中,α为尺度因子,Gx,α为影像在尺度α下水平方向的梯度,Gy,α为影像在尺度α下垂直方向的梯度,log为对数计算符号,M1,α(i)为影像在尺度α下第i个像素点5×5邻域上半部分的均值,M2,α(i)为影像在尺度α下第i个像素点5×5邻域下半部分的均值,M3,α(i)为影像在尺度α下第i个像素点5×5邻域左半部分的均值,M4,α(i)为影像在尺度α下第i个像素点5×5邻域右半部分的均值。
4.2计算协方差矩阵:
Figure BDA0003023245740000133
其中,GSH(x,y,α)为像素点(x,y)在尺度α下的协方差矩阵,
Figure BDA0003023245740000134
为均值为0,标准差为α的高斯核,*为卷积运算符号。
4.3计算角点响应值:
RSH(x,y,α)=det(GSH(x,y,α))-d·(trace(GSH(x,y,α)))2 (8)
其中RSH(x,y,α)为像素点(x,y)在尺度α下的响应值,d为根据现有技术取值,本实施例可以取d=0.04。
4.4求取关键点坐标
RSH(x,y,α)>0.8 (9)
通过选取5×5邻域内的最大角点响应值RSH(x,y,α),并满足式(9)的点的位置,即为求得的关键点坐标。
5)计算特征描述子
根据步骤4)检测出的关键点,采用梯度直方图的思想计算其特征描述子。
具体步骤如下:
5.1计算星载SAR影像的每个像素点的梯度与幅值:
Figure BDA0003023245740000141
Figure BDA0003023245740000142
其中Gt,α为所在像素点的梯度值,Gn,α为所在像素点的梯度方向。
5.2以步骤4.4所得关键点为中心,统计大小为16*16的邻域范围内的梯度方向直方图,将16*16大小的区域分为16块4*4大小的小区域,统计每个小区域8个梯度方向的幅值直方图,合并所有小区域形成128维的特征向量即为该关键点的特征描述子。
6)关键点匹配
根据公式(11)所得欧式距离的相似性度量,采用最近邻原则进行暴力匹配,得初始匹配结果;
Figure BDA0003023245740000143
其中ml为第一景影像分割后的小块影像的特征描述子向量,mr为第二景影像分割后的对应的小块影像的特征描述子向量,Dm为两匹配点描述子之间的欧式距离。
7)误匹配点剔除
初步匹配结果中含有较多误匹配点对,为了得到较为准确的同名点对,通过多次采样迭代的方案,找到最多匹配点子集即为剔除误匹配点后的准确的同名点集。具体步骤如下:
7.1随机从步骤6)所得匹配结果的同名点集中随机抽出4个样本数据(此4个样本之间不能共线),计算出变换矩阵:
Figure BDA0003023245740000151
记为模型M,其中T为求得的透视变换矩阵,t为透视变换矩阵中的具体参数。
7.2计算步骤6)所得匹配结果的同名点集中所有数据与模型M的投影误差,若误差小于阈值,加入内点集I,
7.3如果步骤7.2所得内点集I中元素个数大于最优内点集I_best,则更新I_best=I,同时更新迭代次数K,
7.4如果迭代次数大于K=1000,则退出;否则迭代次数加1,并重复上述步骤7.1~7.3。
图7、图8为采取步骤4)到步骤7)匹配的其中两对小块影像对的匹配结果。
8)对所有小块影像对进行匹配
步骤4)到步骤7)仅是一次匹配的过程,在本发明中,由于采用分块匹配的策略,因此需要重复步骤4)到步骤7),直到得到所有划分的小块影像对的匹配结果,结果图如图9所示(图9中的线条为匹配后同名点的连线)。
9)合并所有同名点对,再次剔除误匹配点
虽然每对小块影像的匹配结果较为准确,但其仅针对于小块影像所在的范围,而合并后的同名点集,其针对实际匹配的星载影像的重叠区域,其中某些同名点的误差也较大,因此需要再次剔除误差较大的匹配点,具体操作步骤同步骤7),得最终匹配所得同名点分布图(图10中的线条为匹配后所有同名点的连线),如图10。
实施例2
本实施例与实施例1的差别仅在于,删除步骤3、8、9,此方法为直接匹配而不需要对两景星载SAR影像进行分割,对计算机要求过高,可按需进行操作。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (10)

1.一种无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,包括以下步骤:
步骤a)计算一景星载SAR影像的四个顶点的地理坐标;
步骤b)根据步骤a)所得地理坐标计算其在另一景星载SAR影像上的像平面坐标;
步骤c)根据步骤b)所得像平面坐标形成的四边形,选取其最大内接矩形,作为两景星载SAR影像的重叠区域;
步骤d)判断两景星载SAR影像成像轨道与观测方向之间的关系,并将步骤c)所得两景星载SAR影像的重叠区域图翻转为同轨同侧;
步骤e)进行两景星载SAR影像的匹配。
2.根据权利要求1所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤a)的具体操作步骤为:根据一景星载SAR影像的顶点坐标,设置初始高程值H0=100,代入公式(1)和公式(2)计算一景星载SAR影像该顶点对应的经纬度坐标,再根据全球高程模型获取该经纬度位置的实际高程值H1判断△H=|H1-H0|的大小,当△H≥0.001时,再将高程值H1代入公式(1)和公式(2)进行迭代,直至△H=|Hn-Hn-1|<0.001时求得其所对应的P、L后停止迭代,即得一景星载SAR影像的顶点坐标对应的地理坐标(P,L,Hn)。
3.根据权利要求2所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤a)中公式(1)具体为:
Figure FDA0003023245730000011
其中,X表示像平面的横坐标,Y表示像平面的纵坐标,H表示像平面坐标所在实际地理空间的高程,P表示像平面坐标所在实际地理空间的经度;且,
NumP(X,Y,H)=a1+a2Y+a3X+a4H+a5YX+a6YH+a7XH+a8Y2+a9X2+a10H2+a11XYH+a12Y3+a13YX2+a14YH2+a15Y2X+a16X3+a17XH2+a18Y2H+a19X2H+a20H3
DenP(X,Y,H)=b1+b2Y+b3X+b4H+b5YX+b6YH+b7XH+b8Y2+b9X2+b10H2+b11XYH+b12Y3+b13YX2+b14YH2+b15Y2X+b16X3+b17XH2+b18Y2H+b19X2H+b20H3
其中,ai,bi;i=1,2,…,20为星载SAR影像的RPC参数。
4.根据权利要求2所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤a)中公式(2)具体为:
Figure FDA0003023245730000021
其中,X表示像平面的横坐标,Y表示像平面的纵坐标,H表示像平面坐标所在实际地理空间的高程,L表示像平面坐标所在实际地理空间的纬度;且,
NumL(X,Y,H)=c1+c2Y+c3X+c4H+c5YX+c6YH+c7XH+c8Y2+c9X2+c10H2+c11XYH+c12Y3+c13YX2+c14YH2+c15Y2X+c16X3+c17XH2+c18Y2H+c19X2H+c20H3
DenL(X,Y,H)=d1+d2Y+d3X+d4H+d5YX+d6YH+d7XH+d8Y2+d9X2+d10H2+d11XYH+d12Y3+d13YX2+d14YH2+d15Y2X+d16X3+d17XH2+d18Y2H+d19X2H+d20H3
其中,ci,di;i=1,2,…,20为星载SAR影像的RPC参数。
5.根据权利要求1所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤b)的具体操作步骤为:根据步骤a)所得地理坐标,和另一景星载SAR影像的RPC参数,代入公式(3)和公式(4)求解一景星载SAR影像的顶点对应于另一景星载SAR影像上的像平面坐标。
6.根据权利要求5所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤b)中公式(3)具体为:
Figure FDA0003023245730000022
其中,X′表示像平面的横坐标,H表示像平面坐标所在实际地理空间的高程,P表示像平面坐标所在实际地理空间的经度,L表示像平面坐标所在实际地理空间的纬度;且,
NumX(P,L,H)=a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2+a18L2H+a19P2H+a20H3
DemX(P,L,H)=b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2+b18L2H+b19P2H+b20H3
其中,ci,di;i=1,2,…,20为星载SAR影像的RPC参数。
7.根据权利要求5所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤b)中公式(4)具体为:
Figure FDA0003023245730000031
其中,Y′表示像平面的纵坐标,H表示像平面坐标所在实际地理空间的高程,P表示像平面坐标所在实际地理空间的经度,L表示像平面坐标所在实际地理空间的纬度;且,
NumY(P,L,H)=c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2+c18L2H+c19P2H+c20H3
DemY(P,L,H)=d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2+d18L2H+d19P2H+d20H3
其中,ci,di;i=1,2,…,20为星载SAR影像的RPC参数。
8.根据权利要求1所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤c)计算最大内接矩形的具体操作为:分别对解算的横纵坐标按从小到大的顺序排序,得从小到大的横坐标顺序和纵坐标顺序,再根据所得横坐标顺序和纵坐标顺序组成其最大内接矩形。
9.根据权利要求1所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,所述方法还可包括步骤d’)将两景星载SAR影像分割成若干不重叠的小块影像对,以便节省计算过程中所需计算机的算力;所述步骤d’)可添加在步骤d)前或后。
10.根据权利要求1所述的无需先验升降轨信息的星载SAR影像高精度匹配方法,其特征在于,步骤e)包括如下步骤:
步骤e1)将已翻转的两景星载SAR影像或分割后的若干组小块影像对进行关键点检测,计算其特征描述子;
步骤e2)根据步骤e1)计算所得特征描述子,对检测出的关键点进行暴力匹配;
步骤e3)剔除步骤e2)匹配结果中的误匹配点。
CN202110408520.6A 2021-04-16 2021-04-16 一种无需先验升降轨信息的星载sar影像高精度匹配方法 Active CN113096163B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110408520.6A CN113096163B (zh) 2021-04-16 2021-04-16 一种无需先验升降轨信息的星载sar影像高精度匹配方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110408520.6A CN113096163B (zh) 2021-04-16 2021-04-16 一种无需先验升降轨信息的星载sar影像高精度匹配方法

Publications (2)

Publication Number Publication Date
CN113096163A true CN113096163A (zh) 2021-07-09
CN113096163B CN113096163B (zh) 2022-09-27

Family

ID=76677915

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110408520.6A Active CN113096163B (zh) 2021-04-16 2021-04-16 一种无需先验升降轨信息的星载sar影像高精度匹配方法

Country Status (1)

Country Link
CN (1) CN113096163B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102213762A (zh) * 2011-04-12 2011-10-12 中交第二公路勘察设计研究院有限公司 基于rfm模型的多源星载sar影像自动匹配方法
CN105093222A (zh) * 2015-07-28 2015-11-25 中国测绘科学研究院 一种sar影像区域网平差连接点自动提取方法
CN107392951A (zh) * 2017-06-06 2017-11-24 上海卫星工程研究所 遥感图像高精度快速配准方法
CN109523585A (zh) * 2018-11-19 2019-03-26 武汉大学 一种基于方向相位一致性的多源遥感影像特征匹配方法
CN110111274A (zh) * 2019-04-28 2019-08-09 张过 一种星载推扫式光学传感器外方位元素定标方法
US20190251678A1 (en) * 2017-09-30 2019-08-15 Institute Of Remote Sensing And Digital Earth, Chinese Academy Of Sciences Automatic cross-platform geometric correction method for moon-based earth observation image

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102213762A (zh) * 2011-04-12 2011-10-12 中交第二公路勘察设计研究院有限公司 基于rfm模型的多源星载sar影像自动匹配方法
CN105093222A (zh) * 2015-07-28 2015-11-25 中国测绘科学研究院 一种sar影像区域网平差连接点自动提取方法
CN107392951A (zh) * 2017-06-06 2017-11-24 上海卫星工程研究所 遥感图像高精度快速配准方法
US20190251678A1 (en) * 2017-09-30 2019-08-15 Institute Of Remote Sensing And Digital Earth, Chinese Academy Of Sciences Automatic cross-platform geometric correction method for moon-based earth observation image
CN109523585A (zh) * 2018-11-19 2019-03-26 武汉大学 一种基于方向相位一致性的多源遥感影像特征匹配方法
CN110111274A (zh) * 2019-04-28 2019-08-09 张过 一种星载推扫式光学传感器外方位元素定标方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
C. WERNER 等: "SAR geocoding and multi-sensor image registration", 《IEEE INTERNATIONAL GEOSCIENCE AND REMOTE SENSING SYMPOSIUM》 *
E. SANSOSTI 等: "Geometrical SAR image registration", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
李莹莹 等: "基于星载光学和SAR 影像的多视立体定位误差分析", 《光学学报》 *
陈尔学 等: "尺度不变特征变换法在SAR影像匹配中的应用", 《自动化学报》 *

Also Published As

Publication number Publication date
CN113096163B (zh) 2022-09-27

Similar Documents

Publication Publication Date Title
CN111028277B (zh) 基于伪孪生卷积神经网络的sar和光学遥感图像配准方法
CN110097093B (zh) 一种异源图像精确匹配方法
CN110569861B (zh) 一种基于点特征和轮廓特征融合的图像匹配定位方法
CN103822616B (zh) 一种图分割与地形起伏约束相结合的遥感影像匹配方法
WO2019042232A1 (zh) 一种快速鲁棒的多模态遥感影像匹配方法和系统
CN111145228B (zh) 基于局部轮廓点与形状特征融合的异源图像配准方法
CN111507901B (zh) 基于航带gps及尺度不变约束的航拍图像拼接定位方法
CN114936971A (zh) 一种面向水域的无人机遥感多光谱图像拼接方法及系统
CN111369495B (zh) 一种基于视频的全景图像的变化检测方法
CN113222820B (zh) 一种位姿信息辅助的航空遥感图像拼接方法
CN105550994A (zh) 一种基于卫星影像的无人机影像快速概略拼接方法
CN112419374A (zh) 一种基于图像配准的无人机定位方法
CN107862319B (zh) 一种基于邻域投票的异源高分光学影像匹配误差剔除方法
CN115187798A (zh) 一种多无人机高精度匹配定位方法
CN109671109B (zh) 密集点云生成方法及系统
Koch et al. A new paradigm for matching UAV-and aerial images
Huang et al. SAR and optical images registration using shape context
Li et al. Adaptive regional multiple features for large-scale high-resolution remote sensing image registration
CN113096163B (zh) 一种无需先验升降轨信息的星载sar影像高精度匹配方法
CN115205564B (zh) 基于无人机的船体维护巡检方法
CN104484647B (zh) 一种高分辨率遥感图像云高度检测方法
CN116206139A (zh) 一种基于局部自卷积的无人机图像升尺度匹配方法
Ren et al. SAR image matching method based on improved SIFT for navigation system
CN114565653A (zh) 一种存在旋转变化和尺度差异的异源遥感图像匹配方法
CN114255398A (zh) 一种卫星视频图像的特征提取与匹配的方法及装置

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