CN103218789A - 一种基于几何形变分步解算的自动化正射纠正方法 - Google Patents

Abstract

本发明提供一种自动化的中分辨率卫星影像的正射纠正方法，特别是对地形起伏明显的地区，能够快速精确地实现影像的正射纠正过程。该方法是建立在对遥感影像成像过程中几何形变产生原因的分析基础上，将不同的几何形变类型分离，对一致性的成像系统造成的畸变和地形起伏造成的畸变进行分步解算。首先通过模拟成像的手段模拟地形起伏造成的不规律形变，将生成的模拟影像作为自动配准参考影像。随后通过自动配准方法，用多项式模型纠正具有一致性的成像系统畸变。最后解算目标影像配准后存在的地形偏差。该方法可适用于不同传感器拍摄的影像，针对不同的传感器，地形解算部分采用查找表记录的其地形形变，大大提高了整个纠正流程的计算速度，自动高效地实现了影像的正射纠正。

Description

一种基于几何形变分步解算的自动化正射纠正方法

技术领域

本发明涉及遥感影像处理技术，具体地说，涉及对遥感影像中几何形变的校正技术以及影像的自动配准技术。本发明适用于推扫式中分辨率卫星影像的正射校正过程。

背景技术

遥感影像的几何处理是遥感影像预处理的一个关键步骤，是多源影像融合、遥感影像处理及信息提取等后续应用正确进行的前提。正射校正可以选择的方法主要包括严格物理模型和通用经验模型两种。严格物理模型以传感器构像方程为依据，为了获得较高的精度需要已知传感器的轨道参数和姿态参数等等，其实现过程相对复杂；经验模型应用灵活，只要有足够数量的控制点就可以获得正射影像，但是其精度往往受到地形和控制点的限制，在地形起伏较大的地区精度较差。因此，需要一种快速高效的校正方法对不同地区的形变进行校正。相应的正射校正方法的参考文献包括D.

R.Richter，Geo-atmospheric processing ofairborne imaging spectrometry data，International Journal of Remote Sensing，200223(13)2609-2630、R.B.Huseby，O.M.Halck，R.Solberg，A model-based approach for geometricalcorrection of optical satellite images，International Journal of Remote Sensing，200526(15)3205-3223.

卫星影像在成像的过程中，受到透视投影、摄影轴倾斜、大气折光、地球曲率及地形起伏等诸多因素影响，致使影像中各像点产生不同程度的几何变形而失真。对遥感影像的几何形变进行分类可更好地针对不同误差进行纠正。其中，从形变来源的角度出发，可分为影像获取系统引起的形变和被观测物引起的形变。但目前常用的几何校正模型中及多源影像配准技术中很少考虑形变的类型，导致校正结果在整体平差效果良好但局部误差仍无法有效拟合的情况。相应的分析遥感影像几何形变的参考文献包括T.TOUTIN，Review article：Geometricprocessing of remote sensing images：models，algorithms and methods，International Journal ofRemote Sensing，200425(10)1893-1924、Glasbey C.A.，Mardia K.V.A review of image-warpingmethods[J].Journal of Applied Statistics，1998，25(2)，155-171.

自动配准技术是遥感影像几何纠正自动化过程中常用的手段，多源影像配准也是遥感数据生产及应用的关键步骤。但由于整体变形的不一致性，特征点的自动选取效果往往不好，因此需要一种合理的自动选取特征点的方法。而近几年来出现的SIFT算子、SURF算子虽然提高了特征点提取的精度，但其实现算法复杂度高，计算时间较长，且无法从根本上解决配准中模型拟合上的局部误差。因此需要根据遥感数据的特点确定一种合理高效的配准方法。相应的影像配准的参考文献包括Barbara Zitova，Jan Flusser，Image registration methods：A survey，Image and Vision Computing，2003(21)，977-1000.

在卫星影像正射校正流程方面，目前可见专利/文献中相互独立的研究较多，但从整体流程的角度提出的正射校正的方案不多，且往往在自动化程度及局部校正精度上有所欠缺。

发明内容

本发明的目的是提供一种针对中分辨率卫星的自动化正射纠正流程。该流程通过对成像过程中的几何形变分离解算，能够有效地纠正不同地形遥感影像的局部形变，并通过自动化的方法实现，能够应用于目前亟需的国产卫星数据生产过程。

