一种光学遥感卫星正射影像拉花区域自动检测方法

摘要

本发明公开了一种光学遥感卫星正射影像拉花区域自动检测方法，首先加载原始卫片、RPC参数和数字高程模型DEM，并迭代计算正射校正后正射影像的大小和范围；然后利用反解法逐像素反算每个像素对应原始影像上像素的坐标；最后逐像素统计并判断是否为拉花像素点，形成拉花区域二值标记图像并进行矢量化，得到拉花区域矢量文件，达到光学遥感卫星正射影像上拉花区域自动检测的目的。其显著效果是：实现了光学遥感卫星正射影像中拉花区域的自动检测，解决了传统人工目视辨别拉花区域费时费力，效率低下的问题，大大提高了光学遥感卫星正射影像的质量检查效率。

权利要求

1.一种光学遥感卫星正射影像拉花区域自动检测方法，其特征在于按照以下步骤进行：步骤1：加载原始卫片、RPC模型参数和数字高程模型DEM，并迭代计算正射校正后正射影像的大小和范围；

步骤2：利用反解法逐像素反算每个像素对应原始影像上像素的坐标；

步骤3：逐像素统计并判断是否为拉花像素点，形成拉花区域二值标记图像并进行矢量化，得到拉花区域矢量文件。

2.根据权利要求1所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：所述步骤1的具体处理步骤如下：

步骤1.1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；

步骤1.2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；

步骤1.3：由地面点坐标根据RPC模型反算出该地面点在原始影像上的位置；

步骤1.4：基于数字高程模型DEM迭代求解，计算正射校正后影像的大小。

3.根据权利要求2所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：步骤1.2中所述地面点高程值的计算步骤为：步骤1.2.1：根据目标点的坐标和数字高程模型DEM，计算得到目标点在数字高程模型DEM格网中的行列数；

步骤1.2.2：根据得到的行列数计算数字高程模型DEM格网中最邻近的四个高程点，并内插出目标点的高程值。

4.根据权利要求2所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：所述步骤1.3中反算地面点在原始影像上位置的步骤为：步骤1.3.1：将地面点坐标按照RPC模型参数中的空间坐标标准化参数计算得到标准化之后的空间坐标；

步骤1.3.2：根据RPC模型计算该点在原始影像上的标准化坐标；

步骤1.3.3：根据像点坐标标准化参数对标准化坐标进行处理，得到原始影像上的像点坐标。

5.根据权利要求2所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：所述步骤1.4中正射校正后影像的大小计算步骤为：步骤1.4.1：根据步骤1.1计算得到原始影像左上角像点的初始地理坐标，根据步骤1.2利用双线性内插方法在DEM上内插得到初始地理坐标的高程值，根据步骤1.3得到地面点在原始影像上的像点坐标；

步骤1.4.2：计算像点坐标对应的初始地理坐标，内插数字高程模型DEM并迭代求解，得到原始影像左上角像点对应的地面点坐标；

步骤1.4.3：按照步骤1.4.1与步骤1.4.2迭代求解，直至前后两次计算中内插DEM的高程值的差值小于阈值，则最后一次计算得到的地面点坐标为左上角像素点对应的地面点坐标;步骤1.4.4：按照步骤1.4.1至步骤1.4.3分别计算原始影像四个角点对应地面点的坐标，并计算得到X、Y方向的极值；

步骤1.4.5：根据所得的X、Y方向的极值，计算得到正射校正后正射影像的行数和列数。

6.根据权利要求1所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：步骤2中所述反解法的具体步骤如下：

步骤2.1：根据当前像素在正射校正后影像中的行列坐标计算得到当前像素的地面点坐标；

步骤2.2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点的高程值；

步骤2.3：根据RPC模型反算出地面点坐标在原始影像上的像点坐标；

步骤2.4：将所有像素按照步骤2.1 2.3计算得到其对应的原始影像上的像点坐标。

~

7.根据权利要求1所述的光学遥感卫星正射影像拉花区域自动检测方法，其特征在于：步骤3中拉花区域矢量文件的获得步骤如下：

步骤3.1：以当前像素为中心，建立一个适当大小的方形窗口；

步骤3.2：遍历和统计方形窗口中所有像素对应原始影像上像点坐标和和当前中心像素对应原始影像上像点坐标的重叠次数；

步骤3.3：当重叠次数大于阈值时，在检测结果影像上将当前中心像素标记为拉花像素，否则标记为非拉花像素；

步骤3.4：逐像素统计每个像素对应原始影像像点坐标与周边像素对应原始影像像点坐标的关系，找出纠正后正射影像上的所有拉花像素与其构成的拉花区域二值标记图像；

