CN102968631B - 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法 - Google Patents

Abstract

本发明公开了一种山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，包括以下步骤：基准影像的选择与拼接；同名地物点的自动选择；同名地物点的筛选、纠正影像的几何纠正与正射校正；精度评价；影像输出。山区多光谱遥感卫星影像如中国国产环境减灾卫星(型号：HJ-1A/B)影像，具有几何畸变严重的特点，本发明通过引入同名地物点的概念进行自动几何纠正与正射校正方法，能够实现批量化完成山区多光谱遥感卫星影像的正射纠正，并生成相关影像的残差报告，具有高效率、高精度的优点，对海量数据的山区多光谱遥感卫星影像处理非常有效，显著节省了传统遥感影像处理过程所需的人力物力，尤其对中国国产环境减灾卫星影像的处理非常有利。

Description

山区多光谱遥感卫星影像的自动几何纠正与正射校正方法

技术领域

本发明涉及一种山区多光谱遥感卫星影像的几何纠正与正射校正方法，尤其涉及一种航空航天领域所需的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，属于山区多光谱遥感卫星影像的计算机处理领域。

背景技术

遥感卫星影像的几何纠正与正射校正是遥感影像预处理的基础性工作，也是后期应用遥感卫星影像开展如生态环境监测、灾害监测与评估等方面的关键步骤。实际应用中，遥感卫星影像的几何畸变比较严重，图1示出了典型遥感卫星影像的几何畸变的原理，该图可广泛见于教科书中，如高等教育出版社出版、梅安新主编的《遥感导论》，科学出版社出版一书的第二版(2001年7月)，赵英时主编的《遥感应用分析原理与方法》一书第一版(2003年6月)第36页及第178页就列出了该图。从图中可以看出，这些几何畸变有的是由于卫星的姿态、轨道、地球的运动和形状等外部因素引起的；有的是由于传感器本身结构性能和扫描镜的不规则运动、检测器采样延迟、探测器的配置等外部因素所引起的；也有的则由于纠正上述误差而进行一系列换算和模拟而产生的处理误差。这些误差有的是系统的，有的是随机的，有的是连续的，有的是非连续的，十分复杂，它们构成了遥感卫星影像所固有的结合特性。遥感卫星图像的定位，一般是由卫星跟踪系统所提供的卫星轨道和姿态参数，根据卫星轨道公式进行的。由于所提供的卫星参数尚不足以精确的确定每个象元的地理位置，所以影像中的几何定位残差非常大。人们若需要得到高定位精度的遥感数据，需通过地面控制点信息，研究和改进遥感图像的几何精纠正方法，纠正影像中的几何误差。

然而，图1仅仅描述了卫星遥感影像的几何畸变，由于偏离星下观测地形导致的定位畸变没有给出描述。在山区，除了以上的几何畸变外，还存在着由于偏离卫星星下观测的起伏地形造成的山体定位误差，即山区多光谱遥感卫星影像还存在着由于偏离卫星星下观测的起伏地形造成的山体定位误差。卫星偏离星下观测导致的几何定位误差国内文献较少报道。美国J.Storey于2006年发表的LANDSAT7ImageAssessmentSystem(IAS)geometricAlgorithmTheoreticalBasisDocument一文的第117页中关于Landsat(地球资源卫星)上搭载的TM(ThematicMapper，专题制图仪)传感器的文章里提到正射校正的算法，并给出了由于卫星偏离星下观测导致的几何畸变规律，见图2的遥感卫星影像地形校正的几何原理示意图。然而，该原理仅适用于卫星垂直向下观测的规律。中国最新发射的环境减灾小卫星星座具有重复访问周期短、影像覆盖面积大的优势，该星座具有2颗光学卫星(型号：HJ1A/B)，卫星上各搭载两颗同底安装的CCD相机，编号分别为CCD1与CCD2，相机之间呈30°夹角。由于相机间存在30°夹角，因此，环境减灾卫星获取的影像为侧视影像。其成像规律与传统的垂直向下观测TM(ThematicMapper，专题制图仪)、CBERS(中巴地球资源卫星)等影像不同，影像的几何定位误差分布规律也有很大区别。中国国产环境减灾卫星同时获取的CCD1和CCD2影像共享星下基线，因此，星下基线位于CCD1与CCD2获取影像的重叠区域。沿着星下基线，降轨工作的CCD1影像位于右侧，地物坐标受地形影响自西向东误差逐渐增大，CCD2的误差分布规律则相反，而传统的单传感器影像畸变规律是由影像中心向两侧逐渐增大，见图3所示中国国产环境减灾卫星拍摄多光谱卫星遥感影像的获取影像示意。同时，国产环境减灾卫星(HJ-1A/B)的设计幅宽为360公里，而TM的设计幅宽为185公里，由于TM的星下基线位于影像中心，而国产环境减灾卫星(HJ-1A/B)位于影像边缘，因此，国产环境减灾卫星(HJ-1A/B)由中心至影像远端边界的距离为360公里，而TM仅93公里左右，环境减灾卫星的畸变更大。根据图4所示的地形及星下点基线对象元坐标精度影响拟合图进行理论分析，在环境减灾卫星影像边界，当象元点海拔为2000m时，环境减灾卫星影像由于偏离星下观测而造成的误差将达到1km，约合33个象元。因此，在校正之前，需要确认获取的国产环境减灾卫星(HJ-1A/B)影像的升降轨类型及其CCD相机编号，以了解其影像的基本几何误差分布规律。侧视影像的景中心位置与垂直观测的影像不同，环境减灾卫星降轨影像CCD1相机影像景中心位置偏左，CCD2偏右，见图3所示。因此，采用传统的算法原理开展几何畸变很大、导致原因复杂的山区侧视多光谱遥感卫星影像如中国国产环境减灾卫星拍摄的多光谱遥感卫星影像的几何纠正与正射校正，无法准确的纠正影像中的畸变，尤其是地形引起的几何畸变。