本发明的思路为：通过对遥感影像成像过程中的几何形变的分析，将形变类型分为成像系统的整体变形和地形起伏在全景投影下的局部变形，并通过一定的手段将两种形变分阶段校正。首先通过由基准影像模拟成像的方法消除成像过程中地形起伏的投影变形；然后将模拟成像的影像与待纠正影像选取合适的流程进行自动配准以消除成像系统的整体变形。对于同类型传感器的地形起伏形变，通过建立查找表的方式获得效率上的提升。

本发明的技术方案提供了一种通过将不同几何形变分离解算的自动化正射校正流程，所需辅助数据包括能够覆盖目标影像范围的基准影像及DEM影像，具体包含以下实施步骤：

1)数据准备过程，系统校正后的影像通常存在一定像素的粗偏轨，故首先将目标影像与基准影像进行FFT配准，粗纠正目标影像中的平移形变，根据目标影像的地理坐标准备相应范围内的基准影像和DEM影像，并重采样至相同分辨率；

2)生成查找表，查找表记录对应传感器在不同高程下地形的投影偏差，每一种传感器对应一个查找表；

3)基准影像结合DEM提供的高程信息从步骤2)生成的查找表中获取相应的地形偏差，进行成像模拟，模拟结果保存为模拟影像，模拟影像相当于在真实坐标下添加了局部投影形变；

4)在模拟影像上提取控制点，以相关系数作为目标影像与模拟影像的相似性测度，使目标影像与模拟影像相匹配；

5)匹配后的影像根据步骤3)进行地形偏差解算，消除步骤3)中添加的地形起伏偏差，生成最终正射产品。

上述步骤的特征在于：

步骤1)中目的在于粗纠正以确定目标影像大概对应的地理范围，避免平移形变较大的情况，从而使后续步骤4)的配准过程在较小的搜索窗口中就能搜索到同名点。FFT配准方法只需针对平移形变，不需要变换至极坐标空间计算旋转和缩放形变，另外该步骤的目的为粗纠正，因此可以只裁剪出小块区域进行配准操作，可以大大减少实际运算量。

步骤2)中计算投影偏差的只跟成像时的几何位置相关，因此记录对应影像上所有点在可能高程范围的偏差即可满足该类传感器所有地形偏差的解算，提高了流程效率。

步骤3)与步骤5)都是根据步骤2)中生成的查找表查找偏差，避免了不同景影像纠正时的重复解算，不同的是步骤3)从正向模拟，由现有坐标解算出模拟影像坐标；而步骤5)采用反向解算，即根据待求解的正射影像坐标反算出重采样点的位置再进行重采样。

步骤4)中提取特征点可采用Moravec算子、Frostner算子或Harris算子等经典算子，由于步骤1)和步骤3)中剔除了平移畸变和地形投影的偏差，因此在该配准过程中采用以上算子就可以达到快速准确的效果。

本发明与现有技术相比具有以下特点：从影像形变的效果出发，将不一致的地形形变与整体一致的成像系统形变分解，良好地解决了用模拟模型不适用地形起伏较大地区的问题；并在此基础上简化了影像自动配准的过程，使得正射校正自动化实现的精度大大提高；另外，以查找表的方式减少了计算复杂度，且查找表与传感器一一对应使得该流程能够适用于不同传感器拍摄的影像。

附图说明

图1是正射纠正方法流程图。

图2是FFT配准算法流程图。

图3是成像过程形变示意图。

图4是地形偏差几何关系示意图，其中h为地物高程，Re为地球半径，Alt为传感器高度，LOS和LOS’分别为正确视距线和偏差视距线，dS为成像偏差。

图5是模拟/解算地形偏差流程图，其中模拟地形偏差由上至下正向处理，解算地形偏差由下至上逆向解算。

图6是自动配准流程图。

图7是特征点提取效果图。

图8是自动配准局部效果图，其中左图为配准前对比图，右图为配准后对比图。

图9是正射纠正整体效果图，其中上图为校正前对比图，下图为校正后对比图。

具体实施方式

图1是整个正射纠正流程实现的流程图：该流程适用的影像为中分辨率系统纠正级卫星影像，实际实现过程中以近几年获取的CBERS02B及HJ星数据作为实验数据。所需要的辅助数据包括覆盖目标数据地理范围的基准影像和DEM影像，基准影像可以用相近分辨率下的正射影像，在本流程中整个流程分为FFT配准粗纠正、模拟成像、自动配准及地形偏差解算四个主要步骤，具体实现过程如下：