步骤3.5：对拉花区域二值标记图像进行矢量化，得到拉花区域矢量文件。

说明书

一种光学遥感卫星正射影像拉花区域自动检测方法

技术领域

[0001] 本发明涉及到光学遥感卫星正射影像处理技术领域，具体地说，是一种光学遥感卫星正射影像拉花区域自动检测方法。

背景技术

[0002] 数字正射影像图(DOM)是利用光学遥感影像结合地形数据(数字高程模型DEM)根据相应的数学模型经过逐像素计算地面点与原始光学遥感影像像素的关系，并进行灰度重采样得到的既有正确的位置信息又有丰富的纹理信息的影像图。制作数字正射影像图的原始光学遥感影像根据获取方式可以分为两类：光学卫星遥感影像和低空航摄影像，两种光学影像的成像方式、数学模型和处理方法均不相同。由于光学卫星遥感影像具有获取速度快、影像框幅大、覆盖范围广的特点被广泛应用于国土资源监测、地理国情监测更新等国家重大战略项目中，在城市建设用地监测、城市违法建筑监测、城市规划管理、生态环境监测保护等领域也发挥着越来越重要的作用。

[0003] 光学卫星遥感影像的正射校正是制作数字正射影像图的关键一环，正射校正是利用光学卫星遥感影像及其对应的数学模型，如有理多项式模型(rational polynomial coefficients,RPC)和数字高程模型DEM数据共同消除原始卫星影像中的各种畸变(如投影差)而得到一幅既有地理坐标位置信息又有纹理信息的新影像的过程。地形起伏和卫星影像成像时光学镜头中心投影的成像方式不能保证所有地面点都能在卫星影像上成像，如坡势较陡的山坡可能会被山顶所遮挡。而在卫星影像正射校正时根据地面点位置计算原始影像上对应的像素，然后进行灰度重采样，在被遮挡的区域或是成像信息匮乏的区域重采样将会过于稠密或是重复采样，这样就会造成正射校正后的影像出现拉伸现象，如果拉伸过度就会造成纹理失真的现象，我们称之为“拉花”现象，纹理失真的区域我们称之为“拉花区域”。拉花区域造成的纹理失真直接影像着数字正射影像图的质量，特别是山地区域拉花现象特别严重，直接影像着数字正射影像图制作效率。

[0004] 目前，在国内外文献和专利中，还没有针对光学卫星正射影像拉花区域的自动检测方法，在目前的数字正射影像制作中，如果出现拉花区域，需要人工目视辨别和寻找，然后再通过修改DEM再纠正的方式进行处理。而光学卫星遥感影像框幅较大，如Worldview2卫星一景影像的覆盖范围可达280平方公里，依靠人工目视辨别查找拉花区域的方法将非常地耗时，同时也可能因为人工原因造成遗漏，因此光学卫星正射影像拉花区域的查找效率目前非常低下。

发明内容

[0005] 针对现有技术的不足，本发明的目的是提供一种光学遥感卫星正射影像拉花区域自动检测方法，利用原始光学卫星影像及其RPC模型参数和测区数字高程模型DEM，通过逐像素反解其对应原始卫星影像上的坐标，并进行统计比对，判定其是否为拉花像素，最后将生成的拉花区域二值图像进行矢量化得到拉花区域矢量数据，达到光学遥感卫星正射影像上拉花区域自动检测的目的。

[0006] 为达到上述目的，本发明采用的技术方案如下：

[0007] 一种光学遥感卫星正射影像拉花区域自动检测方法，其关键在于按照以下步骤进行：

[0008] 步骤1：加载原始卫片、RPC参数和数字高程模型DEM，并迭代计算正射校正后正射影像的大小和范围；

[0009] 步骤2：利用反解法逐像素反算每个像素对应原始影像上像素的坐标；

[0010] 步骤3：逐像素统计并判断是否为拉花像素点，形成拉花区域二值标记图像并进行矢量化，得到拉花区域矢量文件。

[0011] 进一步的，所述步骤1的具体处理步骤如下：

[0012] 步骤1.1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；

[0013] 步骤1.2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；

[0014] 步骤1.3：由地面点坐标根据RPC模型反算出该地面点在原始影像上的位置；

[0015] 步骤1.4：基于数字高程模型DEM迭代求解，计算正射校正后影像的大小。

[0016] 更进一步的，步骤1.2中所述地面点高程值的计算步骤为：

[0017] 步骤1.2.1：根据目标点的坐标和数字高程模型DEM，计算得到目标点在数字高程模型DEM格网中的行列数；