发明内容

本发明的目的就在于为了解决上述问题而提供一种能够高效率、高精度处理的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法。

为了达到上述目的，本发明采用了以下技术方案：

本发明所述山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，包括以下步骤：

(1)基准影像的选择与拼接：a、基准影像的选择：基准影像满足以下条件：整体云覆盖量小于5％；不同平台影像及邻近影像接边处配准误差小于0.5个象元；b、基准影像的拼接：拼接基准影像的近红外波段；

(2)同名地物点的自动选择：a、在基准影像和纠正影像上分别布设核窗口和搜索窗口，所述纠正影像为待纠正和校正的山区多光谱遥感卫星影像，下述亦然；b、将基准影像上核窗口内的关联点投影至纠正影像，然后以关联点为中心，计算纠正影像与基准影像的核窗口内的相关系数；c、在待纠正影像的搜索窗口内移动该核窗口，并计算相应核窗口内的相关系数；d、当搜索窗口内的核窗口全部计算后，若存在关联点核窗口的相关系数大于设定的阈值条件，且搜索窗口满足条件的关联点小于关联点数量上限阈值条件时，将相关系数最大的点认为是该窗口内的同名地物点；否则该搜索窗口内未找到同名地物点；

(3)同名地物点的筛选、纠正影像的几何纠正与正射校正：a、同名地物点的筛选：对纠正影像进行边缘检测，估算出星下基线斜率，然后根据从纠正影像头文件获取的影像中心点坐标与卫星的设计高度，估算出纠正影像中心点坐标对应的星下点坐标，进而根据点斜式线性方程估算出纠正影像星下基线方程截距；在获取星下基线方程后，将该基线方程投影至基准影像，结合基准影像上同名地物点对应数字高程模型DEM的高程信息与卫星高度，计算每个基准影像上同名地物点偏离星下点的地形偏移，叠加上每个同名地物点的对应的偏移以恢复成为与纠正影像一致的观测条件；b、纠正影像的几何纠正：所有同名地物点拟合获取的多项式均方根误差不能超过0.75个象元，否则去除拟合误差最大的点重新拟合多项式系数，并再次计算拟合多项式的均方根误差，同名地物点的预测几何定位误差大于1个象元时，去除该点，当所有同名地物点均满足纠正要求时，根据满足筛选条件控制点拟合获取多项式系数，对纠正影像的每个象元位置进行重新计算；c、纠正影像的正射校正：结合纠正影像的星下基线方程、卫星高度及数字高程模型DEM的数据，根据每个象元经过几何纠正后新的坐标信息，计算每个象元距星下基线的垂直距离，再计算每个象元由于偏离星下观测而导致的坐标偏差，去除由于偏离星下观测导致的偏差而完成正射校正；d、计算完毕后，对原始影像的光谱特征进行重采样；

(4)精度评价：通过光谱匹配，在基准影像和输出的正射校正结果之间搜索同名地物点作为验证点，将验证点按照影像中心分为左上、右上、左下、右下四个象限，分别判断每个象限内同名地物点配准误差，仅当四个区域的验证点比例均大于60％时，报告当前影像通过验证，否则增加配准多项式阶数配准后返回步骤(3)重新进行几何纠正与正射校正；

(5)影像输出：当达到程序收敛条件后，以纠正影像的原始影像范围输出经过几何纠正与正射校正的影响；

所述步骤(3)中，“估算纠正影像中心点坐标对应的星下点坐标”、“对纠正影像的每个象元位置进行重新计算”和“计算每个象元由于偏离星下观测而导致的坐标偏差”均采用以下公式：

$\left\{\begin{array}{c}\mathrm{dis}=\mathrm{Alt}*\tan \mathrm{\θ}\\ \mathrm{\Δx}=\mathrm{dis}*\sin \mathrm{\α}\\ \mathrm{\Δy}=\mathrm{dis}*\cos \mathrm{\α}\\ \mathrm{wcx}=\mathrm{cx}-\mathrm{\Δx}\\ \mathrm{wcy}=\mathrm{cy}-\mathrm{\Δy}\end{array}\right.$ 式（III)

式III中，dis为景中心位置距星下基线的距离，Alt为卫星高度，θ为卫星侧视角度，α为星下基线倾斜角度，Δx和Δy分别为景中心位置坐标偏移量，cx和cy为影像景中心位置，该中心位置通过与影像一起分发的XML头文件信息scenecenterlong和sceneeenterlat获取，(wcx，wcy)为纠正影像中心点坐标对应的星下点坐标。

作为优选，所述步骤(2)中，将基准影像上核窗口内的关联点按下述公式的坐标投影至纠正影像：

$\left\{\begin{array}{c}\mathrm{wx}=\frac{\mathrm{bx}*\mathrm{re}{s}_{\mathrm{base}}+({\mathrm{ulx}}_{\mathrm{base}}-\mathrm{ul}{x}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\\ \mathrm{wy}=\frac{\mathrm{by}*\mathrm{re}{s}_{\mathrm{base}}+({\mathrm{uly}}_{\mathrm{base}}-{\mathrm{uly}}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\end{array}\right.$ (式I)

式I中，(bx，by)为基准影像的坐标，(wx，wy)为纠正影像的估计坐标，(ulx_base，uly_base)为基准影像左上角象元坐标，(ulx_warp，uly_warp)为纠正影像左上角坐标，res_base和res_warp为基准影像和纠正影像的空间分辨率。

所述步骤(2)中，按下述公式计算纠正影像与基准影像的核窗口内的相关系数：

$r_{\mathrm{mn}}=\frac{\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}(A_{\mathrm{ij}}-\stackrel{\‾}{A})(B_{\mathrm{ij}}-\stackrel{\‾}{B})}{\sqrt{\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(A_{\mathrm{ij}}-\stackrel{\‾}{A})}^2\right)\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(B_{\mathrm{ij}}-\stackrel{\‾}{B})}^2\right))}}$ (式II)

