CN115423696B - 一种自适应线程参数的遥感正射影像并行生成方法 - Google Patents

Abstract

本发明公开了一种自适应线程参数的遥感正射影像并行生成方法，属于摄影测量技术领域。该方法包括：获取原始立体影像，根据其角点的像方坐标和DEM通过RFM模型迭代计算获得四至地理坐标，建立模拟纠正影像；根据GPU的性能参数和CUDA核函数确定线程块和线程网格参数；使用GPU并行计算，获得模拟纠正影像像元的地理坐标及地面高程值；通过RFM模型计算其在原始立体影像中的行列，并以其插值计算的灰度值作为纠正影像的灰度值，获得正射纠正影像。本发明使用自适应线程块适应不同的GPU，具有较好的性能及更加自动化。同时使用GPU并行计算，能够极大地缩短正射纠正的时间，满足对立体影像批量进行正射纠正处理的需要。

Description

一种自适应线程参数的遥感正射影像并行生成方法

技术领域

本发明涉及摄影测量技术领域，特别涉及一种自适应线程参数的遥感正射影像并行生成方法。

背景技术

随着遥感技术的发展，卫星遥感影像获取技术正趋向于三多(多平台、多传感器和多角度)和三高(高空间分辨率、高光谱分辨率和高时相分辨率)，逐渐成为地面信息采集最为有效的手段。摄影测量的研究对象也不再局限于航空影像，其理论与算法也被广泛应用于高分辨卫星，高分辨率卫星遥感影像进入国际市场以后，为航天摄影测量开辟了新的应用研究领域，成像几何模型作为遥感影像几何处理和地球空间信息提取的理论基础，已成为亟待研究解决的关键技术。

正射纠正是遥感影像数据处理的关键环节，可以消除成像过程中的各种因素导致的影像畸变，也是整个遥感数据处理过程中计算量最大、耗时最长的步骤之一。正射纠正的主要方法有共线性模型、基于放射变换的严格几何模型、改进多项式模型、有理函数模型和神经网络纠正等。上述方法可分为通用经验模型和严格物理模型，通用经验模型以控制点进行正射纠正，受地形和控制点的限制；严格物理模型需要已知传感器的轨道参数和姿态参数等等，可获得较高的精度，但是多数卫星影像不再提供传感器各项参数，而是提供RPC参数文件。

图形处理器GPU专门为计算密集型、高度并行化的处理任务设计，有效地解决海量计算的速度问题，使用GPU对正射纠正过程进行加速，能够极大地缩短对立体影像正射纠正的时间。使用基于CUDA的GPU并行计算时，需要在调用核函数时设置核函数的线程参数即线程块和线程网格的参数。对影像等数据使用二维线程参数时，常规的选择是线程块使用16*16大小，然而这种常规设置并不一定是最优的选择，以往的研究大多是测试不同线程块大小的运行时间来选择较好的一个，较为繁琐。

发明内容

本发明为了解决上述问题，提出了一种能自适应不同GPU的方法确定线程参数，基于数字高程模型和遥感影像几何纠正RPC模型对立体影像进行间接正射纠正，并使用GPU并行计算加速正射纠正过程，使用不同线程数、寄存器数和共享内存大小计算内核占用率，通过最大占用率的线程确定线程块和线程网格参数调用正射校正核函数CUDA进行GPU并行计算，对高分辨率卫星遥感立体影像高效地进行正射纠正，生成DOM影像产品，能够满足摄影测量、地理测绘等实际生产需要。

为了达到上述目的，提出了一种自适应线程参数的遥感正射影像并行生成方法，包括以下步骤：

(1)获取原始立体影像，根据其角点的像方坐标、影像的RPC文件的参数和DEM通过RFM模型进行迭代计算获得四至地理坐标，建立模拟纠正影像；

(2)根据GPU的性能参数和CUDA核函数使用的寄存器及共享内存确定线程块、线程网格参数；

(3)使用所述GPU并行计算，对所述模拟纠正影像每个像元进行坐标变换获得其地理坐标(B,L)，通过坐标变换计算对应所述DEM中的行列，插值计算获得地面高程值H；