[0018] 步骤1.2.2：根据得到的行列数计算数字高程模型DEM格网中最邻近的四个高程点，并内插出目标点的高程值。

[0019] 再进一步的，所述步骤1.3中反算地面点在原始影像上位置的步骤为：

[0020] 步骤1.3.1：将地面点坐标按照RPC参数中的空间坐标标准化参数计算得到标准化之后的空间坐标；

[0021] 步骤1.3.2：根据RPC有理多项式模型计算该点在原始影像上的标准化坐标；

[0022] 步骤1.3.3：根据像点坐标标准化参数对标准化坐标进行处理，得到原始影像上的像点坐标。

[0023] 再进一步的，所述步骤1.4中正射校正后影像的大小计算步骤为：

[0024] 步骤1.4.1：根据步骤1.1计算得到原始影像左上角像点的初始地理坐标，根据步骤1.2利用双线性内插方法在DEM上内插得到初始地理坐标的高程值，根据步骤1.3得到地面点在原始影像上的像点坐标；

[0025] 步骤1.4.2：计算像点坐标对应的初始地理坐标，内插数字高程模型DEM并迭代求解，得到原始影像左上角像点对应的地面点坐标；

[0026] 步骤1.4.3：按照步骤1.4.1与步骤1.4.2迭代求解，直至前后两次计算中内插DEM的高程值的差值小于阈值，则最后一次计算得到的地面点坐标为左上角像素点对应的地面点坐标；

[0027] 步骤1.4.4：按照步骤1.4.1至步骤1.4.3分别计算原始影像四个角点对应地面点的坐标，并计算得到X、Y方向的极值；

[0028] 步骤1.4.5：根据所得的X、Y方向的极值，计算得到正射校正后正射影像的行数和列数。

[0029] 进一步的，步骤2中所述反解法的具体步骤如下：

[0030] 步骤2.1：根据当前像素在正射校正后影像中的行列坐标计算得到当前像素的地面点坐标；

[0031] 步骤2.2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点的高程值；

[0032] 步骤2.3：根据RPC模型反算出地面点坐标在原始影像上的像点坐标；

[0033] 步骤2.4：将所有像素按照步骤2.1～2.3计算得到其对应的原始影像上的像点坐标。

[0034] 进一步的，步骤3中拉花区域矢量文件的获得步骤如下：

[0035] 步骤3.1：以当前像素为中心，建立一个适当大小的方形窗口；

[0036] 步骤3.2：遍历和统计方形窗口中所有像素对应原始影像上像点坐标和当前中心像素对应原始影像上像点坐标的重叠次数；

[0037] 步骤3.3：当重叠次数大于阈值时，在检测结果影像上将当前中心像素标记为拉花像素，否则标记为非拉花像素；

[0038] 步骤3.4：逐像素统计每个像素对应原始影像上像点坐标与周边像素对应原始影像上像点坐标的关系，找出纠正后正射影像上的所有拉花像素与其构成的拉花区域二值标记图像；

[0039] 步骤3.5：对拉花区域二值标记图像进行矢量化，得到拉花区域矢量文件。

[0040] 本发明的显著效果是：本方案利用原始光学卫星影像及其RPC模型参数和测区数字高程模型DEM，通过逐像素反解其对应原始卫星影像上的坐标，并进行统计比对，判定其是否为拉花像素，最后将生成的拉花区域二值图像进行矢量化得到拉花区域矢量数据，达到光学遥感卫星正射影像上拉花区域自动检测的目的。实现了光学遥感卫星正射影像中拉花区域的自动检测，解决了传统人工目视辨别拉花区域费时费力，效率低下的问题，大大提高了光学遥感卫星正射影像的质量检查效率。

附图说明

[0041] 图1是本发明的方法流程图；

[0042] 图2是步骤1的具体处理流程图；

[0043] 图3是步骤2的具体处理流程图；

[0044] 图4是步骤3的具体处理流程；

[0045] 图5是卫星正射影像上检测出的局部拉花区域1；

[0046] 图6是卫星正射影像上检测出的局部拉花区域2。

具体实施方式

[0047] 下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。

[0048] 本实施例结合某多山地区一景分辨率为0.5米的WorldView2卫星影像，原始卫片大小为35283*15257，大小为1.33GB，正射校正后影像的大小为36085*14669，大小为1.03GB，覆盖地面范围约115Km2，对本发明方法进行详细说明。

[0049] 如图1所示，一种光学遥感卫星正射影像拉花区域自动检测方法，具体按照以下步骤进行：