式II中，(m，n)为核窗口中心坐标，A_ij为基准影像的象元光谱值；B_ij为纠正影像的象元光谱值，为基准影像核窗口内光谱平均值；为纠正影像核窗口内光谱平均值，cs代表核窗口大小。

所述步骤(2)中，设定的核窗口大小为7×7个象元，搜索窗口大小为50×50个象元。

所述步骤(2)中，阈值为0.75，设定的搜索窗口内满足阈值条件的最大个数为3个。

所述步骤(3)中，采用以下公式计算每个同名地物点偏离星下点的地形偏移：

$\left\{\begin{array}{c}S=\mathrm{pixel}\_\mathrm{size}*\mathrm{off}\_\mathrm{nadir}\_\mathrm{pixel}\\ s=S/\mathrm{Re}\\ \mathrm{LOS}=\mathrm{sqrt}({\mathrm{Re}}^2+{(\mathrm{Re}+\mathrm{Alt})}^2-2*\mathrm{Re}(\mathrm{Re}+\mathrm{Alt})*\cos \left(s\right))\\ d=a\sin (\mathrm{Re}*\sin \left(s\right)/\mathrm{LOS})\\ a\tan \{(\mathrm{Re}+\mathrm{Alt})*\sin \left(d\right)*(1-(\mathrm{Re}+h)/\mathrm{Re})/[\mathrm{Re}+h]*\mathrm{sqrt}(1-{(\mathrm{Re}+\mathrm{Alt})}^2\\ *\mathrm{si}{n}^2\left(d\right))/{\mathrm{Re}}^2)-(\mathrm{Re}+\mathrm{Alt})*\cos \left(d\right)\left]\right\}\\ z^{\′\′}=a\sin [(\mathrm{Re}+\mathrm{Alt})*\sin (d+\mathrm{dd})/\mathrm{Re}]\\ \mathrm{ds}=z^{\′\′}-s-(d+\mathrm{dd})\\ \mathrm{dS}=\mathrm{Re}*\mathrm{ds}\\ \mathrm{terrain}\_\mathrm{offset}=\mathrm{dS}/\mathrm{pixel}\_\mathrm{size}\end{array}\right.$ (式IV)

式IV中，S表示象元点距离星下基线的距离，pixel_size代表象元大小，off_nadir_pixel代表待校正象元距离星下的象元个数，s代表地球中心角度，Re代表地球参考椭球体在影像中心位置的地球半径，Alt为卫星的椭球体上空高度，terrain_offet为地形引起的位移。

所述步骤(3)中，根据以下公式拟合获取的多项式：