首先应对辅助数据进行一定的预处理，需裁剪出与目标影像相同的地理范围的部分并重投影至相同投影坐标系下。但由于卫星成像过程中传感器外方位元素的变化，实际成像区域与设定成像区域往往存在一定距离的偏差，在国产卫星影像上这类粗偏差尤为明显。因此为保证步骤5)自动配准能够正确进行，需要在数据准备时就剔除此类偏差。图2是FFT配准算法的流程图，该步骤利用影像在傅里叶空间下的特征关系，可以快速确定基准影像与目标影像的坐标平移偏差。

FFT配准实现算法如下：

1)将基准影像与目标影像分别进行FFT，得到频率域影像F₁(δ，θ)，F₂(δ，θ)；

2)计算频率域下的比例因子R，计算公式如下

$R=\frac{F_1(\mathrm{\δ},\mathrm{\θ})*\mathrm{conj}\left(F_2(\mathrm{\δ},\mathrm{\θ})\right)}{\mathrm{abs}\left(F_1(\mathrm{\δ},\mathrm{\θ})\right)*\mathrm{abs}\left(F_2(\mathrm{\δ},\mathrm{\θ})\right)}$

其中conj表示复共轭，abs表示复数的绝对值，δ和θ表示频率影像上的坐标；

3)统计比例因子最大值对应的δ和θ坐标，即在x、y方向上的偏移值。

由于FFT变换计算时空间复杂度较高，在实际的实现过程中可只裁剪1024*1024大小的目标影像进行计算。

图3是成像过程形变示意图，是本发明的理论基础。从地理学的角度而言，几何校正的目的是确定影像坐标点与真实地面点的坐标一一对应的过程。如图3所示，遥感卫星成像形变包括：

1)传感器外方位元素变化，传感介质不均匀，传感器机械误差等一系列来源于成像系统的畸变，该变化公式

P(x，y)＝I(X，Y，Z)*F

其中I(X，Y， Z)表示真实地面亮度值，F表示系统畸变函数；

2)在有地形起伏情况下正射投影与全景投影间的投影形变。该类形变并非成像过程中的错误，而是不同投影方式之间的差异。该变化公式

O(x，y)＝P(x，y)*G

其中P(x，y)表示全景投影影像，O(x，y)表示正射投影影像。

其中，第一类形变为成像系统形变，在影像上体现出整体性；第二类形变与地面点的高程相关，在影像上不一致，影像两侧的变形较为明显。本发明根据成像中畸变的产生顺序，首先通过成像模拟来消除目标影像和基准影像之间的第二类畸变，再通过多项式模型拟合消除第一类畸变。

图4是地形偏差几何关系示意图，地形偏差的解算参考该示意图中的几何关系，计算地形偏差dS。计算公式如下：

LOS＝sqrt(Re²+(Re+Alt)²-2*Re*(Re+Alt)*cos(s))

d＝asin(Re*sin(s)/LOS)

$\mathrm{dd}=a\tan \left(\frac{(\mathrm{Re}+\mathrm{Alt})*\sin \left(d\right)*(1-(\mathrm{Re}+h))/\mathrm{Re})}{(\mathrm{Re}+h)*\sqrt{1-{(\mathrm{Re}+\mathrm{Alt})}^2*{\sin }^{2}d/{\mathrm{Re}}^2}-(\mathrm{Re}+\mathrm{Alt})*\cos d}\right)$

$z^{\′\′}={\sin }^{-1}\left(\frac{(\mathrm{Re}+\mathrm{Alt})*\sin (d+\mathrm{dd})}{\mathrm{Re}}\right)$

ds＝Re*(z″-s-d-dd)

按照上述公式，计算传感器幅宽范围内所有采样点的地形形变，实际实现时采样点的高度范围设定为[-300，9000]，可以满足地球上所有高程。对于线阵推扫式传感器而言，地形偏差只存在于垂直运行轨道的方向。因此，查找表的记录方式为一个二维数组，行列分别表示与中心轨道的距离(像素距离，即实际距离/分辨率)和该位置处的高程。这样的记录方式能够满足查找该类传感器的所有高程下的起伏情况，在解算地形畸变时可减少逐象元计算的过程。

图5为模拟/解算地形偏差流程图，模拟与解算地形偏差为互逆的处理，因此操作流程相同，只是输入输出数据及解算过程相反。以解算地形偏差流程图为例，实现步骤如下：

1)逐像素依次将待生成的正射影像上的坐标由像平面坐标系(x，y)转换为像轨道坐标系坐标(x₀，y₀)，此处像轨道坐标系定义为像平面上的二维坐标系，以像平面上轨道线为Y轴，轨道垂线为X轴，通过计算影像有效范围的中垂线确定；