[0050] 步骤1：加载原始卫片、RPC参数和数字高程模型DEM，并迭代计算正射校正后正射影像的大小和范围，具体处理步骤参见附图2：

[0051] 步骤1.1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；

[0052] 光学卫星遥感影像不仅有对应的RPC参数模型，同时也自带了原始影像上地理坐标与影像行列数之间变换的6个仿射变换参数，存储在数组adfGeoTransform[6]中，按照公式 可以通过仿射变换参数从一个像元的行列数得到其对应的地理坐标；

[0053] 式中，(row,col)是像素的行列坐标，(Xgeo，Ygeo)为该像素对应的地理坐标，其中，Xgeo为经度，Ygeo为维度,(adfGeoTransform[0],adfGeoTransform[3])为该影像左上角像素对应的地理坐标， 是行列坐标到地理坐标的四个仿射变换矩阵。

[0054] 步骤1.2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；

[0055] 数字高程模型DEM是将地形数据按照一定的间隔的格网进行规则存储，DEM中每个点的行列数与该点的地理坐标相关，数值则是该点的高程值Z，可以根据一个点的地理坐标(Xgeo，Ygeo)和DEM来内插该点对应的高程值。DEM内插方法分为三类：最邻近内插法，双线性内插法和三次卷积内插法，三种方法的计算量依次变大，本例中采用双线性内插法进行处理，具体步骤为：

[0056] 步骤1.2.1：根据目标点的坐标(Xgeo，Ygeo)和DEM坐标变换参数按照公式得到目标点在DEM格网中的行列数(RDEM,CDEM)，式中(X0DEM,Y0DEM)是DEM数据左上角的地理坐标，CellsizeDEM是DEM数据相邻两点之间的间隔。

[0057] 步骤1.2.2：根据行列数(RDEM,CDEM)得到DEM格网中最邻近的四个高程点Zzsh,Zysh,Zyx,Zzx，然后根据公式 内插出该点的高程值Z，其中dx,dy为该点距离高程点Zzsh在x、y方向上距离；

[0058] 步骤1.3：由地面点坐标(Xgeo，Ygeo，Z)根据RPC模型反算出该地面点在原始影像上的位置：

[0059] 步骤1.3.1：将地面点坐标按照RPC参数中的空间坐标标准化参数，按照公式计算得到标准化之后的空间坐标(U,V,W)，其中 λ0,h0为标准化平移参数， λs,hs为标准化比例参数，它们与RPC有理多项式模型的80个系数共同保存在卫星厂家提供给用户的RPC文件中；

[0060] 步骤1.3.2：根据RPC有理多项式模型，按照公式 计算该点在原始影像上的标准化坐标(s,l)；其中，Nums(U,V,W),Dens(U,V,W),Numl(U,V,W),Denl(U,V,W)均为如p＝a1+a2V+a3U+a4W+a5VU+a6VW+a7UW+a8V2+a9U2+a10W2+a11UVW+a12V3+a13VU2+a14VW2+a15V2U+a16U3+a17UW2+a18V2W+a19U2W+a20W3形式的多项式，RPC参数中一共有四组多项式系数分别对应以上四项数据；

[0061] 步骤1.3.3：根据像点坐标标准化参数，按照公式 对标准化坐标(s,l)进行处理，得到原始影像上的像点坐标(S,L)。其中S0,L0为标准化平移参数，Ss,Ls为标准化比例参数，它们与RPC有理多项式模型的80个系数共同保存在卫星厂家提供给用户的RPC文件中。

[0062] 步骤1.4：基于数字高程模型DEM迭代求解，计算正射校正后影像的大小。

[0063] RPC有理多项式模型虽然建立了从三维地理坐标到原始影像平面坐标的映射关系，但并没有给出从原始影像平面坐标到地理坐标的数学模型，而这一过程是从二维坐标到三维坐标的转换过程，因此只能基于数字高程模型进行迭代求解，具体过程如下：

[0064] 步骤1.4.1：根据步骤1.1计算得到原始影像左上角像点(0,0)的初始地理坐标(X0，Y0)，根据步骤1.2利用双线性内插方法在DEM上内插得到点(X0，Y0)的高程Z0，根据步骤1.3得到地面点(X0,Y0,Z0)在原始影像上的像点坐标(S0,L0)；

[0065] 步骤1.4.2：计算像点坐标(S0,L0)对应的初始地理坐标，对应的初始地理坐标(Xi，Yi)，内插DEM得到高程Zi,得到地面点(Xi,Yi,Zi)在原始影像上的像点坐标(Si,Li)，迭代求解，直至前后两次计算过程中的|Zi-Zi-1|＜δ，结束迭代。则原始影像左上角像点(0,0)对应的地面点坐标为(Xi,Yi,Zi)；

