发明内容
本发明为实现图像匹配策略自动评价遥感影像无控定位精度,本发明提供了一种基于参考底图的遥感影像定位精度评价方法,本发明提供了以下技术方案:
一种基于参考底图的遥感影像定位精度评价方法,包括以下步骤:
步骤1:提取待评价的影像,根据提取待评价的影像确定LandSat参考影像;
步骤2:对待评价的影像和LandSat参考影像预处理,采用配准SURF算法匹配参数,获取备选的控制点;
步骤3:以LandSat参考影像作为参考评价待评价的影像定位精度。
优选地,所述步骤1具体为:提取待评价的影像,根据待评价的影像的经纬度最大范围作为检索条件,查询LandSat参考影像数据库,选取覆盖待评价的影像的经纬度最大范围的LandSat参考影像作为待评价影像的初选控制数据。
优选地,所述步骤2具体为:
步骤2.1:对待评价影像和LandSat参考影像进行预处理,提高待评价影像和LandSat参考影像的信噪比,将待评价影像和LandSat参考影像的分辨率调整至不超过5倍的关系;
步骤2.2:将待评价影像的四角点及中心点通过坐标转换分别在LandSat参考影像上进行定位,设置分块的宽度和高度,宽度和高度分别设为1024个像素,对分块的子区域进行特征点提取与匹配;
步骤2.3:采用SURF算法对待评价影像和LandSat参考影像进行匹配,筛选最优匹配点,基于误匹配的点存在孤立特性,按同名点所占比例提取,剔除误匹配点。
优选地,所述步骤3具体为:
步骤3.1:根据分块区域内的匹配点,计算LandSat影像匹配点在待评价影像中的影像坐标,将LandSat影像的TM坐标转换为WGS84大地坐标,通过下式表示转换过程:
其中,B为转换后的纬度,L为转换后经度,e为椭球的第一偏心率,L0为原点经度,K为比例因子,YE为TM坐标的纵坐标,XN为TM坐标的横坐标;
利用转换后的LandSat影像平面坐标(B,L)在公开的参考DEM数据中内插该点的椭球高程值H,从而得到LandSat影像匹配点在WGS84椭球下的大地坐标(BLH);
步骤3.2:将(BLH)利用RPC正解求得在待评价影像上的图像坐标(X1,Y1),确定待评价影像的同名点,通过下式表示所述同名点:
Pn=(P-P0)/Ps;Ln=(L-L0)/Ls;Hn=(H-H0)/Hs
rn=(r-r0)/rs;cn=(c-c0)/cs
其中,P、L、H分别为WGS84坐标系下的纬度、经度和大地高,(rn,cn)、(Pn,Ln,Hn)分别为平移和缩放后的正则化像点坐标和地面点坐标,取值为[-1,1]之间,rs,cs,Ps,Ls,Hs为正则化的缩放系数;r0,c0,P0,L0,H0为正则化的平移系数,r和c为像点坐标所在的行列号;
步骤3.3:确定待评价影像的误差,通过下式表示所述误差:
其中,Δ为待评价影像的误差;
对所有同名点进行统计,确定所有同名点的误差,通过下式表示所有同名点的误差:
其中,σ为所有同名点的误差,n为同名点数量;
步骤3.4:根据所有同名点的误差对待测影像的定位精度极性评价,通过下式表示待测影像的定位精度:
Position=σ×GSD
其中,Position为待测影像的定位精度,GSD为地面像元分辨率。
本发明具有以下有益效果:
本发明可广泛用于国产光学遥感卫星影像的几何质量评价,该方法实施难度小,可以进行遥感影像自动、快速、高精度的评价具有很大实用价值,可以消除人工质检时刺点的繁琐、解决人工质检时效性差等问题。处理一景多光谱遥感影像耗时<25s,准确率优于90%。
具体实施方式
以下结合具体实施例,对本发明进行了详细说明。
具体实施例一:
步骤1:提取待评价的影像,根据提取待评价的影像确定LandSat参考影像;
所述步骤1具体为:提取待评价的影像,根据待评价的影像的经纬度最大范围作为检索条件,查询LandSat参考影像数据库,,选取覆盖待评价的影像的经纬度最大范围的LandSat参考影像作为待评价影像的初选控制数据。
步骤2:对待评价的影像和LandSat参考影像预处理,采用配准SURF算法匹配参数,获取备选的控制点;
所述步骤2具体为:
步骤2.1:对待评价影像和LandSat参考影像进行预处理,提高待评价影像和LandSat参考影像的信噪比,将待评价影像和LandSat参考影像的分辨率调整至范围内,分辨率最好不超过5倍关系;
步骤2.2:将待评价影像的四角点及中心点通过坐标转换分别在LandSat参考影像上进行定位,设置分块的宽度和高度,宽度和高度分别设为1024个像素,对分块的子区域进行特征点提取与匹配;
步骤2.3:采用SURF算法对待评价影像和LandSat参考影像进行匹配,筛选最优匹配点,误匹配的点会有孤立的特性,因此按同名点所占比例提取,剔除误匹配点;。
步骤3:以LandSat参考影像作为参考评价待评价的影像的定位精度。
所述步骤3具体为:
步骤3.1:根据分块区域内的匹配点,计算LandSat影像匹配点在待评价影像中的影像坐标,将LandSat影像的TM坐标转换为WGS84大地坐标,通过下式表示转换过程:
其中,B为转换后的纬度,L为转换后经度,e为椭球的第一偏心率,L0为原点经度,K为比例因子(也叫尺度变化),YE为TM坐标纵坐标,XN为TM坐标横坐标;
利用转换后的LandSat影像平面坐标(B,L)在公开的参考DEM数据中内插该点的椭球高程值H,从而得到LandSat影像匹配点在WGS84椭球下的大地坐标(BLH);
步骤3.2:将(BLH)利用RPC正解求得在待评价影像上的图像坐标(X1,Y1),确定待评价影像的同名点,通过下式表示所述同名点:
Pn=(P-P0)/Ps;Ln=(L-L0)/Ls;Hn=(H-H0)/Hs
rn=(r-r0)/rs;cn=(c-c0)/cs
式中,P、L、H分别为WGS84坐标系下的纬度、经度和大地高,单位分别为度和米。r、c为像点坐标所在的行列号。(rn,cn)、(Pn,Ln,Hn)分别为平移和缩放后的正则化像点坐标和地面点坐标,取值为[-1,1]之间。rs,cs,Ps,Ls,Hs为正则化的缩放系数;r0,c0,P0,L0,H0为正则化的平移系数,均可以从影像附带的90个RPC系数中获取。目前RPC系数有两种常用的格式,一种是_RPC.txt形式,常见于WordView、QuickBird、资源三号等系列卫星,另一种是RPB格式,这两种格式都是文本格式。多项式Pi(i=1,2,3,4)中每一项的各坐标分量Pn,Ln,Hn的幂次最大不超过3,且每一项各个坐标分量的幂次之和也不超过3。多项式Pi的形式如下:
p1(Pn,Ln,Hn)=a1+a2·Ln+a3·Pn+a4·Hn+a5·Ln·Pn+a6·Ln·Hn+a7·Pn·Hn+a8·Ln 2+a9·Pn 2+a10·Hn 2+a11·Pn·Ln·Hn+a12·Ln 3+a13·Ln·Pn 2+a14·Ln·Hn 2+a15·Ln 2·Pn+a16·Pn 3+a17·Pn·Hn 2+a18·Ln 2·Hn+a19·Pn 2·Hn+a20·Hn 3
p2(Pn,Ln,Hn)=b1+b2·Ln+b3·Pn+b4·Hn+b5·Ln·Pn+b6·Ln·Hn+b7·Pn·Hn+b8·Ln 2+b9·Pn 2+b10·Hn 2+b11·Pn·Ln·Hn+b12·Ln 3+b13·Ln·Pn 2+b14·Ln·Hn 2+b15·Ln 2·Pn+b16·Pn 3+b17·Pn·Hn 2+b18·Ln 2·Hn+b19·Pn 2·Hn+b20·Hn 3
p3(Pn,Ln,Hn)=c1+c2·Ln+c3·Pn+c4·Hn+c5·Ln·Pn+c6·Ln·Hn+c7·Pn·Hn+c8·Ln 2+c9·Pn 2+c10·Hn 2+c11·Pn·Ln·Hn+c12·Ln 3+c13·Ln·Pn 2+c14·Ln·Hn 2+c15·Ln 2·Pn+c16·Pn 3+c17·Pn·Hn 2+c18·Ln 2·Hn+c19·Pn 2·Hn+c20·Hn 3
p4(Pn,Ln,Hn)=d1+d2·Ln+d3·Pn+d4·Hn+d5·Ln·Pn+d6·Ln·Hn+d7·Pn·Hn+d8·Ln 2+d9·Pn 2+d10·Hn 2+d11·Pn·Ln·Hn+d12·Ln 3+d13·Ln·Pn 2+d14·Ln·Hn 2+d15·Ln 2·Pn+d16·Pn 3+d17·Pn·Hn 2+d18·Ln 2·Hn+d19·Pn 2·Hn+d20·Hn 3
式中,ai、bi、ci、di(i=1,2,3…,20)为多项式展开形式的系数即RPC系数。
步骤3.3:确定待评价影像的误差,通过下式表示所述误差:
其中,Δ为待评价影像的误差;
对所有同名点进行统计,确定所有同名点的误差,通过下式表示所有同名点的误差:
其中,σ为所有同名点的误差,n为同名点数量;
步骤4.3:根据所有同名点的误差对待测影像的定位精度极性评价,通过下式表示待测影像的定位精度:
Position=σ×GSD
其中,Position为待测影像的定位精度,GSD为地面像元分辨率。
具体实施例二:
本发明提取吉林一号光学A星21702_202的浙江德清县影像数据作为试验数据,影像数据包括经纬度范围以及RPC参数。该轨影像共获取19景影像,以其中第11景影像为例,该景影像的四角点及中心点经纬度分别为[119.967096,30.732424]、[120.089773,30.732424]、[119.937748,30.648389]、[120.060289,30.624353]、[120.013600,30.690445],以该影像覆盖的范围作为检索条件,搜索符合条件的LandSat参考影像。
步骤2,配准SURF算法匹配参数,获取备选控制点;
待评价和参考影像预处理:
本方法采用的参考影像LandSat8全色波段的分辨率为15m,吉林一号光学A星全色分辨率为0.72m,两者之间的分辨率相差太大,直接处理比较困难,需将二者之间的分辨率调整到一定范围内,因此对吉林一号光学A星全色影像进行降分辨率处理,降为10m。
步骤2.2控制点特征提取与匹配:
影像分块处理操作:
计算光学A星四角点及中心点在LandSat影像的位置,位置分别为(10583,13160),(11349,13328),(10408,13956),(11172,14123),(10878,13642)。对待评价的光学A星及LandSat影像进行分块处理,光学A星分块影像大小为512*512像素,LandSat分块影像大小为100*100。对分块的区域进行特征点提取等操作。
待评价影像与参考影像提取同名点:
利用配置的SURF算法匹配参数,对待评价的光学A星影像以及LandSat参考影像进行特征点提取并进行匹配,共匹配65对同名点,以其中3对为例进行说明,点对1:(1980,3071)与(6175.66,8728.63),点对2:(2620,2783)与(6310.48,8700.17),点对3:(1328,2795)与(6059.76,8648.47)。
最优匹配点筛选:
通过以上步骤获取的匹配特征点会有误匹配,需要进行筛选,最优匹配点的筛选可以通过按一定比例提取出排名靠前的匹配点来实现,最终选取13对排名靠前的同名点对。
先基于步骤2获取单块区域的匹配点,计算LandSat影像匹配点在待评价影像的图像坐标,由于LandSat影像采用地是TM坐标,所以需将TM坐标先转为WGS84大地坐标,点对1的LandSat影像的TM坐标为(6175.66,8728.63),转换(B,L)为(116.91,41.64),利用转换后的LandSat影像平面坐标(B,L)在公开的参考DEM数据中内插该点的椭球高程值H,从而得到LandSat影像匹配点在WGS84椭球下的大地坐标(BLH),(BLH)为(116.91,41.64,998.94)。将(BLH)利用RPC正解求得在待评价影像上的图像坐标(X1,Y1),(X1,Y1)为(2011.62,3056.47)。计算待评价影像同名点与(X1,Y1)的差值,计算该点的误差△,△为34.796,单位为像素。重复筛选出的所有最优匹配点对,计算所有的同名点对的差值的中误差σ,中误差σ为25.892,单位为像素。计算该景影像的定位精度Position,Position为77.41米。
以上所述仅是一种基于参考底图的遥感影像定位精度评价方法的优选实施方式,一种基于参考底图的遥感影像定位精度评价方法的保护范围并不仅局限于上述实施例,凡属于该思路下的技术方案均属于本发明的保护范围。应当指出,对于本领域的技术人员来说,在不脱离本发明原理前提下的若干改进和变化,这些改进和变化也应视为本发明的保护范围。