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

一种无需先验升降轨信息的星载sar影像高精度匹配方法
技术领域
1.本发明涉及通信领域，特别涉及一种无需先验升降轨信息的星载sar影像高精度匹配方法。


背景技术：

2.合成孔径雷达(sar)技术的不断进步，使人们从sar影像中获取更多地表信息成为可能。影像匹配是指通过一定的匹配算法在两幅或多幅影像之间识别同名点的过程。它是图像融合、目标识别、目标变化检测、计算机视觉等问题中的一个重要前期步骤，在遥感、数字摄影测量、计算机视觉、地图学及军事应用等多个领域都有着广泛的应用，对卫星影像的配准、拼接、3d重建等具有重要意义。
3.绝大多数针对星载影像的匹配都是同轨同侧sar影像之间的匹配，针对如异轨异侧等其它情况的星载sar影像匹配方法多是先将其中一景影像通过正射纠正的方式将其纠正到和另一景影像同一视角后再进行匹配，进行正射纠正后进行匹配，匹配到的同名点会增加正射纠正这一步所带来的累积误差。另外，星载影像较大，直接匹配结果较差，耗费的计算资源也较大。


技术实现要素：

4.针对现有技术存在的上述问题，本发明提供一种通过先行判断两景sar影像之间的成像轨道和观测方位之间的关系，然后直接进行匹配的方法，解决了非同轨同侧的两景影像需要正射纠正后匹配的问题，并通过分块匹配的方案对两景影像进行匹配。
5.为实现上述目的，本发明采取以下技术方案：一种无需先验升降轨信息的星载sar影像高精度匹配方法，其具体包括以下步骤：
6.步骤a)计算一景星载sar影像的四个顶点的地理坐标；
7.步骤b)根据步骤a)所得地理坐标计算其在另一景星载sar影像上的像平面坐标；
8.步骤c)根据步骤b)所得像平面坐标形成的四边形，选取其最大内接矩形，作为两景星载sar影像的重叠区域；
9.步骤d)判断两景星载sar影像成像轨道与观测方向之间的关系，并将步骤c)所得两景星载sar影像的重叠区域翻转为同轨同侧；
10.步骤e)进行两景星载sar影像的匹配。
11.优选地，步骤a)的具体操作为：根据一景星载sar影像的顶点坐标，设置初始高程值h0＝100，代入公式(1)和公式(2)计算一景星载sar影像该顶点对应的经纬度坐标，再根据全球高程模型获取该经纬度位置的实际高程值h1判断
△
h＝|h1‑
h0|的大小，当
△
h≥0.001时，再将高程值h1代入公式(1)和公式(2)进行迭代，直至
△
h＝|h
n
‑
h
n
‑1|＜0.001时求得其所对应的p、l后停止迭代，即得一景星载sar影像的顶点坐标对应的地理坐标(p,l,h
n
)。
12.更优选地，公式(1)和公式(2)具体为：
[0013][0014][0015]
其中，x表示像平面的横坐标，y表示像平面的纵坐标，h表示像平面坐标所在实际地理空间的高程，p表示像平面坐标所在实际地理空间的经度，l表示像平面坐标所在实际地理空间的纬度；且，
[0016]
num
p
(x,y,h)＝a1+a2y+a3x+a4h+a5yx+a6yh+a7xh+a8y2+a9x2+a
10
h2+a
11
xyh+a
12
y3+a
13
yx2+a
14
yh2+a
15
y2x+a
16
x3+a
17
xh2+a
18
y2h+a
19
x2h+a
20
h3[0017]
den
p
(x,y,h)＝b1+b2y+b3x+b4h+b5yx+b6yh+b7xh+b8y2+b9x2+b
10
h2+b
11
xyh+b
12
y3+b
13
yx2+b
14
yh2+b
15
y2x+b
16
x3+b
17
xh2+b
18
y2h+b
19
x2h+b
20
h3[0018]
num
l
(x,y,h)＝c1+c2y+c3x+c4h+c5yx+c6yh+c7xh+c8y2+c9x2+c
10
h2+c
11
xyh+c
12
y3+c
13
yx2+c
14
yh2+c
15
y2x+c
16
x3+c
17
xh2+c
18
y2h+c
19
x2h+c
20
h3[0019]
den
l
(x,y,h)＝d1+d2y+d3x+d4h+d5yx+d6yh+d7xh+d8y2+d9x2+d
10
h2+d
11
xyh+d
12
y3+d
13
yx2+d
14
yh2+d
15
y2x+d
16
x3+d
17
xh2+d
18
y2h+d
19
x2h+d
20
h3[0020]
其中，a
i
,b
i
,c
i
,d
i
；i＝1,2,
…
,20为星载sar影像的rpc参数。
[0021]
优选地，步骤b)的具体操作为：根据步骤a)所得地理坐标，和另一景星载sar影像的rpc参数，代入公式(3)和公式(4)求解一景星载sar影像的顶点对应于另一景星载sar影像上的像平面坐标。
[0022]
更优选地，公式(3)和公式(4)具体为：
[0023][0024][0025]
其中，
[0026]
num
x
(p,l,h)＝a1+a2l+a3p+a4h+a5lp+a6lh+a7ph+a8l2+a9p2+a
10
h2+a
11
plh+a
12
l3+a
13
lp2+a
14
lh2+a
15
l2p+a
16
p3+a
17
ph2+a
18
l2h+a
19
p2h+a
20
h3[0027]
dem
x
(p,l,h)＝b1+b2l+b3p+b4h+b5lp+b6lh+b7ph+b8l2+b9p2+b
10
h2+b
11
plh+b
12
l3+b
13
lp2+b
14
lh2+b
15
l2p+b
16
p3+b
17
ph2+b
18
l2h+b
19
p2h+b
20
h3[0028]
num
y
(p,l,h)＝c1+c2l+c3p+c4h+c5lp+c6lh+c7ph+c8l2+c9p2+c
10
h2+c
11
plh+c
12
l3+c
13
lp2+c
14
lh2+c
15
l2p+c
16
p3+c
17
ph2+c
18
l2h+c
19
p2h+c
20
h3[0029]
dem
y
(p,l,h)＝d1+d2l+d3p+d4h+d5lp+d6lh+d7ph+d8l2+d9p2+d
10
h2+d
11
plh+d
12
l3+d
13
lp2+d
14
lh2+d
15
l2p+d
16
p3+d
17
ph2+。d
18
l2h+d
19
p2h+d
20
h3[0030]
优选地，步骤c)计算最大内接矩形的具体操作为：分别对解算的横纵坐标按从小到大的顺序排序，得从小到大的横坐标顺序和纵坐标顺序，再根据所得横坐标顺序和纵坐标顺序组成其最大内接矩形。
[0031]
优选地，步骤d)的判断方法为：根据步骤a和步骤b计算一景星载sar影像的四个顶
点在其“同轨同侧”、“同轨异侧”、“异轨同侧”和“异轨异侧”的星载sar影像的像平面坐标，对一景星载sar影像的四个顶点进行标号并将其对应的像平面坐标同时对应进行标号，得到一景星载sar影像对应其不同成像轨道和观测方向之间的对应关系，再将待测两景星载sar影像的四个顶点和其对应像平面坐标的标号与所得对应关系进行比较，得两景星载sar影像的实际成像轨道和观测方向之间的对应关系。
[0032]
优选地，所述方法还可包括步骤d’)将两景星载sar影像分割成若干不重叠的小块影像对，以便节省计算过程中所需计算机的算力；所述步骤d’)可添加在步骤d)前或后。
[0033]
优选地，步骤e)包括如下步骤：
[0034]
步骤e1)将已翻转的两景星载sar影像(或分割后的若干组小块影像对)进行关键点检测，计算其特征描述子；
[0035]
步骤e2)根据步骤e1)计算所得特征描述子，对检测出的关键点进行暴力匹配；
[0036]
步骤e3)剔除步骤e2)匹配结果中的误匹配点。
[0037]
更优选地，步骤e1)具体为：检测步骤d’)所得一两景小块影像对或步骤d)所得两景星载sar影像的关键点，计算特征描述子。
[0038]
进一步，步骤e1)中检测两景星载sar影像或两景小块影像对的关键点的具体实现方式如下：
[0039]
⑤
计算影像水平方向和垂直方向的比值梯度：
[0040][0041][0042]
⑥
计算协方差矩阵：
[0043][0044]
⑦
计算角点响应值：
[0045]
r
sh
(x,y,α)＝det(g
sh
(x,y,α))
‑
d
·
(trace(g
sh
(x,y,α)))2ꢀꢀꢀ
(8)
[0046]
其中d的取值范围为0.04到0.08；
[0047]
⑧
通过角点响应局部极大值抑制，计算求得关键点坐标。
[0048]
进一步，步骤e1)中计算关键点特征描述子的具体实现方式如下：
[0049]
③
计算影像每个像素点的梯度与幅值：
[0050][0051][0052]
④
以关键点为中心，统计大小为16*16的邻域范围内的梯度方向直方图，将16*16大小的区域分为16块4*4大小的小区域，统计每个小区域8个梯度方向的幅值直方图，合并所有小区域形成128维的特征向量即为该关键点的特征描述子。
[0053]
更优选地，步骤e2)中，通过欧式距离(公式(11))的相似性度量，采用最近邻原则进行暴力匹配，得到初始匹配结果：
[0054][0055]
其中m
l
表示一景影像分割后的小块影像的特征描述子向量，m
r
表示另一景影像分割后的小块影像的特征描述子向量。
[0056]
更优选地，步骤e3)中剔除步骤e2)匹配结果中的误匹配点的具体实现方式如下：
[0057]
①
随机从匹配结果的同名点集中随机抽出4个样本数据，计算出变换矩阵：
[0058][0059]
记为模型m，
[0060]
②
计算同名点集中所有数据与模型m的投影误差，若误差小于阈值，加入内点集i，
[0061]
③
如果当前内点集i元素个数大于最优内点集i_best，则更新i_best＝i，同时更新迭代次数k，
[0062]
④
如果迭代次数大于k，则退出；否则迭代次数加1，并重复上述步骤。
[0063]
进一步，若步骤e3)是匹配若干小块影像对，需合并匹配得到的所有同名点，在全局情况下再次剔除误匹配点，剔除误匹配点的具体实现方式如下：
[0064]
①
随机从匹配结果的同名点集中随机抽出4个样本数据，计算出变换矩阵：
[0065][0066]
记为模型m，
[0067]
②
计算同名点集中所有数据与模型m的投影误差，若误差小于阈值，加入内点集i，
[0068]
③
如果当前内点集i元素个数大于最优内点集i_best，则更新i_best＝i，同时更新迭代次数k，
[0069]
④
如果迭代次数大于k＝1000，则退出；否则迭代次数加1，并重复上述步骤。
[0070]
与现有技术相比，本发明的优势如下：
[0071]
1、本发明无需知晓两景星载sar影像间的成像轨道与观测方向间的关系，通过星载影像附带的rpc参数信息计算可以得到两景影像间的关系。
[0072]
2、采用分块匹配的策略，节省计算过程中所需计算机的算力。
[0073]
3、通过自动判断两景影像间的关系，省略非同轨同侧影像匹配过程中的正射纠正步骤，省去了人为干预过程，大大节省人力成本。
附图说明
[0074]
图1为本发明提供的无需先验升降轨信息的星载sar影像高精度匹配方法的流程
图；
[0075]
图2为本发明提供的用于匹配的两景异轨同侧sar影像；
[0076]
图3为本发明提供的用于匹配的两景异轨同侧sar影像的重叠区域；
[0077]
图4为本发明提供的两景星载影像成像轨道和观测方向判定示例图；
[0078]
图5为本发明提供的用于判断两景sar影像为异轨同侧关系的实例图；
[0079]
图6为本发明提供的用于匹配的两景异轨同侧sar影像划分后的其中一对小块影像对；
[0080]
图7为本发明提供的用于匹配的两景异轨同侧sar影像划分后的其中一对小块影像对剔除误匹配点后的实际匹配效果图；
[0081]
图8为本发明提供的用于匹配的两景异轨同侧sar影像划分后的另一对小块影像对剔除误匹配点后的实际匹配效果图；
[0082]
图9为本发明提供的用于匹配的两景异轨同侧sar影像合并所有小块影像对同名点的匹配效果图；
[0083]
图10为本发明提供的用于匹配的两景异轨同侧sar影像最终剔除误匹配点后同名点的匹配效果图。
具体实施方式
[0084]
下面结合附图和实施例对本发明做进一步具体的说明。
[0085]
如图1所示，本发明提供的无需先验升降轨信息的星载sar影像高精度匹配方法，具体包括以下步骤：
[0086]
1)计算两景星载sar影像重叠区域：
[0087]
由于星载sar影像较大，直接匹配比较耗时且对计算机硬件要求较高，故确定两景星载sar影像的重叠区域有助于减小实际的匹配区域，如图2所示，为本实施例用于匹配的两景异轨同侧星载sar影像。
[0088]
现有技术中，常规的匹配过程需要搜索匹配区域，本发明省去了搜索匹配区域这一操作，通过公式(1)～公式(4)计算直接得到星载sar影像中明确的匹配范围，即通过两景星载sar影像包含的有理多项式系数参数(下称rpc参数)及全球高程模型图像文件，求解第一景星载sar影像已知的四个顶点(0,0)、(w
l
,0)、(w
l
,h
l
)、(0,h
l
)在第二景星载sar影像中对应的像平面坐标(x1′
,y1′
)、(x2′
,y2′
)、(x3′
,y3′
)、(x4′
,y4′
)，其中h
l
、w
l
分别表示第一景星载sar影像的高与宽，将第二景星载sar影像上述像平面坐标所组成的四边形的最大内接矩形作为两景星载sar影像的重叠区域。
[0089]
具体步骤如下：
[0090]
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。
[0091][0092][0093]
其中，x表示像平面的横坐标，y表示像平面的纵坐标，h表示像平面坐标所在实际地理空间的高程，p表示像平面坐标所在实际地理空间的经度，l表示像平面坐标所在实际地理空间的纬度；且，
[0094]
num
p
(x,y,h)＝a1+a2y+a3x+a4h+a5yx+a6yh+a7xh+a8y2+a9x2+a
10
h2+a
11
xyh+a
12
y3+a
13
yx2+a
14
yh2+a
15
y2x+a
16
x3+a
17
xh2+a
18
y2h+a
19
x2h+a
20
h3[0095]
den
p
(x,y,h)＝b1+b2y+b3x+b4h+b5yx+b6yh+b7xh+b8y2+b9x2+b
10
h2+b
11
xyh+b
12
y3+b
13
yx2+b
14
yh2+b
15
y2x+b
16
x3+b
17
xh2+b
18
y2h+b
19
x2h+b
20
h3[0096]
num
l
(x,y,h)＝c1+c2y+c3x+c4h+c5yx+c6yh+c7xh+c8y2+c9x2+c
10
h2+c
11
xyh+c
12
y3+c
13
yx2+c
14
yh2+c
15
y2x+c
16
x3+c
17
xh2+c
18
y2h+c
19
x2h+c
20
h3[0097]
den
l
(x,y,h)＝d1+d2y+d3x+d4h+d5yx+d6yh+d7xh+d8y2+d9x2+d
10
h2+d
11
xyh+d
12
y3+d
13
yx2+d
14
yh2+d
15
y2x+d
16
x3+d
17
xh2+d
18
y2h+d
19
x2h+d
20
h3[0098]
其中，a
i
,b
i
,c
i
,d
i
；i＝1,2,
…
,20为星载sar影像的rpc参数，此处为随高分3号卫星发来的第一星载sar影像所附带的rpc参数。
[0099]
将顶点坐标(0,0)与求解得到的准确高程值h
n
代入公式(1)和公式(2)，即可得到顶点坐标(0,0)对应的实际经纬度坐标(p,l)，与h
n
结合后最终得顶点坐标(0,0)对应的地理坐标(p,l,h
n
)。
[0100]
1.2根据第二景星载sar影像的rpc参数，和步骤1.1所得第一景星载sar影像顶点(0,0)对应的地理坐标(p,l,h
n
)，根据公式(3)和公式(4)求解其对应于第二景星载sar影像上的像平面坐标(x1′
,y1′
)：
[0101][0102][0103]
其中，
[0104]
num
x
(p,l,h)＝a1+a2l+a3p+a4h+a5lp+a6lh+a7ph+a8l2+a9p2+a
10
h2+a
11
plh+a
12
l3+a
13
lp2+a
14
lh2+a
15
l2p+a
16
p3+a
17
ph2+a
18
l2h+a
19
p2h+a
20
h3[0105]
dem
x
(p,l,h)＝b1+b2l+b3p+b4h+b5lp+b6lh+b7ph+b8l2+b9p2+b
10
h2+b
11
plh+b
12
l3+b
13
lp2+b
14
lh2+b
15
l2p+b
16
p3+b
17
ph2+b
18
l2h+b
19
p2h+b
20
h3[0106]
num
y
(p,l,h)＝c1+c2l+c3p+c4h+c5lp+c6lh+c7ph+c8l2+c9p2+c
10
h2+c
11
plh+c
12
l3+c
13
lp2+c
14
lh2+c
15
l2p+c
16
p3+c
17
ph2+c
18
l2h+c
19
p2h+c
20
h3[0107]
dem
y
(p,l,h)＝d1+d2l+d3p+d4h+d5lp+d6lh+d7ph+d8l2+d9p2+d
10
h2+d
11
plh+d
12
l3+d
13
lp2+d
14
lh2+d
15
l2p+d
16
p3+d
17
ph2+d
18
l2h+d
19
p2h+d
20
h3[0108]
其中，a
i
,b
i
,c
i
,d
i
；i＝1,2,
…
,20为星载sar影像的rpc参数，此处为随高分3号卫星发来的第二星载sar影像所附带的rpc参数。
[0109]
1.3重复步骤1.1和1.2分别求得第一景星载sar影像另外三个顶点(w
l
,0)、(w
l
,h
l
)、(0,h
l
)在第二景星载sar影像中的像平面坐标(x2′
,y2′
)、(x3′
,y3′
)、(x4′
,y4′
)。
[0110]
1.4由于步骤1.1～1.3所得第一景星载sar影像的四顶点在第二景星载sar影像中的像平面坐标(x1′
,y1′
)、(x2′
,y2′
)、(x3′
,y3′
)、(x4′
,y4′
)所形成的四边形(两景星载sar影像的实际重叠区域)图形不规则，故选取其最大内接矩形为两景星载sar影像的重叠区域。具体选择方式为：
[0111]
分别对解算的横纵坐标按从小到大的顺序排序，如对第一景星载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′
)。
[0112]
如图3中矩形框为选取的两景星载sar影像通过上述步骤1所述方法得到的实际重叠区域，该重叠区域约为待匹配两景星载sar影像的一半大小，实际匹配范围缩小，匹配效率明显提高。
[0113]
类似地，本实施例也可以通过步骤1.1～步骤1.3求解第二景星载sar影像四顶点(0,0)、(w
r
,0)、(w
r
,h
r
)、(0,h
r
)在第一景星载sar影像中的像平面坐标(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)，其中h
r
、w
r
分别表示第二景星载sar影像的高与宽。再通过步骤1.4的方式得到第二景星载sar影像在第一景星载sar影像中最大内接矩形的四顶点坐标，在此不再赘述。
[0114]
2)判断两景星载sar影像成像轨道与观测方向之间的关系：
[0115]
由于现有技术中多采用正射纠正的方式进行对同轨异侧、异轨同侧和异轨异侧的星载sar影像进行视角纠正，本发明根据预先确定好的两景星载sar影像成像轨道与观测方向之间的关系，即可将不是同轨同侧的两景星载sar影像进行翻转操作，人为将两景星载sar影像处理为“同轨同侧”，以避免正射纠正带来的累积误差。
[0116]
发明人将一星载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)所示。
[0117]
通过上述方法，可得到一景星载sar影像的四个顶点对应于另一景星载sar影像中选取的重叠区域的四个顶点位置的四种可能情况。若如图4(b)，即四顶点呈逆时针旋转且
①
号顶点在左下角，则两景影像为同轨同侧的关系；若如图4(c)，即四顶点呈顺时针旋转且
①
号顶点在右下角，则两景影像为同轨异侧的关系；若如图4(d)，即四顶点呈顺时针旋转且
①
号顶点在左上角，则两景影像为异轨同侧的关系；若如图4(e)，即四顶点呈逆时针旋转且
①
号顶点在右上角，则两景影像为异轨异侧的关系。
[0118]
判断影像成像轨道与观测方向之间的关系后，若为同轨同侧影像，在分割影像对时不进行操作；若为同轨异侧影像，在分割影像时对第二景影像进行左右翻转操作；若为异轨同侧影像，在分割影像时对第二景影像进行上下翻转操作；若为异轨异侧影像，在分割影像时对第二景影像进行旋转180度操作。
[0119]
本实施例中，通过步骤1计算第二景星载sar影像的顶点坐标在第一景星载sar影像中的像平面坐标并对其标号，如图5所示，通过与图2进行对比可知，本实施例的两景影像为异轨同侧，应对第二景影像进行上下翻转操作，以便后续匹配。
[0120]
3)将两景星载sar影像分割成小块影像对：
[0121]
通过步骤2)已经得到本实施例中两景影像为同轨异侧影像，并对第二景影像进行了上下翻转操作。从两景影像中分割出步骤1)选中的区域，记其重叠区域的大小分别为m
l
*n
l
、m
r
*n
r
，将第一景星载sar影像的重叠区域分割成不重叠的若干600*600大小的小块影像，将翻转后的第二景星载sar影像的重叠区域分割成不重叠的同等数量大小的小块影像，分割后的两景影像的小块影像根据小块影像在重叠区域中的位置一一对应，得到若干分割后的小块影像对，如图6所示。
[0122]
4)检测步骤3所得小块影像对关键点
[0123]
关键点检测是图像匹配过程的核心步骤之一，通过角点检测的思路进行关键点检测，具体步骤为：
[0124]
4.1计算小块影像对的各影像水平方向和垂直方向的比值梯度：
[0125][0126][0127]
其中，α为尺度因子，g
x,α
为影像在尺度α下水平方向的梯度，g
y,α
为影像在尺度α下垂直方向的梯度，log为对数计算符号，m
1,α
(i)为影像在尺度α下第i个像素点5
×
5邻域上半部分的均值，m
2,α
(i)为影像在尺度α下第i个像素点5
×
5邻域下半部分的均值，m
3,α
(i)为影像在尺度α下第i个像素点5
×
5邻域左半部分的均值，m
4,α
(i)为影像在尺度α下第i个像素点5
×
5邻域右半部分的均值。
[0128]
4.2计算协方差矩阵：
[0129][0130]
其中，g
sh
(x,y,α)为像素点(x,y)在尺度α下的协方差矩阵，为均值为0，标准差为α的高斯核，*为卷积运算符号。
[0131]
4.3计算角点响应值：
[0132]
r
sh
(x,y,α)＝det(g
sh
(x,y,α))
‑
d
·
(trace(g
sh
(x,y,α)))2ꢀꢀ
(8)
[0133]
其中r
sh
(x,y,α)为像素点(x,y)在尺度α下的响应值，d为根据现有技术取值，本实施例可以取d＝0.04。
[0134]
4.4求取关键点坐标
[0135]
r
sh
(x,y,α)＞0.8
ꢀꢀ
(9)
[0136]
通过选取5
×
5邻域内的最大角点响应值r
sh
(x,y,α)，并满足式(9)的点的位置，即为求得的关键点坐标。
[0137]
5)计算特征描述子
[0138]
根据步骤4)检测出的关键点，采用梯度直方图的思想计算其特征描述子。
[0139]
具体步骤如下：
[0140]
5.1计算星载sar影像的每个像素点的梯度与幅值：
[0141][0142][0143]
其中g
t,α
为所在像素点的梯度值，g
n,α
为所在像素点的梯度方向。
[0144]
5.2以步骤4.4所得关键点为中心，统计大小为16*16的邻域范围内的梯度方向直方图，将16*16大小的区域分为16块4*4大小的小区域，统计每个小区域8个梯度方向的幅值直方图，合并所有小区域形成128维的特征向量即为该关键点的特征描述子。
[0145]
6)关键点匹配
[0146]
根据公式(11)所得欧式距离的相似性度量，采用最近邻原则进行暴力匹配，得初始匹配结果；
[0147][0148]
其中m
l
为第一景影像分割后的小块影像的特征描述子向量，m
r
为第二景影像分割后的对应的小块影像的特征描述子向量，d
m
为两匹配点描述子之间的欧式距离。
[0149]
7)误匹配点剔除
[0150]
初步匹配结果中含有较多误匹配点对，为了得到较为准确的同名点对，通过多次采样迭代的方案，找到最多匹配点子集即为剔除误匹配点后的准确的同名点集。具体步骤如下：
[0151]
7.1随机从步骤6)所得匹配结果的同名点集中随机抽出4个样本数据(此4个样本之间不能共线)，计算出变换矩阵：
[0152][0153]
记为模型m，其中t为求得的透视变换矩阵，t为透视变换矩阵中的具体参数。
[0154]
7.2计算步骤6)所得匹配结果的同名点集中所有数据与模型m的投影误差，若误差小于阈值，加入内点集i，
[0155]
7.3如果步骤7.2所得内点集i中元素个数大于最优内点集i_best，则更新i_best＝i，同时更新迭代次数k，
[0156]
7.4如果迭代次数大于k＝1000，则退出；否则迭代次数加1，并重复上述步骤7.1～7.3。
[0157]
图7、图8为采取步骤4)到步骤7)匹配的其中两对小块影像对的匹配结果。
[0158]
8)对所有小块影像对进行匹配
[0159]
步骤4)到步骤7)仅是一次匹配的过程，在本发明中，由于采用分块匹配的策略，因
此需要重复步骤4)到步骤7)，直到得到所有划分的小块影像对的匹配结果，结果图如图9所示(图9中的线条为匹配后同名点的连线)。
[0160]
9)合并所有同名点对，再次剔除误匹配点
[0161]
虽然每对小块影像的匹配结果较为准确，但其仅针对于小块影像所在的范围，而合并后的同名点集，其针对实际匹配的星载影像的重叠区域，其中某些同名点的误差也较大，因此需要再次剔除误差较大的匹配点，具体操作步骤同步骤7)，得最终匹配所得同名点分布图(图10中的线条为匹配后所有同名点的连线)，如图10。
[0162]
实施例2
[0163]
本实施例与实施例1的差别仅在于，删除步骤3、8、9，此方法为直接匹配而不需要对两景星载sar影像进行分割，对计算机要求过高，可按需进行操作。
[0164]
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。