2)由x₀和DEM影像上对应位置的高程值z查找成像偏差Δx，即在模拟影像上对应的点坐标x_p＝x₀+Δx；

3)将(x₀，y₀)还原到像平面坐标系上，得到目标影像上重采样点(x’，y’)。

根据之前成像过程中的形变分析，模拟影像与目标影像之间的形变只剩下整体性系统畸变。图6是通过自动配准方法消除系统畸变的过程，具体算法如下：

1)根据设定的特征点总个数对基准影像进行均匀分块。特征点的总个数由影像的实际大小决定，确定原则为尽量均匀覆盖整幅影像且不宜过多，实验中对于7000*7000的影像特征点在设定为400个左右；

2)在每一块区域中心划分一个256*256大小的窗口，用Frostner算子计算该窗口内每个象元特征值，设定阈值Tq＝0.7，Tw＝1.05*Ew。其中Ew表示特征值w的平均值。特征值超过阈值的点选取为待选点，从待选点中选取特征值w最大的点为该窗口的最佳特征点；

3)每个特征点在目标影像上对应坐标处设定搜索窗口寻找匹配点，由于之前进行了平移畸变的粗纠正，可以保证在不大的搜索窗口就能搜索到匹配点。设定搜索窗口为49*49，移动窗口为15*15，采用两幅影像移动窗口之间的相关系数作为相似性测度，阈值设定为0.85。大于该阈值时认为两个窗口匹配良好，在所有大于阈值的点中选取最大值作为匹配点，若不存在满足要求的匹配点，则认为该区域没有匹配点；

4)将选取的特征点对用多项式模型拟合，使用迭代最小二乘法平差，最大平差阈值设定为1个像素；

5)按照求解的模型对目标影像进行几何变换与重采样。

图7为特征点提取效果图，其中背景影像(基准模拟影像)为拼接裁剪的Landsat正射影像，前景目标影像为CBERS02B影像。红色叉状符号为基准模拟影像上提取的特征点，黄色十字符号为目标影像上的配准点。从图中可以看出，匹配点对的距离均在较近范围内，且均匀分布在影像上。图8为两幅影像配准结果的局部放大效果图，该自动配准整体平差精度在1个象元内，多项式模型拟合效果良好。

图9为正射纠正流程整体效果图，选定地理范围为新疆天山山脉地区，平均海拔4000米。图中效果显示，本发明在地形起伏较为明显的地区也能取得良好的纠正效果。

本发明的实例在PC平台上实现，通过自动化的方法实现了正射纠正流程。整个流程在7000*7000的影像上耗时不到300秒，并达到了较为理想的纠正效果。本发明中所提及的方法可广泛用于中分辨率卫星影像的正射校正处理，如国产卫星数据(HJ星等)生产等应用。

Claims (3)

1.本发明提供了一种基于几何形变分步解算的自动化正射纠正方法，具体包含以下步骤：

步骤1：数据准备过程，首先将目标影像与基准影像进行FFT配准，粗纠正目标影像中的平移形变，根据目标影像的地理坐标准备相应范围内的基准影像和DEM影像，并重采样至相同分辨率；

步骤2：生成查找表，查找表用于记录对应地形起伏的形变，地形形变由像素距星下点距离及像素点高程同时确定，每一种传感器对应一个查找表；

步骤3：基准影像的二维坐标结合DEM提供的高程信息，从步骤2中生成的查找表中获取相应的地形偏差，进行成像模拟，模拟结果保存为模拟影像，模拟影像相当于在真实坐标下添加了局部投影形变；

步骤4：在模拟影像上提取控制点，以相关系数作为目标影像与模拟影像的相似性测度，在目标影像选取匹配点，将匹配点对作为控制点，用迭代最小二乘法解算多项式模型，将目标影像配准到模拟影像上；

步骤5：匹配后的影像根据步骤3进行地形偏差解算，将之前添加的局部地形形变消除，生成最终正射产品。

2.根据权利要求1所述的模拟成像的方法，其特征在于该方法基于几何形变分步解算的原理，在自动配准前消除不规律的地形形变，从而使得后一步的自动配准达到良好的效果。

3.根据权利要求1中所述的查找表生成方法，其特征在于将传感器参数与正射纠正计算流程相独立，使得该流程的计算过程可适用于不同的传感器影像，计算过程中只需根据不同的传感器类型查找不同的查找表；另外，查找表的存在可大大减少整个流程的计算时间。

Images