(4)根据所述模拟纠正影像像元的地理坐标(B,L)和所述地面高程值H，通过所述RFM模型计算在所述原始立体影像中的行列，并以该行列进行插值计算的灰度值作为所述纠正影像的灰度值，获得正射纠正影像。

进一步的，所述步骤1建立模拟纠正影像包括：

(1.1)根据所述原始立体影像的左上、右上、左下和右下四个角角点像方坐标(r₀,c₀)以及所述RPC文件的参数建立RFM模型；

(1.2)以所述DEM的平均高程H₀作为初值，通过所述RFM模型计算获得初始经纬度(B₁,L₁)，坐标变换得到(B₁,L₁)在所述DEM中的行列(x₁,y₁)，插值求得高程H₁；

(1.3)将所述像方坐标(r₀,c₀)和高程H₁输入所述RFM模型，获取新的经纬度(B₂,L₂)及高程H₂，重复迭代计算，直到计算的经纬度满足误差要求，将此经纬度作为角点的地理坐标，并将四个角点的地理坐标进行对比，获得影像的最大、最小维度和最大、最小经度作为四至，计算所述模拟纠正影像的行列数，建立所述模拟纠正影像。

进一步的，所述RFM模型公式为：

其中：

Num_L(P,L,H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²

+a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³

+a₁₇PH²+a₁₈L²H+a₁₉P²H+a₂₀H³

Den_L(P,L,H)＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²

+b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³

+b₁₇PH²+b₁₈L²H+b₁₉P²H+b₂₀H³

Num_S(P,L,H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²

+c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³

+c₁₇PH²+c₁₈L²H+c₁₉P²H+c₂₀H³

Den_S(P,L,H)＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²

+d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³

+d₁₇PH²+d₁₈L²H+d₁₉P²H+d₂₀H³

(P,L,H)为标准化的地面坐标；(X,Y)为标准化的影像坐标，其标准化原理如下：

其中：LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的标准化参数；SAMP_OFF、SAMP_SCALE、LINE_OFF、LINE_SCALE为影像坐标的标准化参数，上述参数均在所述RPC文件中获得。

进一步的，所述步骤2中确定线程块和线程网格参数包括如下步骤：

(2.1)根据所述GPU的性能参数和CUDA核函数使用的寄存器计算每个线程的寄存器数限制流处理器使用的线程束数值Regs_LimitperSM为：

其中：Threads_warp为线程束包含的线程数量，Warp_alloc为线程束的分配粒度，Regs_allocsize为寄存器分配单元大小，Regs_perSM每个流处理器包含的寄存器数；Regs_perB为CUDA核函数中一个线程使用的寄存器数量；

(2.2)根据所述核函数使用的共享内存计算每个线程块使用的实际线程束Warps_perB为：

threads_perB＝64*n

其中：threads_perB为线程块大小，最少为64，最大为1024，n＝1,2,…,16；

(2.3)所述核函数使用的寄存器及共享内存计算对此线程数的GPU线程束占用率：

其中：

Blocks_warps为线程束限制流处理器可分配的线程块数：

Blocks_regs为寄存器限制流处理器可分配的线程块数：

Blocks_Smemory为共享内存限制流处理器可分配的线程块数：

其中：Smemory_perB为所述核函数中每个线程使用的共享内存，Smemory_allocsize为共享内存的分配单元大小，Smemory_LimitperSM为流处理器的共享内存，Warps_LimitperSM为线程束数量；

(2.4)遍历计算设定的线程数对应的线程束占用率，Max(Occupancy_n),n＝1,2,…,16，当占用率最大时获取最小的n值，确定所述线程块的大小threads_perB，选择二维线程索引的Block_x、Block_y由以下公式计算：

对应的二维线程网格为：

其中：Img_rows，Img_cols为所述原始立体影像的行列数；m＝1,2,…,n。

进一步的，所述步骤3中获得地理坐标(B,L)及地面高程值H包括：

将所述原始立体影像及所述四至地理坐标输入所述GPU，利用步骤2中获得的所述线程块和线程网格参数调用所述CUDA核函数对每个所述模拟纠正影像像元通过左上角地理坐标和像元像方坐标进行坐标变换获得像元的地理坐标(B,L),并根据所述DEM的左上角坐标坐标变换在DEM中的行列(dem_r,dem_c)，使用双线性插值进行重采样计算，获得所述地面高程H。

进一步的，所述步骤4中获得正射纠正影像包括：

根据模拟纠正影像像元的地理坐标(B,L)和地面高程H通过所述RFM模型计算得到所述模拟纠正影像像元对应在所述原始立体影像像方中的行列(r,c)，使用插值进行重采样计算，获得灰度值H作为纠正影像的灰度值，每个像元的地理坐标为(B,L,H)，纠正影像的行列坐标为(r,c)，使用纠正后数据生成所述正射纠正影像。

本发明的有益效果：

1、基于DEM和RPC模型进行间接纠正，对大范围、高分辨卫星遥感影像不受地形和控制点的限制进行正射纠正。

2、能自适应不同的GPU确定使用CUDA核函数时的线程块、线程网格参数，性能更好且自动化。

3、使用GPU并行计算对正射纠正过程进行加速，极大地缩短了对立体影像正射纠正的时间，同时满足对影像批量快速、实时进行处理的需要。

附图说明

图1为本发明实施例自适应线程参数的遥感正射影像并行生成方法流程示意图。

图2位本发明实施例自适应GPU确定线程块、线程网格参数的示意图。

图3为本发明实施例对高分辨率卫星影像正射纠正前后对比图。

具体实施方式

下面结合附图及较佳实施例详细说明本发明的具体实施方式。

如图1所示，本发明提供了一种自适应线程参数的遥感正射影像并行生成方法，具体包括如下步骤：

S101、获取原始立体影像，根据其角点的像方坐标、影像的RPC文件的参数和DEM通过RFM模型进行迭代计算获得四至地理坐标，建立模拟纠正影像；

根据原始影像的左上、右上、左下和右下四个角角点像方坐标(r₀,c₀)和RPC文件的参数建立RFM模型，以DEM的平均高程H₀作为初值，通过RFM模型计算获得初始经纬度(B₁,L₁)，坐标变换得到(B₁,L₁)在DEM中的行列(x₁,y₁),插值求得高程H₁，将像方坐标(r₀,c₀)和高程H₁输入RFM模型，计算得到经纬度(B₂,L₂)及高程H₂，重复迭代计算，直到计算的经纬度满足误差要求，将此经纬度作为角点的地理坐标，并将四个角点的地理坐标进行对比，获得影像的最大、最小维度和最大、最小经度作为四至，计算模拟影像的行列数，建立模拟纠正影像。

RFM模型公式如下：

Num_L(P,L,H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²+a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³+a₁₇PH²+a₁₈L²H+a₁₉P²H+b₂₀H³

Den_L(P,L,H)＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²+b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³+b₁₇PH²+b₁₈L²H+b₁₉P²H+a₂₀H³

Num_S(P,L,H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²+c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³+c₁₇PH²+c₁₈L²H+c₁₉P²H+c₂₀H³

Den_S(P,L,H)＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²+d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³+d₁₇PH²+d₁₈L²H+d₁₉P²H+d₂₀H³

(P,L,H)为标准化的地面坐标；(X,Y)为标准化的影像坐标，其标准化原理如下：

式中，LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的标准化参数；SAMP_OFF、SAMP_SCALE、LINE_OFF、LINE_SCALE为影像坐标的标准化参数，上述参数在影像的RPC文件中获得。

S102、根据GPU的性能参数和CUDA核函数使用的寄存器及共享内存确定线程块、线程网格参数；

如图2所示，从GPU的参数中可知每个流处理器使用的线程束数量为Warps_LimitperSM，一个线程束包Threads_warp个线程，线程束分配粒度为Warp_alloc，寄存器分配单元大小Regs_allocsize，每个流处理器大小的寄存器数Regs_perSM,正射纠正核函数中一个线程使用的寄存器数Regs_perB，可计算每个线程的寄存器数限制流处理器使用的线程束数值为：

核函数中每个线程使用的共享内存为Smemory_perB，共享内存的分配单元大小Smemory_allocsize，流处理器的共享内存为Smemory_LimitperSM。已有经验可知线程块的线程数最少为64，最大为1024，可设线程块大小为threads_perB，每个线程块使用的实际线程束为：

Warps_perB＝threads_perB/32

threads_perB＝64*n

其中，n＝1,2,…,16。

由线程束限制流处理器可分配的线程块数:

由寄存器限制流处理器可分配的线程块数:

由共享内存限制流处理器可分配的线程块数:

可计算对此线程数的GPU线程束占用率：

通过对设定的线程数遍历计算对应的线程束占用率，Max(Occupancy_n),n＝1,2,…,16,当占用率最大时获取最小的n值，即可确定线程块的大小threads_perB。启动纠正的核函数时CUDA通过设定线程块和网格，通过上述计算的线程块中的线程数，选择二维线程索引的Block_x、Block_y应尽可能平均，可由以下公式计算：

对应的二维线程网格为：

其中：Img_rows，Img_cols为原始影像的行列数；m＝1,2,…,n。即可确定CUDA核函数的线程块、线程网格参数。

S103、使用所述GPU并行计算，对所述模拟纠正影像每个像元进行坐标变换获得其地理坐标，通过坐标变换计算对应所述DEM中的行列，插值计算获得地面高程值；

向GPU中传输影像数据及参数，使用步骤S102中确定的线程块、线程网格参数调用CUDA核函数对每个模拟纠正影像像元通过左上角地理坐标和像元像方坐标进行坐标变换获得像元的地理坐标(B,L),并根据DEM的左上角坐标坐标变换在DEM中的行列(dem_r,dem_c)，使用双线性插值进行重采样计算，获得地面高程H。

S104、根据所述模拟纠正影像像元的地理坐标和所述地面高程值通过所述RFM模型计算在所述原始立体影像中的行列，并以该行列进行插值计算的灰度值作为所述纠正影像的灰度值，获得正射纠正影像。

根据模拟纠正影像像元的地理坐标(B,L)和地面高程H，通过RFM模型计算得到模拟纠正影像像元对应在原始立体影像像方中的行列(r,c)，使用插值进行重采样计算，获得灰度值H作为纠正影像的灰度值，每个像元的地理坐标为(B,L,H)，纠正影像的行列坐标为(r,c)，将纠正后数据传输回CPU中并生成正射纠正影像，

表1为本发明与ERDAS LPS商用软件对影像进行正射纠正的水平精度结果，对两个纠正影像提取一定数量的同名点，以ERDAS LPS纠正影像的经纬度为真值进行均方根误差计算，获得误差结果为平面精度在经度偏差1.3m、纬度偏差1.19m、综合为1.77m，效果较好，可用于对高分辨率遥感卫星影像的正射纠正使用。

表1

虽然以上描述了本发明的具体实施方式，但是本领域的技术人员应当理解，这些仅是举例说明，在不背离本发明的原理和实质的前提下，可以对这些实施方式做出多种变更或修改，因此，本发明的保护范围由所附权利要求书限定。

Claims (4)

1.一种自适应线程参数的遥感正射影像并行生成方法，其特征在于，包括以下步骤：

(1)获取原始立体影像，根据其角点的像方坐标、RPC文件的参数和DEM通过RFM模型进行迭代计算获得四至地理坐标，建立模拟纠正影像；

(2)根据GPU的性能参数和CUDA核函数使用的寄存器及共享内存确定线程块和线程网格参数；

(2.1)根据所述GPU的性能参数和CUDA核函数使用的寄存器计算每个线程的寄存器数限制流处理器使用的线程束数值Regs_LimitperSM为：

其中：Threads_warp为线程束包含的线程数量，Warp_alloc为线程束的分配粒度，Regs_allocsize为寄存器分配单元大小，Regs_perSM每个流处理器包含的寄存器数；Regs_perB为所述CUDA核函数中一个线程使用的寄存器数量；

(2.2)根据所述CUDA核函数使用的共享内存计算每个线程块使用的实际线程束Warps_perB为：

threads_perB＝64*n

其中：threads_perB为线程块大小，最少为64，最大为1024，n＝1,2,…,16；

(2.3)根据所述CUDA核函数使用的寄存器及共享内存计算对此线程数的GPU线程束占用率：

其中：

Blocks_warps为线程束限制流处理器可分配的线程块数：

Blocks_regs为寄存器限制流处理器可分配的线程块数：

Blocks_Smemory为共享内存限制流处理器可分配的线程块数：

其中：Smemory_perB为所述核函数中每个线程使用的共享内存，Smemory_allocsize为共享内存的分配单元大小，Smemory_LimitperSM为流处理器的共享内存，Warps_LimitperSM为线程束数量；

(2.4)遍历计算设定的线程数对应的线程束占用率，Max(Occupancy_n),n＝1,2,…,16，当占用率最大时获取最小的n值，确定所述线程块的大小threads_perB，选择二维线程索引的Block_x、Block_y，由以下公式计算：

对应的二维线程网格为：

其中：Img_rows，Img_cols为所述原始立体影像的行列数；m＝1,2,…,n；

(3)使用所述GPU并行计算，对所述模拟纠正影像每个像元进行坐标变换获得其地理坐标(B,L)，通过坐标变换计算对应所述DEM中的行列，插值计算获得地面高程值H；

将所述原始立体影像及所述四至地理坐标输入所述GPU，利用步骤(2)中获得的所述线程块和线程网格参数调用所述CUDA核函数对每个所述模拟纠正影像像元通过左上角地理坐标和像元像方坐标进行坐标变换获得像元的地理坐标(B,L),并根据所述DEM的左上角坐标坐标变换在DEM中的行列(dem_r,dem_c)，使用双线性插值进行重采样计算，获得所述地面高程H；

(4)根据所述模拟纠正影像像元的地理坐标(B,L)和所述地面高程值H，通过所述RFM模型计算在所述原始立体影像中的行列，并以该行列进行插值计算的灰度值作为所述纠正影像的灰度值，获得正射纠正影像。

2.根据权利要求1所述的一种自适应线程参数的遥感正射影像并行生成方法，其特征在于，所述步骤(1)建立模拟纠正影像包括：

(1.1)根据所述原始立体影像的左上、右上、左下和右下四个角角点像方坐标(r₀,c₀)以及所述RPC文件的参数建立RFM模型；

(1.2)以所述DEM的平均高程H₀作为初值，通过所述RFM模型计算获得初始经纬度(B₁,L₁)，坐标变换得到(B₁,L₁)在所述DEM中的行列(x₁,y₁)，插值求得高程H₁；

(1.3)将所述像方坐标(r₀,c₀)和高程H₁输入所述RFM模型，获取新的经纬度(B₂,L₂)及高程H₂，重复迭代计算，直到计算的经纬度满足误差要求，将此经纬度作为角点的地理坐标，并将四个角点的地理坐标进行对比，获得影像的最大、最小维度和最大、最小经度作为四至，计算所述模拟纠正影像的行列数，建立所述模拟纠正影像。

3.根据权利要求2所述的一种自适应线程参数的遥感正射影像并行生成方法，其特征在于，所述RFM模型公式为：

其中：

Num_L(P,L,H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²

+a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³

+a₁₇PH²+a₁₈L²H+a₁₉P²H+a₂₀H³

Den_L(P,L,H)＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²

+b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³+b₁₇PH²+b₁₈L²H+b₁₉P²H+b₂₀H³

Num_S(P,L,H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²

+c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³

+c₁₇PH²+c₁₈L²H+c₁₉P²H+c₂₀H³

Den_S(P,L,H)＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²

+d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³+d₁₇PH²+d₁₈L²H+d₁₉P²H+d₂₀H³

(P,L,H)为标准化的地面坐标；(X,Y)为标准化的影像坐标，其标准化原理如下：

其中：LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的标准化参数；

SAMP_OFF、SAMP_SCALE、LINE_OFF、LINE_SCALE为影像坐标的标准化参数，上述参数均在所述RPC文件中获得。

4.根据权利要求1所述的一种自适应线程参数的遥感正射影像并行生成方法，其特征在于，所述步骤(4)中获得正射纠正影像包括：根据所述模拟纠正影像像元的地理坐标(B,L)和地面高程H，通过所述RFM模型计算得到所述模拟纠正影像像元对应在所述原始立体影像像方中的行列(r,c)，使用插值进行重采样计算，获得灰度值Y作为纠正影像的灰度值，每个像元的地理坐标为(B,L,H)，纠正影像的行列坐标为(r,c)，使用纠正后数据生成所述正射纠正影像。