[0066] 步骤1.4.3：按照步骤1.4.1与步骤1.4.2迭代求解，直至前后两次计算中内插DEM的高程值的差值小于阈值0.1米,则最后一次计算得到的地面点坐标为左上角像素点对应的地面点坐标(Xmin,Ymax)；

[0067] 步骤1.4.4：按照步骤1.4.1至步骤1.4.3分别计算原始影像四个角点对应地面点的坐标为(Xzsh,Yzsh),(Xysh,Yysh),(Xyx,Yyx),(Xzx,Yzx)；

[0068] 之后，按照公式 计算得到X、Y方向的极值；

[0069] 然后按照公式 计算得到正射矫正后正射影像的行数widthdst和列数heightdst，其中GSD是正射校正后影像的分辨率。

[0070] 步骤2：利用反解法逐像素反算每个像素对应原始影像上像素的坐标；

[0071] 根据步骤1解算得到的正射校正后正射影像的大小和范围，根据RPC参数模型利用反解法逐像素解算每个像素对应原始影像上像素的位置，过程如图3所示：

[0072] 步骤2.1：根据当前像素在正射校正后影像中的行列坐标(row,col)，按照公式计算得到当前像素的地面点坐标(X，Y)；其中GSD为正射影像的分辨率，(row,col)为当前线程索引对应像素所在的行列坐标，(Xzsh,Yzsh)为正射校正后影像左上角点坐标，heightdst为正射校正后影像的行数。

[0073] 步骤2.2：根据地面点坐标和数字高程模型DEM，采用双线性内插法，按照步骤1.2从DEM上内插出点(X，Y)的高程值Z；

[0074] 步骤2.3：根据RPC模型，按照步骤1.3反算出地面点坐标(X,Y,Z)在原始影像上的像点坐标(S,L)；并根据公式 判断坐标(S,L)是否在原始影像内部，其中widthsrc和heightsrc分别为原始卫星影像的列数和行数。如果在影像内部则保留当前像素(row,col)对应的原始影像像素坐标(S,L)，如果在外部则将S、L的值其设置为-1；

[0075] 步骤2.4：将所有像素按照步骤2.1～2.3计算得到其对应的原始影像上的像点坐标(Si,Lj)。

[0076] 步骤3：逐像素统计并判断是否为拉花像素点，形成拉花区域二值标记图像并进行矢量化，得到拉花区域矢量文件。

[0077] 卫星影像正射校正时根据地面点位置计算原始影像上对应的像素，然后进行灰度重采样，在被遮挡的区域或是成像信息匮乏的区域重采样将会过于稠密或是重复采样，这样就会造成正射校正后的影像出现拉伸现象，如果拉伸过度就会造成纹理失真的现象，我们称之为“拉花”现象，重复采样的像素即为拉花像素。我们通过统计像素与周边像素的重叠关系来进行拉花像素的判定，详细步骤如图4所示：

[0078] 步骤3.1：以像素点(row,col)为中心，建立一个适当大小方形窗口Win[Size,Size]，此处Size的值取为30，当前像素对应原始影像上的像点坐标为(S,L)；

[0079] 步骤3.2：设近似像素个数num＝0，遍历窗口中每一个像素(ri,ci)∈Win[Size,Size]对应的原始影像上像点坐标(Si,Li)与当前像素对应的像点坐标(S,L)对比，当满足公式 时，令num＝num+1，此时阈值δ＝0.8；

[0080] 步骤3.3：当num大于阈值λ＝3时，在检测结果影像上将其标记为拉花像素，在检测结果影像上将其标记为拉花像素(灰度值为255)，否则标记为非拉花像素(灰度值为0)；

[0081] 步骤3.4：按照步骤3.1～3.3，逐像素统计其对应原始影像像点坐标(S,L)与周边像素对应原始影像像点坐标(Si,Li)的关系，找出纠正后正射影像上的所有拉花像素与其构成的拉花区域二值标记图像；

[0082] 步骤3.5：对拉花区域二值标记图像进行矢量化，得到拉花区域矢量文件，如图5与图6所示。

[0083] 本发明所述方案利用原始光学卫星影像及其RPC模型参数和测区数字高程模型DEM，通过逐像素反解其对应原始卫星影像上的坐标，并进行统计比对，判定其是否为拉花像素，最后将生成的拉花区域二值图像进行矢量化得到拉花区域矢量数据，达到光学遥感卫星正射影像上拉花区域自动检测的目的，实现了光学遥感卫星正射影像中拉花区域的自动检测。