CN108230326B - 基于gpu-cpu协同的卫星影像拉花变形快速检测方法 - Google Patents

基于gpu-cpu协同的卫星影像拉花变形快速检测方法 Download PDF

Info

Publication number
CN108230326B
CN108230326B CN201810126230.0A CN201810126230A CN108230326B CN 108230326 B CN108230326 B CN 108230326B CN 201810126230 A CN201810126230 A CN 201810126230A CN 108230326 B CN108230326 B CN 108230326B
Authority
CN
China
Prior art keywords
pixel
garland
gpu
coordinate
image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810126230.0A
Other languages
English (en)
Other versions
CN108230326A (zh
Inventor
李朋龙
丁忆
胡艳
张泽烈
徐永书
李静
罗鼎
张灿
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing geographic information and Remote Sensing Application Center (Chongqing surveying and mapping product quality inspection and testing center)
Original Assignee
CHONGQING GEOGRAPHICAL INFORMATION CENTER
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CHONGQING GEOGRAPHICAL INFORMATION CENTER filed Critical CHONGQING GEOGRAPHICAL INFORMATION CENTER
Priority to CN201810126230.0A priority Critical patent/CN108230326B/zh
Publication of CN108230326A publication Critical patent/CN108230326A/zh
Application granted granted Critical
Publication of CN108230326B publication Critical patent/CN108230326B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/60—Memory management
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/28—Indexing scheme for image data processing or generation, in general involving image processing hardware
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Abstract

本发明公开了一种基于GPU‑CPU协同的卫星影像拉花变形快速检测方法，包括CPU端加载原始卫片、RPC参数和数字高程模型DEM，迭代求解校正后正射影像的大小和范围；CPU端设计GPU线程格网的分配，并将相应数据从内存拷贝入对应显存；GPU端按线程格网并行计算每个像素对应原始卫片上的像素坐标；GPU端统计当前像素与周边像素重叠次数，并进行拉花像素判断；GPU端按线程格网对二值图像进行腐蚀处理与膨胀处理；CPU端将二值图像结果从全局存储器拷到内存并矢量化，获得拉花变形区域矢量数据等步骤。其显著效果是：实现了卫星影像中拉花变形区域的快速自动检测，大大提高了光学遥感卫星正射影像中拉花变形的查找和质检效率。

Description

基于GPU-CPU协同的卫星影像拉花变形快速检测方法
技术领域
本发明涉及到光学遥感卫星影像处理技术领域，具体地说，是一种基于GPU-CPU协同的卫星影像拉花变形快速检测方法。
背景技术
数字正射影像图(DOM)是利用光学遥感影像结合地形数据(数字高程模型DEM)根据相应的数学模型经过逐像素计算地面点与原始光学遥感影像像素的关系，并进行灰度重采样得到的既有正确的位置信息又有丰富的纹理信息的影像图。制作数字正射影像图的原始光学遥感影像根据获取方式可以分为两类：光学卫星遥感影像和低空航摄影像。两种光学影像的成像方式、数学模型和处理方法均不相同，同时由于卫星遥感影像具有获取速度快、影像框幅大、覆盖范围广的特点被广泛应用于国土资源监测、地理国情监测更新等国家重大战略项目中，在城市建设用地监测、城市违法建筑监测、城市规划管理、生态环境监测保护等领域也发挥着越来越重要的作用。
光学卫星遥感影像的正射校正是制作数字正射影像图的关键一环，正射纠正是利用光学卫星影像及其对应的数学模型，如有理多项式模型(rational polynomialcoefficients,RPC)和数字高程模型DEM数据共同消除原始卫星影像中的各种畸变(如投影差)而得到一幅既有地理坐标位置信息又有纹理信息的新影像的过程。地形起伏和卫星影像成像时光学镜头中心投影的成像方式不能保证所有地面点都能在卫星影像上成像，如坡势较陡的山坡可能会被山顶所遮挡。而在卫星影像正射校正时根据地面点位置计算原始影像上对应的像素，然后进行灰度重采样，在被遮挡的区域或是成像信息匮乏的区域重采样将会过于稠密或是重复采样，这样就会造成正射校正后影像出现拉伸现象，如果拉伸过度就会造成纹理失真的现象，我们称之为“拉花变形”，纹理失真的区域我们称之为“拉花变形区域”。拉花变形造成的纹理失真直接影响着数字正射影像图的质量，特别是山地区域拉花变形现象特别严重，直接影响着数字正射影像图制作的效率。
目前，在国内外文献和专利中，还没有针对光学卫星正射影像拉花变形的自动检测方法，在正常的数字正射影像制作中，如果出现拉花变形，需要人工目视辨别和寻找，然后再通过修改DEM再纠正的方式进行处理。而光学卫星遥感影像框幅较大，如Worldview2卫星一景影像的覆盖范围可达280平方公里，依靠人工目视辨别查找拉花区域的方法将非常地耗时，同时也可能因为人工原因造成遗漏，因此光学遥感卫星正射影像拉花区域的查找效率非常低下。
发明内容
针对现有技术的不足，本发明的目的是提供一种基于GPU-CPU协同的卫星影像拉花变形快速检测方法，基于光学卫星影像RPC模型和测区数字高程模型DEM，并利用GPU-CPU协同处理技术能够快速准确地检测出纠正后正射影像上的拉花变形区域，达到拉花变形自动检测的目的。
为达到上述目的，本发明采用的技术方案如下：
一种基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其关键在于按照以下步骤进行：
步骤1：CPU端加载原始卫片、RPC参数和数字高程模型DEM，并迭代求解正射校正后正射影像的大小和范围，计算步骤为：
步骤1-1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；
步骤1-2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；
步骤1-3：由地面点坐标根据RPC模型反算出该地面点在原始影像上的位置；
步骤1-4：基于数字高程模型DEM进行迭代求解，计算正射校正后影像的大小；
步骤2：CPU端根据校正后正射影像大小来设计GPU线程格网，开辟相应的显存，并将相应数据从内存拷贝入对应显存；
步骤3：GPU端按线程格网从校正后影像像素反算其对应原始卫片上的像素坐标；
步骤4：GPU端按照线程格网统计当前像素与周边像素重叠次数，并判断其是否为拉花像素；
步骤5：GPU端按线程格网对二值图像进行腐蚀处理与膨胀处理；
步骤6：CPU端将检测后的二值图像结果从全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据。
进一步的，步骤1-2中所述地面点高程值的计算步骤为：
步骤1-2-1：根据目标点的坐标和数字高程模型DEM，计算得到目标点在数字高程模型DEM格网中的行列数；
步骤1-2-2：根据得到的行列数计算数字高程模型DEM格网中最邻近的四个高程点，并内插出目标点的高程值。
再进一步的，所述步骤1-3中反算地面点在原始影像上位置的步骤为：
步骤1-3-1：将地面点坐标按照RPC参数中的空间坐标标准化参数计算得到标准化之后的空间坐标；
步骤1-3-2：根据RPC有理多项式模型计算该点在原始影像上的标准化坐标；
步骤1-3-3：根据像点坐标标准化参数对标准化坐标进行处理，得到原始影像上的像点坐标。
更进一步的，所述步骤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方向的极值，计算得到正射校正后正射影像的行数和列数。
进一步的，步骤3中所述原始卫片上的像素坐标求解步骤为：
步骤3-1：GPU端多线程同时执行，根据当前线程索引计算当前像素的地面点坐标；
步骤3-2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点的高程值；
步骤3-3：计算地面点在原始影像上的像点坐标，并判断该坐标是否在原始影像内部；
步骤3-4：等待所有线程计算完毕后，将整幅正射影像反算得到的每个像素对应原始影像像点位置存储在全局存储器中。
进一步的，步骤4中拉花像素的判断步骤为：
步骤4-1：按照线程索引以当前像素为中心，建立一个适当大小的方形窗口；
步骤4-2：遍历和统计方形窗口中所有像素对应原始影像上像点坐标和当前中心像素对应原始影像上像点坐标的重叠次数；
步骤4-3：当重叠次数大于阈值时，在检测结果影像上将当前中心像素标记为拉花像素，否则标记为非拉花像素；
步骤4-4：等待所有线程计算完毕，并将拉花像素判定结果存储到全局存储器中。
进一步的，所述步骤5中所述腐蚀处理的具体步骤为：
步骤A：按照线程索引以当前像素为中心，建立一个适当大小的方形腐蚀窗口，若当前中心像素灰度值为g1时进入步骤B；
步骤B：遍历腐蚀窗口中每一个像素，统计灰度值为g1的像素个数，若其小于腐蚀窗口内像素总个数的一半时，将当前中心像素赋值为0；
步骤C：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
进一步的，步骤5中所述膨胀处理的具体步骤为：
步骤S1：按照线程索引以当前像素为中心，建立一个适当大小的方形膨胀窗口；
步骤S2：将当前像素对应的灰度值为与灰度阈值g2进行比较，若相等则遍历膨胀窗口中每一个像素，并将其对应的灰度值赋值为灰度阈值g2；
步骤S3：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
进一步的，步骤6中所述矢量数据的获得步骤为：
步骤6-1：CPU端将拉花检测后的二值图像结果从全局存储器拷贝到内存中，同时释放已经开辟的GPU存储器；
步骤6-2：对二值图像中拉花变形区域进行边缘检测及提取，并存储为矢量数据。
本发明首先CPU端加载原始卫星影像及其RPC参数和测区DEM，并迭代计算正射校正后正射影像的大小和范围；然后CPU端根据计算任务分配GPU线程格网，并开辟多种显存将DEM数据、RPC参数等数据从内存拷入相应的显存；接着GPU端根据RPC有理多项式模型利用反解法解算正射影像上每个像素对应原始卫片的像素的像点坐标；其次，GPU端统计正射影像上每个像素与周围像素所对应原始卫片上像素的重叠次数，并判断其是否为拉花像素，生成拉花区域二值图像；然后GPU端对拉花变形区域二值图像进行包括腐蚀处理和膨胀处理的图像形态学处理；最后CPU端将GPU端处理得到的二值图像从显存拷出至内存，释放显存资源，并在内存中对拉花变形区域二值图像进行矢量化处理，得到拉花变形区域矢量数据。
本发明的显著效果是：利用CPU-GPU并行处理技术实现了卫星影像中拉花变形区域的快速自动检测，效率是传统CPU串行计算方法的55倍，解决了传统人工目视辨别拉花变形区域费时费力，效率低下的问题，大大提高了光学遥感卫星正射影像中拉花变形的查找和质检效率。
附图说明
图1是本发明的方法流程图；
图2是卫星正射影像拉花变形检测结果；
图3是局部拉花变形检测结果1；
图4是局部拉花变形检测结果2。
具体实施方式
下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。
本实施例结合某多山地区一景分辨率为0.5米的WorldView2卫星影像，图幅大小为36821*36132，大小为3.3GB，覆盖地面范围约330Km2，对本发明方法进行详细说明。
如图1所示，一种基于GPU-CPU协同的卫星影像拉花变形快速检测方法，具体步骤如下：
步骤1：CPU端加载原始卫片、RPC参数和数字高程模型DEM，并迭代求解正射校正后正射影像的范围大小；
具体实施时，所述正射校正后正射影像的大小和范围的计算步骤为：
步骤1-1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；
光学卫星遥感影像不仅有对应的RPC模型参数，同时也自带了原始影像上地理坐标与影像行列数之间变换的6个仿射变换参数，存储在数组adfGeoTransform[6]中，按照公式(1)可以通过仿射变换参数从一个像元的行列数得到其对应的地理坐标。
式中，(row,col)是像素的行列坐标，(Xgeo，Ygeo)为该像素对应的地理坐标其中Xgeo为经度，Ygeo为维度,(adfGeoTransform[0],adfGeoTransform[3])为该影像左上角像素对应的地理坐标，是行列坐标到地理坐标的四个仿射变换矩阵。
步骤1-2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；
数字高程模型DEM是将地形数据按照一定的间隔的格网进行规则存储，DEM中每个点的行列数与该点的地理坐标相关，数值则是该点的高程Z，可以根据一个点的地理坐标(Xgeo，Ygeo)和DEM来内插该点对应的高程值。DEM内插方法分为三类：最邻近内插法，双线性内插法和三次卷积内插法，三种方法的计算量依次变大，本例中采用双线性内插法，计算步骤如下：
步骤1-2-1：根据目标点的坐标(Xgeo，Ygeo)和DEM坐标变换参数按照公式(2)得到目标点在DEM格网中的行列数(RDEM,CDEM)，式中(X0DEM,Y0DEM)是DEM数据左上角的地理坐标，CellsizeDEM是DEM数据相邻两点之间的间隔。
步骤1-2-2：根据(RDEM,CDEM)得到DEM格网中最邻近的四个高程点Zzsh,Zysh,Zyx,Zzx，然后根据公式(3)内插出该点的高程值Z，其中dx,dy为该点距离高程点Zzsh在x,y方向上距离。
步骤1-3：由地面点坐标(Xgeo,Ygeo,Z)根据RPC模型反算出该地面点在原始影像上的位置，计算步骤为：
步骤1-3-1：将地面点坐标(Xgeo,Ygeo,Z)按照RPC参数中的空间坐标标准化参数根据公式(4)进行处理得到标准化之后的空间坐标(U,V,W)。
步骤1-3-2：根据RPC有理多项式模型，按照公式(5)计算该点在原始影像上的标准化坐标(s,l)；
其中，Nums(U,V,W),Dens(U,V,W),Numl(U,V,W),Denl(U,V,W)均为如公式(6)所示的多项式，RPC参数中一共有四组多项式系数分别对应以上四项数据：
步骤1-3-3：根据像点坐标标准化参数，按照公式(7)对(s,l)进行处理，得到原始影像上像点坐标(S,L)，
以上公式中，λ0,h0,S0,L0为标准化平移参数，λs,hs,Ss,Ls为标准化比例参数，它们与RPC有理多项式模型的80个系数共同保存在卫星厂家提供给用户的RPC文件中。
步骤1-4：基于数字高程模型DEM进行迭代求解，计算正射校正后影像的大小。
RPC有理多项式模型虽然建立了从三维地理坐标到原始影像平面坐标的映射关系，但并没有给出从原始影像平面坐标到地理坐标的数学模型，而这一过程是从二维坐标到三维坐标的转换过程，因此只能基于数字高程模型进行迭代求解，具体过程如下：
步骤1-4-1：根据步骤1-1计算得到原始影像左上角像点(0,0)的初始地理坐标(X0，Y0)，根据步骤1-2利用利用双线性内插方法在DEM上内插得到点(X0，Y0)的高程值Z0，根据步骤1-3得到地面点(X0,Y0,Z0)在原始影像上的像点坐标(S0,L0)；
步骤1-4-2：再根据以上过程计算像点坐标(S0,L0)，对应的初始地理坐标(Xi，Yi)，内插DEM得到高程Zi,得到地面点(Xi,Yi,Zi)在原始影像上的像点坐标(Si,Li)，迭代求解，直至前后两次计算过程中的|Zi-Zi-1|＜α＝0.1米，结束迭代。则原始影像左上角像点(0,0)对应的地面点坐标为(Xi,Yi,Zi)；
步骤1-4-3：按照步骤1-4-1与步骤1-4-2迭代求解，直至前后两次计算中内插DEM的高程值的差值小于阈值，则最后一次计算得到的地面点坐标为左上角像素点对应的地面点坐标(Xmin,Ymax)；
步骤1-4-4：按照步骤1-4-1到步骤1-4-3分别计算原始影像四个角点对应地面点的坐标为(Xzsh,Yzsh),(Xysh,Yysh),(Xyx,Yyx),(Xzx,Yzx)，然后根据公式(8)得到X,Y方向的极值，正射校正后正射影像左上角对应的地理坐标为(Xzsh,Yzsh)；
步骤1-4-5：根据所得的X、Y方向的极值，按照公式(9)计算得到正射校正后正射影像的行数和列数，其中GSD是正射校正后影像的分辨率。
步骤2：CPU端根据校正后正射影像大小来设计GPU线程格网的分配，开辟相应的显存，并将相应数据从内存拷贝入对应显存；
GPU并行计算是依托GPU多线程之间的并行计算，而GPU线程则是通过线程格网进行组织和调度的。光学卫星影像的正射校正是对影像每个像素进行重采样的过程，因此总的计算任务为widthdst*heightdst，也就是说需要widthdst*heightdst个GPU线程，每个线程对应一个像素。线程格网可以按照一维、二维或是三维组织，假设按照二维线程格网进行组织，线程块的大小为dimBlock(N,N)，即每个线程块中有N*N个线程，格网dimGrid的大小可以根据公式(10)计算得到：
然后，GPU-CPU协同处理技术中，CPU和GPU均有自己的存储地址空间，且相互之间可通信。GPU共有六种存储器：寄存器、局部存储器、共享存储器、全局存储器、常数存储器和纹理存储器，每种存储器的大小和访问速度均不相同，构成了GPU多层次性存储的特点。线程格网中的每条线程拥有仅供自己访问的寄存器与局部存储器，同一个线程块中的线程能够读写块内的共享存储器，全局存储器可以供格网内的所有线程读写，而常数存储器和纹理存储器只能被格网内的所有线程读取，不能进行写操作。
因此对于不同大小的数据及其访问方式，选择不同的存储器可以有效提高数据的访问速度和处理能力，因此本方法中将供每条线程访问且数据量大的DEM存储在纹理存储器中，利用纹理存储器缓存加速的特点提高访问效率，将每条线程读取且数据量较小的RPC参数存储在访问速度较快的常数存储器中，将大量的中间计算结果存放在全局存储器中。
进入步骤3：GPU端按线程格网从校正后影像像素并行反算其对应原始卫片上的像素坐标，计算步骤如下：
步骤3-1：GPU端多线程同时执行，每条线程对应了一个像素，利用公式(11)根据当前线程索引计算当前像素的地面点坐标(X,Y)，其中GSD为正射影像的分辨率，(row,col)为当前线程索引对应像素所在的行列坐标，(Xzsh,Yzsh)为正射校正后左上角点坐标，heightdst为正射校正后影像的行数。
步骤3-2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点(X,Y)的高程值Z，具体方法参见步骤1-2；
步骤3-3：采用步骤1-3所述方法计算地面点(X,Y,Z)在原始影像上的像点坐标(S,L)，并根据公式(12)判断坐标(S,L)是否在原始影像内部，如果在外部则将S、L的值其设置为-1，widthsrc和heightsrc分别为原始卫星影像的列数和行数。
步骤3-4：等待所有线程计算完毕后，将整幅正射影像反算得到的每个像素对应原始影像像点位置存储在全局存储器中。
步骤4：GPU端按照线程格网统计当前像素与周围像素重叠次数，并判断其是否为拉花像素并标记，判断方法如下：
步骤4-1：按照线程索引以当前像素(row,col)为中心，建立一个大小为20*20方形窗口Win1[20,20]，当前像素对应原始影像上的像点坐标为(S,L)；
步骤4-2：设重叠像素个数num＝0，遍历窗口中每一个像素(ri,ci)∈Win1[20,20]对应的原始影像上像点坐标(Si,Li)与当前像素对应的(S,L)对比，此时δ＝0.8，若满足公式(13)则num＝num+1。
步骤4-3：当重叠次数num大于阈值λ＝2时，将该像素视为拉花像素，在检测结果影像上将其标记为拉花像素(灰度值为255)，否则标记为非拉花像素(灰度值为0)；
步骤4-4：等待所有线程计算完毕，并将拉花像素判定结果存储到全局存储器中。
步骤5：GPU端按线程格网对二值图像进行腐蚀处理与膨胀处理；
其中所述腐蚀处理的具体步骤为：
步骤A：按照线程索引以当前像素(row,col)为中心，建立一个大小为3*3的方形腐蚀窗口Win2[3,3]，若当前中心像素灰度值为g1＝255时进入步骤B；
步骤B：遍历腐蚀窗口中每一个像素(ri,ci)∈Win2[3,3]，统计灰度值为g1＝255的像素个数sum，若sum小于4，则该中心像素赋值为0；
步骤C：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
所述膨胀处理的具体步骤为：
步骤S1：按照线程索引以当前像素(row,col)为中心，建立一个大小为5*5的方形膨胀窗口Win3[5,5]，当前像素对应的灰度值为gray0；
步骤S2：将当前像素对应的灰度值为与灰度阈值g2＝255进行比较，若gray0＝g2＝255，则遍历膨胀窗口中每一个像素(ri,ci)∈Win3[5,5]，并将其对应的灰度值赋值为灰度阈值g2＝255；
步骤S3：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
步骤6：CPU端将检测后的二值图像结果从全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据：
步骤6-1：CPU端将拉花检测后的二值图像结果从全局存储器拷贝到内存中，同时释放已经开辟的GPU存储器；
步骤6-2：对二值图像中拉花变形区域进行边缘检测及提取，并存储为矢量数据。
本发明首先CPU端加载原始卫星影像及其RPC参数和测区DEM，并迭代计算正射校正后正射影像的大小和范围；然后CPU端根据计算任务分配GPU线程格网，并开辟多种显存将DEM数据、RPC参数等数据从内存拷入相应的显存；接着GPU端根据RPC有理多项式模型利用反解法解算正射影像上每个像素对应原始卫片的像素的像点坐标；其次，GPU端统计正射影像上每个像素与周围像素所对应原始卫片上像素的重叠次数，并判断其是否为拉花像素，生成拉花区域二值图像；然后GPU端对拉花变形区域二值图像进行包括腐蚀处理和膨胀处理的图像形态学处理；最后CPU端将GPU端处理得到的二值图像从显存拷出至内存，释放显存资源，并在内存中对拉花变形区域二值图像进行矢量化处理，得到拉花变形区域矢量数据。

Claims (9)

1.一种基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于按照以下步骤进行处理：
步骤1：CPU端加载原始影像、RPC参数和数字高程模型DEM，并迭代求解正射校正后正射影像的大小和范围，计算步骤为：
步骤1-1：通过仿射变换参数对原始影像行列数仿射变换，获得其对应的地理坐标；
步骤1-2：根据像点地理坐标和数字高程模型DEM，采用双线性内插法内插出像点对应地面点的高程值；
步骤1-3：由地面点坐标根据RPC模型反算出该地面点在原始影像上的位置；
步骤1-4：基于数字高程模型DEM进行迭代求解，计算正射校正后影像的大小；
步骤2：CPU端根据校正后正射影像大小来设计GPU线程格网，开辟相应的显存，并将相应数据从内存拷贝入对应显存；
步骤3：GPU端按线程格网从校正后影像像素反算其对应原始影像上的像素坐标；
步骤4：GPU端按照线程格网统计当前像素与周边像素重叠次数，并判断其是否为拉花像素；
步骤5：GPU端按线程格网对二值图像进行腐蚀处理与膨胀处理；
步骤6：CPU端将检测后的二值图像结果从全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据。
2.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：步骤1-2中所述地面点高程值的计算步骤为：
步骤1-2-1：根据目标点的坐标和数字高程模型DEM，计算得到目标点在数字高程模型DEM格网中的行列数；
步骤1-2-2：根据得到的行列数计算数字高程模型DEM格网中最邻近的四个高程点，并内插出目标点的高程值。
3.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：所述步骤1-3中反算地面点在原始影像上位置的步骤为：
步骤1-3-1：将地面点坐标按照RPC参数中的空间坐标标准化参数计算得到标准化之后的空间坐标；
步骤1-3-2：根据RPC模型计算该点在原始影像上的标准化坐标；
步骤1-3-3：根据像点坐标标准化参数对标准化坐标进行处理，得到原始影像上的像点坐标。
4.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：所述步骤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方向的极值，计算得到正射校正后正射影像的行数和列数。
5.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：步骤3中所述原始影像上的像素坐标求解步骤为：
步骤3-1：GPU端多线程同时执行，根据当前线程索引计算当前像素的地面点坐标；
步骤3-2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点的高程值；
步骤3-3：计算地面点在原始影像上的像点坐标，并判断该坐标是否在原始影像内部；
步骤3-4：等待所有线程计算完毕后，将整幅正射影像反算得到的每个像素对应原始影像像点位置存储在全局存储器中。
6.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：步骤4中拉花像素的判断步骤为：
步骤4-1：按照线程索引以当前像素为中心，建立一个适当大小的方形窗口；
步骤4-2：遍历和统计方形窗口中每一个像素对应的原始影像上像点坐标和周围像素重叠次数；
步骤4-3：当重叠次数大于阈值时，在检测结果影像上将当前中心像素标记为拉花像素，否则标记为非拉花像素；
步骤4-4：等待所有线程计算完毕，并将拉花像素判定结果存储到全局存储器中。
7.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：所述步骤5中所述腐蚀处理的具体步骤为：
步骤A：按照线程索引以当前像素为中心，建立一个适当大小的方形腐蚀窗口，若当前中心像素灰度值为g1时进入步骤B；
步骤B：遍历腐蚀窗口中每一个像素，统计灰度值为g1的像素个数，若其小于腐蚀窗口内像素总个数的一半时，将当前中心像素赋值为0；
步骤C：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
8.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：步骤5中所述膨胀处理的具体步骤为：
步骤S1：按照线程索引以当前像素为中心，建立一个适当大小的方形膨胀窗口；
步骤S2：将当前像素对应的灰度值为与灰度阈值g2进行比较，若相等则遍历膨胀窗口中每一个像素，并将其对应的灰度值赋值为灰度阈值g2；
步骤S3：等待所有线程计算完毕，并将计算结果存储到全局存储器中。
9.根据权利要求1所述的基于GPU-CPU协同的卫星影像拉花变形快速检测方法，其特征在于：步骤6中所述矢量数据的获得步骤为：
步骤6-1：CPU端将拉花检测后的二值图像结果从全局存储器拷贝到内存中，同时释放已经开辟的显存；
步骤6-2：对二值图像中拉花变形区域进行边缘检测及提取，并存储为矢量数据。
CN201810126230.0A 2018-02-08 2018-02-08 基于gpu-cpu协同的卫星影像拉花变形快速检测方法 Active CN108230326B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810126230.0A CN108230326B (zh) 2018-02-08 2018-02-08 基于gpu-cpu协同的卫星影像拉花变形快速检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810126230.0A CN108230326B (zh) 2018-02-08 2018-02-08 基于gpu-cpu协同的卫星影像拉花变形快速检测方法

Publications (2)

Publication Number Publication Date
CN108230326A CN108230326A (zh) 2018-06-29
CN108230326B true CN108230326B (zh) 2018-11-30

Family

ID=62670992

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810126230.0A Active CN108230326B (zh) 2018-02-08 2018-02-08 基于gpu-cpu协同的卫星影像拉花变形快速检测方法

Country Status (1)

Country Link
CN (1) CN108230326B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109448119B (zh) * 2018-10-19 2022-04-19 深圳市工勘岩土集团有限公司 一种地理信息系统开发中关于数字高程模型（dem）应用的方法
CN110796734B (zh) * 2019-10-31 2024-01-26 中国民航科学技术研究院 基于高分辨率卫星技术的机场净空巡查方法及装置
CN110889949B (zh) * 2019-12-09 2021-08-31 国网湖南省电力有限公司 基于分层gpu计算的输电线路卫星监测数据处理方法及系统
CN112050735B (zh) * 2020-09-09 2022-02-22 中国科学院空天信息创新研究院 光学遥感卫星大数据的递推精化地面位置方法及存储介质
CN113205090B (zh) * 2021-04-29 2023-10-24 北京百度网讯科技有限公司 图片矫正方法、装置、电子设备及计算机可读存储介质
CN113836751B (zh) * 2021-11-22 2022-02-08 武汉峰岭科技有限公司 一种基于等高线拉伸的数字高程模型自动调整方法及系统
CN115423696A (zh) * 2022-07-29 2022-12-02 上海海洋大学 一种自适应线程参数的遥感正射影像并行生成方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104180794A (zh) * 2014-09-02 2014-12-03 西安煤航信息产业有限公司 数字正射影像拉花区域的处理方法
CN106815807A (zh) * 2017-01-11 2017-06-09 重庆市地理信息中心 一种基于gpu‑cpu协同的无人机影像快速镶嵌方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150160371A1 (en) * 2013-12-06 2015-06-11 Schlumberger Technology Corporation Gpu accelerated deflation in geomechanics simulator

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104180794A (zh) * 2014-09-02 2014-12-03 西安煤航信息产业有限公司 数字正射影像拉花区域的处理方法
CN106815807A (zh) * 2017-01-11 2017-06-09 重庆市地理信息中心 一种基于gpu‑cpu协同的无人机影像快速镶嵌方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
CPU和GPU协同处理的光学卫星遥感影像正射校正方法;方留杨等;《测绘学报》;20131031;第42卷(第5期);第668-675页 *
Range Cell Migration Correction Using Texture Mapping on GPU;Bin Liu等;《IEEE 10th International Conference on Signal Processing Proceedings》;20101203;第2172-2175页 *

Also Published As

Publication number Publication date
CN108230326A (zh) 2018-06-29

Similar Documents

Publication Publication Date Title
CN108230326B (zh) 基于gpu-cpu协同的卫星影像拉花变形快速检测方法
CN104574347B (zh) 基于多源遥感数据的在轨卫星图像几何定位精度评价方法
CN106815807B (zh) 一种基于gpu‑cpu协同的无人机影像快速镶嵌方法
CN110569797B (zh) 地球静止轨道卫星影像山火检测方法、系统及其存储介质
CN104299228A (zh) 一种基于精确点位预测模型的遥感影像密集匹配方法
CN112100301A (zh) 一种利用高分遥感技术实现水域岸线动态监测的方法
CN108335261B (zh) 一种光学遥感卫星正射影像拉花区域自动检测方法
CN114022783A (zh) 基于卫星图像的水土保持生态功能遥感监测方法和装置
CN104180794B (zh) 数字正射影像拉花区域的处理方法
CN108919319A (zh) 海岛礁卫星影像无地面控制点定位方法及系统
CN111003214B (zh) 基于云控制的国产陆地观测卫星姿轨精化方法
CN108269228B (zh) 基于gpu并行计算的无人机影像拉花区域自动探测方法
CN104574338A (zh) 基于多角度线阵ccd传感器的遥感图像超分辨率重建方法
CN104200527B (zh) 一种真正射影像的生成方法
CN109579796B (zh) 一种投影后影像的区域网平差方法
CN107945218B (zh) 边缘大畸变影像匹配方法
CN108257130B (zh) 一种航空正射影像全景图拉花区域快速检测方法
Liu et al. Automatic building height estimation with shadow correction over heterogeneous compact cities using stereo Gaofen-7 data at sub-meter resolution
Ye et al. Aerial Image Stitching Method Based on Feature Transfer and Tile Image
CN105096275A (zh) 基于暗通道原理的单幅山脉遥感图像高程值提取方法
CN113446992B (zh) 地形勘测中的地形勘测点布测优化方法
CN106856006A (zh) 二维等间距投影图的数字化方法
Li et al. An Automatic Detection Method for Tensile Deformations in Optical Satellite Orthophoto
Peng et al. Block Orthorectification and Mosaic of 1960s Disp Images in Qujing, Yunnan
Bai et al. Weak Target Detection in High-Resolution Remote Sensing Images by Combining Super-Resolution and Deformable FPN

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CP01 Change in the name or title of a patent holder

Address after: 400020 Jiangbei District, Chongqing electric measuring Village No. 231

Patentee after: Chongqing geographic information and Remote Sensing Application Center (Chongqing surveying and mapping product quality inspection and testing center)

Address before: 400020 Jiangbei District, Chongqing electric measuring Village No. 231

Patentee before: Chongqing Geographical Information Center

CP01 Change in the name or title of a patent holder