$\left\{\begin{array}{c}{x}^{\′}=a_0+a_{1}X+a_{2}Y+a_3{X}^2+a_4\mathrm{XY}+a_5{Y}^2+...\\ y^{\′}=b_0+b_{1}X+b_{2}Y+b_3{X}^2+b_4\mathrm{XY}+b_5{Y}^2+...\end{array}\right.$ (式V)

式V中，(x′，y′)为图像坐标，(X，Y)为参考坐标，a_i，j，b_i，j为多项式的系数。

所述步骤(3)中，重采样方法为最邻近法、二次内插法或三次卷积法。

上述方法中，在通过影像匹配算法获取足够的同名地物点下，当设定对应的CCD相机编号时，多光谱遥感卫星影像将通过控制点拟合得到的多项式进行纠正，整个流程进入到几何校正的状态。当几何校正完成后，根据其获得的星下基线与对应的数字高程模型数据，进一步开展正射校正，当正射校正完成后，程序将通过影像匹配的方法再次在基准影像与初校正结果之间搜索精度验证点，当影像的四个象限的精度均满足要求时，输出校正结果，否则增加多项式的阶数，再次重复进行，直至满足阈值条件为止，以达到对环境减灾卫星多光谱CCD遥感影像精确几何定位的目的。

本发明的有益效果在于：

本发明通过引入同名地物点的概念进行自动几何纠正与正射校正，能够实现批量化自动完成山区多光谱遥感卫星影像的正射纠正，并生成相关影像的残差报告，具有高效率、高精度、低成本的优点，对海量数据的山区多光谱遥感卫星影像处理非常有效，解决了传统多项式几何纠正把遥感卫星图像的总体变形看作是平移、缩放、旋转、扭曲、弯曲等基本变形综合作用的结果而导致的失真问题，显著节省了传统遥感影像处理过程所需的人力物力，尤其在利用中国国产环境减灾卫星影像于国土、资源清查、自然灾害应急响应调查、生态环境健康监测等方面具有十分重要的意义。

附图说明

图1是本发明背景技术所述遥感卫星影像几何畸变的原理示意图；

图2是本发明背景技术所述遥感卫星影像地形校正的几何原理示意图；

图3是本发明背景技术所述中国国产环境减灾卫星拍摄多光谱卫星遥感影像的获取影像示意图；

图4是本发明背景技术所述地形及星下点基线对象元坐标精度影响拟合图；

图5是本发明技术方案所述山区多光谱遥感卫星影像的自动几何纠正与正射校正方法的流程图；

图6是本发明技术方案所述核窗口和移动搜索窗口的关系示意图；

图7是本发明技术方案所述单台CCD相机构象几何特征分析示意图；

图8是本发明技术方案所述环境减灾卫星星下基线的计算原理示意图；

图9A是本发明实施例中的待校正影像示意图；

图9B是本发明实施例中的基准影像示意图；

图10A、图10B和图10C是本发明实施例中的光谱匹配同名点分布示意图；

图11是本发明实施例中的控制点筛选随星下基线位置变化情况示意图；

图12A和图12B是本发明实施例中的影响局部区域的纠正前、后对比图；

图13是本发明实施例中的精度检验控制点直方图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步具体描述：

下面以中国国产环境减灾卫星(型号：HJ-1A/B)影像的处理为例，对本发明进行具体说明：

如图5所示，本发明所述方法包括以下步骤：

(1)基准影像的选择与拼接：a、基准影像的选择：基准影像满足以下条件：整体云覆盖量小于5％；不同平台影像及邻近影像接边处配准误差小于0.5个象元；b、基准影像的拼接：拼接基准影像的近红外波段。

满足上述两个条件的原因如下：基准影像的选择是决定国产环境减灾卫星(HJ-1A/B)影像几何精纠正与正射校正的关键。当前公认地理位置较为准确的基准影像包括美国马里兰大学GLCF(GlobalLandCoverFacility)、USGS(美国地质调查局)、中国科学院对地观测中心等单位的数据平台上经过正射校正的Landsat(地球资源卫星)TM(专题制图仪)数据。然而，这些平台上有部分数据仍然存在较大的几何配准误差，该部分影像几何位置偏移较大，不能作为基准影像。此外，基准影像中云较多时，也将影响到控制点的自动匹配导致控制点空间分布不均匀而严重影响校正结果。在依据上述条件选取合适的基准影像后，由于我国环境减灾卫星的幅宽大于TM影像，因此在选取基准影像时应将环境减灾卫星完全覆盖。本方案要求对基准影像进行拼接处理，形成大幅宽的基准影像数据集。由于国产环境减灾卫星(HJ-1A/B)影像轨道存在一定漂移，在进行同轨道数据批处理时，需选取最大的轨道覆盖范围，同时，由于在几何精校正过程中采用的是灰度匹配算法，因此，本方案要求仅拼接基准影像的第四波段即近红外波段即可。选取近红外波段作为基准影像波段的原因是，该波段受大气影响较小，易于寻找同名地物点。

(2)同名地物点的自动选择：同名地物点的准确选择是几何精校正的基础，其选取精度是影响几何校正精度的关键。传统人工选择方法费时费力，需要大量的人工参与。本方案设计基于光谱匹配控制点搜索算法设计环境星的配准策略，是一种自动化的控制点筛选方法，能够有效节省人力物力。

本步骤的具体方法如下：如图6所示，a、在基准影像和纠正影像上分别布设核窗口和搜索窗口，设定的核窗口大小为7×7个象元，搜索窗口大小为50×50个象元，所述纠正影像为待纠正和校正的山区多光谱遥感卫星影像，下述亦然；b、将基准影像上核窗口内的关联点投影至纠正影像，然后以关联点为中心，计算纠正影像与基准影像的核窗口内的相关系数；c、在待纠正影像的搜索窗口内移动该核窗口，并计算相应核窗口内的相关系数；d、当搜索窗口内的核窗口全部计算后，若存在关联点核窗口的相关系数大于设定的阈值条件，且搜索窗口满足条件的关联点小于关联点数量上限阈值条件时，将相关系数最大的点认为是该窗口内的同名地物点，阈值为0.75，设定的搜索窗口内满足阈值条件的最大个数为3个；否则该搜索窗口内未找到同名地物点。

本步骤中，将基准影像上核窗口内的关联点按下述公式的坐标投影至纠正影像：

$\left\{\begin{array}{c}\mathrm{wx}=\frac{\mathrm{bx}*\mathrm{re}{s}_{\mathrm{base}}+({\mathrm{ulx}}_{\mathrm{base}}-\mathrm{ul}{x}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\\ \mathrm{wy}=\frac{\mathrm{by}*\mathrm{re}{s}_{\mathrm{base}}+({\mathrm{uly}}_{\mathrm{base}}-{\mathrm{uly}}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\end{array}\right.$ (式I)

式I中，(bx，by)为基准影像的坐标，(wx，wy)为纠正影像的估计坐标，(ulx_base，uly_base)为基准影像左上角象元坐标，(uls_warp，uly_warp)为纠正影像左上角坐标，res_base和res_warp为基准影像和纠正影像的空间分辨率。

按下述公式计算纠正影像与基准影像的核窗口内的相关系数：

$r_{\mathrm{mn}}=\frac{\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}(A_{\mathrm{ij}}-\stackrel{\‾}{A})(B_{\mathrm{ij}}-\stackrel{\‾}{B})}{\sqrt{\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(A_{\mathrm{ij}}-\stackrel{\‾}{A})}^2\right)\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(B_{\mathrm{ij}}-\stackrel{\‾}{B})}^2\right))}}$ (式II)

式II中，(m，n)为核窗口中心坐标，A_ij为基准影像的象元光谱值；B_ij为纠正影像的象元光谱值，为基准影像核窗口内光谱平均值；为纠正影像核窗口内光谱平均值，cs代表核窗口大小。

(3)同名地物点的筛选、纠正影像的几何纠正与正射校正：由于本方案选取的基准影像已经做过正射校正，因此，选取的同名地物点在基准影像上也已经做过了正射校正，而同名地物点在待校正影像上有地形误差，不能直接用于几何精纠正多项式系数的拟合，需要将基准影像中正射校正后的同名地物点投影成与待校正影像的观测角度、卫星高度一致条件下的地形误差的坐标点，以与待校正影像相匹配。因此，几何精校正与正射校正的前提要求是需要准确计算出每景影像的轨道信息，筛选出精度满足要求的同名地物点。本方案解决该问题的具体方法是根据卫星的飞行高度、地形海拔高度及传感器的侧视角度形成的几何关系进行基准影像控制点坐标位置的恢复，其几何关系如图6所示。

本步骤具体方法如下：a、如图8所示，同名地物点的筛选：对纠正影像进行边缘检测，估算出星下基线斜率，然后根据从纠正影像头文件获取的影像中心点坐标(cx，cy)(经纬度)与卫星的设计高度，估算出纠正影像中心点坐标对应的星下点坐标(wcx，wcy)，进而根据点斜式线性方程估算出纠正影像星下基线方程截距；在获取星下基线方程后，将该基线方程投影至基准影像，结合基准影像上同名地物点对应数字高程模型DEM(DigitalElevationModel的简称)的高程信息与卫星高度，计算每个基准影像上同名地物点偏离星下点的地形偏移，叠加上每个同名地物点的对应的偏移以恢复成为与纠正影像一致的观测条件；b、纠正影像的几何纠正：所有同名地物点拟合获取的多项式均方根误差不能超过0.75个象元，否则去除拟合误差最大的点重新拟合多项式系数，并再次计算拟合多项式的均方根误差，同名地物点的预测几何定位误差大于1个象元时，去除该点，当所有同名地物点均满足纠正要求时，根据满足筛选条件控制点拟合获取多项式系数，对纠正影像的每个象元位置进行重新计算；c、纠正影像的正射校正：结合纠正影像的星下基线方程、卫星高度及数字高程模型DEM的数据，根据每个象元经过几何纠正后新的坐标信息，计算每个象元距星下基线的垂直距离，再计算每个象元由于偏离星下观测而导致的坐标偏差，去除由于偏离星下观测导致的偏差而完成正射校正；d、计算完毕后，对原始影像的光谱特征进行重采样，重采样方法为最邻近法、二次内插法或三次卷积法。

本步骤的具体操作为，根据影像的头文件中影像景中心位置信息(sceneLong，SceneLati)与卫星坐标系统信息，将其转化为中心位置的行列号；由于HJ影像上是双CCD相机，单幅CCD数据是侧视成像，同期获取的2景CCD影像共享星下基线。因此，需根据影像相机位置信息，结合影像景中心位置，加入对应角度信息计算国产环境减灾卫星(HJ-1A/B)影像的星下基线。根据中心位置行列号及高度信息，结合传感器编号及其侧视角度，设定CCD1的侧视角度为15°，CCD2的侧视角度为-15°，估算景中心对应的影像星下观测的坐标位置点；通过扫描影像的有效边界，获取影像边界斜率，结合星下坐标位置点，准确拟合卫星在影像飞行的星下基线位置，进而获取每景影像上的星下基线方程(见图8)；在获取影像的星下基线方程后，将该星下基线方程投影至基准影像上。

本步骤中，“估算纠正影像中心点坐标对应的星下点坐标”、“对纠正影像的每个象元位置进行重新计算”和“计算每个象元由于偏离星下观测而导致的坐标偏差”均采用以下公式：

$\left\{\begin{array}{c}\mathrm{dis}=\mathrm{Alt}*\tan \mathrm{\θ}\\ \mathrm{\Δx}=\mathrm{dis}*\sin \mathrm{\α}\\ \mathrm{\Δy}=\mathrm{dis}*\cos \mathrm{\α}\\ \mathrm{wcx}=\mathrm{cx}-\mathrm{\Δx}\\ \mathrm{wcy}=\mathrm{cy}-\mathrm{\Δy}\end{array}\right.$ (式III)

式III中，dis为景中心位置距星下基线的距离，Alt为卫星高度，θ为卫星侧视角度，α为星下基线倾斜角度，Δx和Δy分别为景中心位置坐标偏移量，cx和cy为影像景中心位置，该中心位置通过与影像一起分发的XML头文件信息scenecenterlong和scenecenterlat获取，(wcx，wcy)为纠正影像中心点坐标对应的星下点坐标。

采用以下公式计算每个同名地物点偏离星下点的地形偏移：

$\left\{\begin{array}{c}S=\mathrm{pixel}\_\mathrm{size}*\mathrm{off}\_\mathrm{nadir}\_\mathrm{pixel}\\ s=S/\mathrm{Re}\\ \mathrm{LOS}=\mathrm{sqrt}({\mathrm{Re}}^2+{(\mathrm{Re}+\mathrm{Alt})}^2-2*\mathrm{Re}(\mathrm{Re}+\mathrm{Alt})*\cos \left(s\right))\\ d=a\sin (\mathrm{Re}*\sin \left(s\right)/\mathrm{LOS})\\ a\tan \{(\mathrm{Re}+\mathrm{Alt})*\sin \left(d\right)*(1-(\mathrm{Re}+h)/\mathrm{Re})/[\mathrm{Re}+h]*\mathrm{sqrt}(1-{(\mathrm{Re}+\mathrm{Alt})}^2\\ *\mathrm{si}{n}^2\left(d\right))/{\mathrm{Re}}^2)-(\mathrm{Re}+\mathrm{Alt})*\cos \left(d\right)\left]\right\}\\ z^{\′\′}=a\sin [(\mathrm{Re}+\mathrm{Alt})*\sin (d+\mathrm{dd})/\mathrm{Re}]\\ \mathrm{ds}=z^{\′\′}-s-(d+\mathrm{dd})\\ \mathrm{dS}=\mathrm{Re}*\mathrm{ds}\\ \mathrm{terrain}\_\mathrm{offset}=\mathrm{dS}/\mathrm{pixel}\_\mathrm{size}\end{array}\right.$ (式IV)

式IV中，S表示象元点距离星下基线的距离，pixel_size代表象元大小，off_nadir_pixel代表待校正象元距离星下的象元个数，s代表地球中心角度，Re代表地球参考椭球体在影像中心位置的地球半径，Alt为卫星的椭球体上空高度，terrain_offset为地形引起的位移。

根据以下公式拟合获取的多项式：

$\left\{\begin{array}{c}{x}^{\′}=a_0+a_{1}X+a_{2}Y+a_3{X}^2+a_4\mathrm{XY}+a_5{Y}^2+...\\ y^{\′}=b_0+b_{1}X+b_{2}Y+b_3{X}^2+b_4\mathrm{XY}+b_5{Y}^2+...\end{array}\right.$ (式V)

式V中，(x′y′)为图像坐标，(X，Y)为参考坐标，a_i，j，b_i，j为多项式的系数。

(4)精度评价：通过光谱匹配，在基准影像和输出的正射校正结果之间搜索同名地物点作为验证点，将验证点按照影像中心分为左上、右上、左下、右下四个象限，分别判断每个象限内同名地物点配准误差，仅当四个区域的验证点比例均大于60％时，报告当前影像通过验证，否则增加配准多项式阶数配准后返回步骤(3)重新进行几何纠正与正射校正；迭代次数小于阈值，程序默认为1，即迭代1次。

(5)影像输出：当达到程序收敛条件后，以纠正影像的原始影像范围输出经过几何纠正与正射校正的影响。

为具体说明本发明所述方法对环境减灾卫星的校正方式，下面选取9景不同时段的环境减灾卫星多光谱CCD影像进行实例演示。选用轨道号为17-76，

时间介于2009年至2010年的共计9景不同时相的HJ影像作为待校正影像(表1)。开展正射校正试验。

表1研究采用的HJ及基准影像轨道号与获取日期

(1)基准影像的选择与拼接

图9A为本实例选取的2009年8月31日HJ1BCCD1相机获取的影像覆盖情况。该影像覆盖范围包括横断山区及青藏高原东部。海拔高度在455m至6225m之间，15°坡度以上面积占全部影像覆盖面积的75％，是我国典型高山峡谷区。图9B为用于光谱匹配的的TM基准影像，采用近红外波段拼接。

采用TM/ETM影像作为基准影像开展HJ影像的正射校正。由于HJ影像幅宽较TM大，且同一轨道号下重访时存在一定的轨道漂移。为了使基准影像与待矫正影像匹配，本实例在马里兰大学GLCF(http：//glcf.umiacs.umd.edu/data/)上下载对应于2010年HJ轨道号的经过几何精校正与正射校正的全部landsatTM/ETM影像，最终选取6景云量较少覆盖范围与HJ对应的TM作为基准影像，其影像平均地理位置误差小于50m。TM影像基准信息如表1示，拼接效果见图9B。

正射校正基准DEM选用马里兰大学GLCF填充后的SRTM填充后的版本。选择WRS-2分幅方式以便与基准TM影像对应。

(2)同名地物点的自动选取

由图4可以看出，象元点随着距离星下基线位置的增加，同一海拔高度下的HJ影像几何定位误差逐渐增大。在星下距离相同的条件下，海拔越高几何定位误差越大。HJ影像的标称幅宽为360公里，由图中可以看出，影像远离星下基线的边缘地区，海拔在2000m的区域地理位置误差将超过1000m，约33个象元，随着海拔高度的增加，地理误差逐渐增大。由此可见，地形对坐标误差的影响是十分严重的。

本实例选取2009年8月31日的HJ1B-CCD1影像对控制点选取情况进行分析。通过光谱匹配算法共搜索出1254个同名地物点。图10A、图10B和图10C给出了程序搜索得出的同名地物点分布情况。可以看出，程序搜索得到的同名地物点空间分布范围较为均匀。

(3)同名地物点的筛选、几何精纠正与正射校正

根据前文分析，CCD1星下基线位于影像左侧边缘，本实例通过设定其侧视角度为15度，结合头文件信息中的景中心位置，计算得出的其星下基线方程为i＝-4.6368j+18231。

为了反映星下基线对影像校正误差的影响规律，本实例将星下基线由影像中心向星下基线方向以基线方程截距增加1个象元的方式进行移动，观察每次控制点筛选情况。每次星下基线移动后，重新计算待校正影像星下基线投射于基准影像上的位置，计算基准影像上控制点的地形坐标偏差，输出控制点的筛选情况。图像坐标系下控制点个数随截距变化的规律如下。图11中，截距为48626时为影像中心象元位置，虚线为景中心位置。由图中可以看出，随着HJ影像假设星下基线由影像中心向卫星边缘移动，有效控制点的个数逐渐增加。至影像真实星下基线时，有效控制点个数达到最大。经筛选后程序最终选择的有效控制点个数为105个，其拟合的多项式RMS为0.75个象元。

在程序校正时，利用控制点拟合的多项式RMS阈值设定为0.75个象元，但其仅能反应控制点的误差特征，难以反应影像整体的误差分布情况。本方法在正射校正完成后，通过光谱匹配，在基准影像和校正结果之间自动搜索验证点进行精度验证。

(4)精度验证

2009年8月31日HJB-CCD1影像共搜索得到了2272个验证点。程序根据影像中心，将验证点分为左上、右上、左下、右下四个区域进行验证。当各个区域验证点误差小于阈值条件的象元个数小于2个象元的60％时，影像通过误差验证。图12A和图12B示出了最终的校正结果。图12A为a、b、c区域所在大区域的总体图。图12B中的a区域位于影像星下基线附近，由对比影像可以看出，HJ校正前影像与基准影像一致，山体变形视觉效果不明显，校正后变化不明显。图12B中的b区域位于影像右侧中心左右位置，远离星下点240公里，可以看出，影像中山体已发生明显位移，西北坡区域被明显拉长，东南坡被压缩。纠正后影像与基准影像一致，地形位移得到了很好的改正；图12B中的c位于影像远离星下点的边缘区域348公里处，由图中可以看出，纠正前山体与基准影像相比变形较图12b更为明显，西北坡明显被拉长，东北坡被压缩，影像纠正后山体变形得到了明显改善，影像整体视觉效果与基准影像相同，整体符合实际情况。

采用批处理的方法，全部9景影像的纠正结果的验证情况如表2所示。

表2各影像验证点分布

所有验证点的误差直方图如图13示。可以看出，本方法在本实例中总体纠正精度为37m，小于2个象元，满足山区2个象元的制图需求。

(5)算法效率

本方案测试的硬件平台为CPU：Intel(R)Core(TM)2DuoCPUE84003.00GHZ，RAM：2GB，软件平台Windows732bit。测试结果为，单景影像的处理时间为20.3分钟，而传统人工控制点选择与筛选纠正1景影像往往需要6-8个小时左右，而且不同操作人员的精度有很大差异。本方案能够大大提高工作效率，可快速完成多幅卫星影像的处理工作，减少人力、物力，在海量环境减灾卫星的数据处理中有十分重要的意义。

Claims (8)

1.一种山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：包括以下步骤：

(1)基准影像的选择与拼接：a、基准影像的选择：基准影像满足以下条件：整体云覆盖量小于5％；不同平台影像及邻近影像接边处配准误差小于0.5个象元；b、基准影像的拼接：拼接基准影像的近红外波段；

(2)同名地物点的自动选择：a、在基准影像和纠正影像上分别布设核窗口和搜索窗口，所述纠正影像为待纠正和校正的山区多光谱遥感卫星影像，下述亦然；b、将基准影像上核窗口内的关联点投影至纠正影像，然后以关联点为中心，计算纠正影像与基准影像的核窗口内的相关系数；c、在待纠正影像的搜索窗口内移动该核窗口，并计算相应核窗口内的相关系数；d、当搜索窗口内的核窗口全部计算后，若存在关联点核窗口的相关系数大于设定的阈值条件，且搜索窗口满足条件的关联点小于关联点数量上限阈值条件时，将相关系数最大的点认为是该窗口内的同名地物点；否则该搜索窗口内未找到同名地物点；

(3)同名地物点的筛选、纠正影像的几何纠正与正射校正：a、同名地物点的筛选：对纠正影像进行边缘检测，估算出星下基线斜率，然后根据从纠正影像头文件获取的影像中心点坐标与卫星的设计高度，估算出纠正影像中心点坐标对应的星下点坐标，进而根据点斜式线性方程估算出纠正影像星下基线方程截距；在获取星下基线方程后，将该基线方程投影至基准影像，结合基准影像上同名地物点对应数字高程模型DEM的高程信息与卫星高度，计算每个基准影像上同名地物点偏离星下点的地形偏移，叠加上每个同名地物点的对应的偏移以恢复成为与纠正影像一致的观测条件；b、纠正影像的几何纠正：所有同名地物点拟合获取的多项式均方根误差不能超过0.75个象元，否则去除拟合误差最大的点重新拟合多项式系数，并再次计算拟合多项式的均方根误差，同名地物点的预测几何定位误差大于1个象元时，去除该点，当所有同名地物点均满足纠正要求时，根据满足筛选条件控制点拟合获取多项式系数，对纠正影像的每个象元位置进行重新计算；c、纠正影像的正射校正：结合纠正影像的星下基线方程、卫星高度及数字高程模型DEM的数据，根据每个象元经过几何纠正后新的坐标信息，计算每个象元距星下基线的垂直距离，再计算每个象元由于偏离星下观测而导致的坐标偏差，去除由于偏离星下观测导致的偏差而完成正射校正；d、计算完毕后，对原始影像的光谱特征进行重采样；

(4)精度评价：通过光谱匹配，在基准影像和输出的正射校正结果之间搜索同名地物点作为验证点，将验证点按照影像中心分为左上、右上、左下、右下四个象限，分别判断每个象限内同名地物点配准误差，仅当四个区域的验证点比例均大于60％时，报告当前影像通过验证，否则增加配准多项式阶数配准后返回步骤(3)重新进行几何纠正与正射校正；

(5)影像输出：当达到程序收敛条件后，以纠正影像的原始影像范围输出经过几何纠正与正射校正的影响；

所述步骤(3)中，“估算纠正影像中心点坐标对应的星下点坐标”、“对纠正影像的每个象元位置进行重新计算”和“计算每个象元由于偏离星下观测而导致的坐标偏差”均采用以下公式：

$\left\{\begin{array}{c}\mathrm{dis}=\mathrm{Alt}*\tan \mathrm{\θ}\\ \mathrm{\Δx}=\mathrm{dis}*\sin \mathrm{\α}\\ \mathrm{\Δy}=\mathrm{dis}*\cos \mathrm{\α}\\ \mathrm{wcx}=\mathrm{cx}-\mathrm{\Δx}\\ \mathrm{wcy}=\mathrm{cy}-\mathrm{\Δy}\end{array}\right.$ (式III)

式III中，dis为景中心位置距星下基线的距离，Alt为卫星高度，θ为卫星侧视角度，α为星下基线倾斜角度，Δx和Δy分别为景中心位置坐标偏移量，cx和cy为影像景中心位置，该中心位置通过与影像一起分发的XML头文件信息scenecenterlong和scenecenterlat获取，(wcx，wcy)为纠正影像中心点坐标对应的星下点坐标。

2.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(2)中，将基准影像上核窗口内的关联点按下述公式的坐标投影至纠正影像：

$\left\{\begin{array}{c}\mathrm{wx}=\frac{\mathrm{bx}*{\mathrm{res}}_{\mathrm{base}}+({\mathrm{ulx}}_{\mathrm{base}}-{\mathrm{ulx}}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\\ \mathrm{wy}=\frac{\mathrm{by}*{\mathrm{res}}_{\mathrm{base}}+({\mathrm{uly}}_{\mathrm{base}}-{\mathrm{uly}}_{\mathrm{warp}})}{{\mathrm{res}}_{\mathrm{warp}}}\end{array}\right.$ (式I)

式I中，(bx，by)为基准影像的坐标，(wx，wy)为纠正影像的估计坐标，(ulx_base，uly_base)为基准影像左上角象元坐标，(ulx_warp，uly_warp)为纠正影像左上角坐标，res_base和res_warp为基准影像和纠正影像的空间分辨率。

3.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(2)中，按下述公式计算纠正影像与基准影像的核窗口内的相关系数：

$r_{\mathrm{mn}}=\frac{\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}(A_{\mathrm{ij}}-\stackrel{\‾}{A})(B_{\mathrm{ij}}-\stackrel{\‾}{B})}{\sqrt{\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(A_{\mathrm{ij}}-\stackrel{\‾}{A})}^2\right)\left(\underset{i=m-\mathrm{cs}/2}{\overset{m+\mathrm{cs}/2}{\mathrm{\Σ}}}\underset{j=n-\mathrm{cs}/2}{\overset{n+\mathrm{cs}/2}{\mathrm{\Σ}}}{(B_{\mathrm{ij}}-\stackrel{\‾}{B})}^2\right))}}$ (式II)

式II中，(m，n)为核窗口中心坐标，A_ij为基准影像的象元光谱值；B_ij为纠正影像的象元光谱值，为基准影像核窗口内光谱平均值；为纠正影像核窗口内光谱平均值，cs代表核窗口大小。

4.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(2)中，设定的核窗口大小为7×7个象元，搜索窗口大小为50×50个象元。

5.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(2)中，阈值为0.75，设定的搜索窗口内满足阈值条件的最大个数为3个。

6.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(3)中，采用以下公式计算每个同名地物点偏离星下点的地形偏移：

$\left\{\begin{array}{c}S=\mathrm{pixel}\_\mathrm{size}*\mathrm{off}\_\mathrm{nadir}\_\mathrm{pixel}\\ s=S/\mathrm{Re}\\ \mathrm{LOS}=\mathrm{sqrt}({\mathrm{Re}}^2+{(\mathrm{Re}+\mathrm{Alt})}^2-2*\mathrm{Re}(\mathrm{Re}+\mathrm{Alt})*\cos \left(s\right))\\ d=a\sin (\mathrm{Re}*\sin \left(s\right)/\mathrm{LOS})\\ a\tan \{(\mathrm{Re}+\mathrm{Alt})*\sin \left(d\right)*(1-(\mathrm{Re}+h)/\mathrm{Re})/[(\mathrm{Re}+h)*\mathrm{sqrt}(1-{(\mathrm{Re}+\mathrm{Alt})}^2\\ *{\sin }^2\left(d\right))/{\mathrm{Re}}^2)-(\mathrm{Re}+\mathrm{Alt})*\cos \left(d\right)\left]\right\}\\ z^{\′\′}=a\sin [(\mathrm{Re}+\mathrm{Alt})*\sin (d+\mathrm{dd})/\mathrm{Re}]\\ \mathrm{ds}=z^{\′\′}-s-(d+\mathrm{dd})\\ \mathrm{dS}=\mathrm{Re}*\mathrm{ds}\\ \mathrm{terrain}\_\mathrm{offset}=\mathrm{dS}/\mathrm{pixel}\_\mathrm{size}\end{array}\right.$ (式IV)

式IV中，S表示象元点距离星下基线的距离，pixel_size代表象元大小，off_nadir_pixel代表待校正象元距离星下的象元个数，s代表地球中心角度，Re代表地球参考椭球体在影像中心位置的地球半径，Alt为卫星的椭球体上空高度，terrain_offset为地形引起的位移。

7.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(3)中，根据以下公式拟合获取的多项式：

$\left\{\begin{array}{c}{x}^{\′}=a_0+a_{1}X+a_{2}Y+a_3{X}^2+a_4\mathrm{XY}+a_5{Y}^2+...\\ y^{\′}=b_0+b_{1}X+b_{2}Y+b_3{X}^3+b_4\mathrm{XY}+b_5{Y}^2+...\end{array}\right.$ (式V)

式V中，(x′，y′)为图像坐标，(X，Y)为参考坐标，a_i，j，b_i，j为多项式的系数。

8.根据权利要求1所述的山区多光谱遥感卫星影像的自动几何纠正与正射校正方法，其特征在于：所述步骤(3)中，重采样方法为最邻近法、二次内插法或三次卷积法。