基于SAR影像强度信息的单体沙丘和沙丘场迁移速度测量方法

技术领域

本发明属于遥感和大地测量技术领域，具体涉及一种基于SAR影像强度信息的单体沙丘和沙丘场迁移速度测量方法，可以测量单体沙丘迁移速度和方向，进一步获得沙丘场迁移速度矢量场图。

背景技术

在风力作用下，迎风坡沙粒被吹蚀，背风沙粒堆积，实现了沙丘迁移。沙丘迁移使得沙漠向外扩张、掩埋耕地、形成沙尘暴甚至破坏基础设施等。多个单体沙丘和沙漠平原可组成沙丘场，沙丘场迁移是沙漠扩散和风沙危害的主要表现形式。因此，量化沙丘的迁移模式和沙丘场迁移的速度矢量场是防沙治沙工程的关键。传统人工野外实地测量面临工作量大、监测时间长、测量误差大的难题。

遥感测量技术可以大范围、周期性的测量地表形变信息。基于遥感影像测量沙丘场迁移速度主要有光学遥感测量和雷达遥感测量两种方法。光学遥感通过圈定沙脊线或沙丘轮廓线估算沙丘场迁移速度，目视解译受人为主观判定影响大且难以测量沙丘场迁移的微小变化。雷达遥感干涉测量技术依据SAR影像的相位信息可高精度的测量大范围地表形变，但单体沙丘和沙丘场迁移速度快、水平向扩散范围广，导致雷达遥感SAR影像的相干性损失大，沙丘迁移难以通过常规的干涉测量实现。因此，当前遥感测量方法不能满足单体沙丘迁移速度高精度测量、沙丘场大范围迁移速度测量的需求。

发明内容

针对现有技术的不足，本发明的目的是针对沙漠地区单体沙丘和沙丘场迁移速度快、扩散范围广、以水平向迁移运动为主的特征，开发了针对此特征的基于SAR影像强度信息的单体沙丘和沙丘场迁移速度测量方法，可测量单体沙丘迁移速度，进而获得大范围沙丘场的迁移速度矢量场图。

为实现以上目的，本发明通过以下技术方案予以实现：

一种基于SAR影像强度信息的单体沙丘和沙丘场迁移速度测量方法，该方法包括以下步骤：主从影像配准：根据影像时间顺序对日期临近的两幅影像配对并指定日期较早的影像为主影像，日期较晚的影像为从影像，对主从SAR影像进行配准，并将配准后的影像裁剪至研究区；单体沙丘的质心坐标获取：对SAR影像进行过采样处理和加权归一化处理，并根据像素灰度值的加权平均法计算单体沙丘的质心坐标；单体沙丘迁移速度求解：计算单体沙丘的质心偏移信息，并根据质心偏移信息计算单体沙丘的质心迁移速度及方向，生成单体沙丘迁移速度矢量图；研究区的地面控制点网格建立：在配准后的SAR影像选定一个地面控制点，通过归一化互相关法确定地面控制点在主从影像中的位置，将地面控制点连接，建立研究区的地面控制点网格；沙丘场迁移速度矢量场计算：以地面控制点网格为基础，计算控制点迁移速度矢量场，集合所有控制点迁移速度矢量场，获得沙丘场迁移速度矢量场图；单体沙丘迁移速度和沙丘场迁移速度矢量场获取：根据SAR影像获取的时间顺序，逐一配对两幅SAR影像，重复上述步骤，获得时间序列的单体沙丘迁移速度矢量图和沙丘场迁移速度矢量场图。

在一些实施例中，所述对主从SAR影像进行配准还包括：对主从SAR影像应用精密轨道数据降低卫星轨道误差，并借助开源数字高程模型(DEM)对主从影像进行配准。

在一些实施例中，所述对SAR影像进行过采样处理具体为：通过最近邻插值法对低分辨率影像进行重构，在不增加感知信息量的情况下提高图像的分辨率，所述最近邻插值法基于距离待求位置最近的四个领域位置上的像素点，从该四个像素点中选择与待求位置距离最近的像素点，并将该像素点的灰度值赋予待求位置的像素点。

在一些实施例中，所述对SAR影像进行加权归一化处理具体为：通过对SAR影像强度信息的像素值在红(R)、绿(G)、蓝(B)三个通道进行加权归一化处理，获得强度影像像素(i，j)的归一化灰度值(img(i，j))，将归一化灰度值用作加权平均的权重系数，计算沙丘的质心坐标(x，y)。

在一些实施例中，所述计算单体沙丘的质心坐标(x，y)通过公式(1)～(2)实现：

img(i，j)＝(0.30R(i，j)+0.59G(i，j)+0.11B(i，j))/255——(2)

公式(2)中0.30，0.59和0.11均为固定系数，为根据人眼对绿色敏感度最高、对蓝色敏感度最低确定；R(i，j)、G(i，j)、B(i.j)分别表示影像中点(i，j)的红、绿、蓝三个通道的像素值。

在一些实施例中，所述计算单体沙丘的质心偏移信息，并根据质心偏移信息计算单体沙丘的质心迁移速度及方向，具体为：计算沙丘在两幅影像时间间隔内的质心坐标沿雷达视线方向(LOS)的偏移距离和方向，根据雷达的入射角转换为单体沙丘迁移速度，迁移速度包括迁移速率和迁移方向。

在一些实施例中，根据单体沙丘质心在主从影像中的位置坐标(x

公式(3)～(5)中，θ表示雷达的入射角，对于固定卫星θ为固定值，单位为度，t表示两幅SAR影像获得的时间间隔，单位为天。

在一些实施例中，所述确定地面控制点在主从影像中的位置包括：

a.在两幅影像的控制点为中心设置一个搜寻窗口；

b.对搜寻窗口内像素的灰度值进行计算；

c.通过计算SAR影像中控制点的相关性建立控制点间一一对应的关系。

在一些实施例中，所述计算SAR影像中控制点的相关性具体为：以主影像(U

在一些实施例中，所述的进行归一化互相关值NCC(P

公式(6)～(8)中，I(x，y)、I(x

在一些实施例中，所述计算控制点迁移速度矢量场，获得沙丘场迁移速度矢量场图包括：

a.基于SAR影像强度信息计算控制点网格的像素偏移量；

b.通过局部平均法填补异常地面控制点的偏移量；

c.通过对网格上计算的速度进行插值，生成研究区的迁移速度矢量场图。

本发明相对于现有技术的有益效果是：本发明提出一种基于SAR影像强度信息的沙丘迁移速度测量方法，根据沙丘迁移前后两景SAR影像的强度信息，计算SAR影像中地面控制点在距离向和方位向上的迁移信息，计算单体沙丘迁移速度，获得沙丘场迁移的速度矢量场；本发明适用于雷达SAR影像覆盖的沙漠地区。本发明的具体优点至少包括如下一个或多个：

(1)本发明的测量方法，采用雷达遥感SAR影像的强度信息测量方法，能够测量迁移速度快、扩散范围广、以水平向迁移运动为主的单体沙丘和沙丘场迁移速度；

(2)本发明的测量方法，根据像素灰度值的加权平均法确定单体沙丘的质心坐标(x，y)，避免了目视解译确定沙丘轮廓线在计算中造成的误差；

(3)本发明的测量方法，基于SAR影像强度信息的测量方法同时适用于高精度单体沙丘迁移速度和大面积沙丘场迁移速度矢量场的计算；

(4)本发明的测量方法，通过SAR影像强度信息测量方法相较于传统遥感测量技术可以更好的克服天气、季节、地形起伏等因素造成的测量误差；

(5)本发明的测量方法，通过SAR影像的幅度值测量水平向的大变形，突破了合成孔径雷达干涉测量技术中的相位解缠难以测量大形变的局限。

应当理解，本发明任一实施方式的实现并不意味要同时具备或达到上述有益效果的多个或全部。

附图说明

为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引申获得其它的实施附图。

本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容涵盖的范围内。

图1示例性示出本发明一种实施例测量方法的流程图；

图2示例性示出某单体沙丘的遥感影像识别与迁移速度测量图；

图3A～3D示例性示出单体沙丘迁移图，图中蓝色、绿色、红色点分别表示该沙丘在2019、2020和2021年的质心坐标；

图4A～4C示例性示出沙丘场时间序列迁移速度场图，图中的箭头方向表示研究周期内的迁移方向，箭头颜色代表迁移距离，其中a～l分别为2020年1月至2021年1月的沙丘场时间序列迁移速度场图。

在各个附图中，相同或对应的标号表示相同或对应的部分。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明实施例作进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

应当理解，术语“包括/包含”、“由……组成”或者任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的产品、设备、过程或方法不仅包括那些要素，而且需要时还可以包括没有明确列出的其他要素，或者是还包括为这种产品、设备、过程或方法所固有的要素。在没有更多限制的情况下，由语句“包括/包含……”、“由……组成”限定的要素，并不排除在包括所述要素的产品、设备、过程或方法中还存在另外的相同要素。

还需要理解，术语“上”、“下”、“前”、“后”、“左”、“右”、“顶”、“底”、“内”、“外”等指示方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置、部件或结构必须具有特定的方位、以特定的方位构造或操作，不能理解为对本发明的限制。

此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。

本发明涉及一种基于SAR影像强度信息的单体沙丘和沙丘场迁移速度测量方法，根据沙丘迁移前后两景SAR影像的强度信息，计算SAR影像中地面控制点在距离向和方位向上的迁移信息，计算单体沙丘迁移速度，获得沙丘场迁移的速度矢量场。本发明适用于雷达SAR影像覆盖的沙漠地区。

如图1所示的流程图，该方法包括以下步骤：步骤1，主从影像配准：根据影像时间顺序对日期临近的两幅影像配对并指定日期较早的影像为主影像，日期较晚的影像为从影像，对主从SAR影像进行配准，并将配准后的影像裁剪至研究区；步骤2，单体沙丘的质心坐标获取：对SAR影像进行过采样处理和加权归一化处理，并根据像素灰度值的加权平均法计算单体沙丘的质心坐标；步骤3，单体沙丘迁移速度求解：计算单体沙丘的质心偏移信息，并根据质心偏移信息计算单体沙丘的质心迁移速度及方向，生成单体沙丘迁移速度矢量图；步骤4，研究区的地面控制点网格建立：在配准后的SAR影像选定一个地面控制点，通过归一化互相关法确定地面控制点在主从影像中的位置，将地面控制点连接，建立研究区的地面控制点网格；步骤5，沙丘场迁移速度矢量场计算：以地面控制点网格为基础，计算控制点迁移速度矢量场，集合所有控制点迁移速度矢量场，获得沙丘场迁移速度矢量场图；步骤6，单体沙丘迁移速度和沙丘场迁移速度矢量场获取：根据SAR影像获取的时间顺序，逐一配对两幅SAR影像，重复上述步骤，获得时间序列的单体沙丘迁移速度矢量图和沙丘场迁移速度矢量场图。

为了更好的理解上述技术方案，下面将结合附图以及具体的实施方式对上述技术方案进行详细的说明。

本发明实施例以新疆塔克拉玛干沙漠研究区为例，具体阐述本发明在实施例中的具体实施步骤。

步骤1：主从影像配准。

下载2020～2021年的开源Sentinel-1卫星数据及精密轨道数据，选定主从影像并进行初配准。根据影像时间顺序对日期临近的两幅影像逐一配对，指定日期较早的为主影像，日期较晚的为从影像。对主从SAR影像应用精密轨道数据，借助精密轨道文件提供的准确卫星位置和速度信息，更新SAR影像元数据的轨道状态向量，降低卫星轨道误差。

本实施例共使用遥感数据13景，时间跨度为2020年1月1日-2021年1月7日，共计362天，相邻影像间隔时间最大36天，最小28天。影像数据范围完全覆盖研究区数据，成像模式为渐进扫描模式(IW)，极化方式为双极化(VV、VH)，入射角为38.9°，影像数据类型为GRD数据，研究区Sentinel-1A卫星数据的信息如表1所示。

借助开源数字高程模型(DEM)数据弥补地形引起的偏移，进行影像精配准。在垂直基线较大的情况下，地面的起伏会在主从影像之间引起额外的偏移，高程每变化1000m就会出现约2m的局部偏移；DEM模型辅助的计算考虑了由地形产生的偏移量误差。下载分辨率为30m的开源SRTM 1Sec HGT数字高程模型(DEM)数据，进行精确配准，空间参考坐标系为GCS_WGS_1984。

表1 Sentinel-1A卫星数据信息说明

将配准后的影像裁剪至研究区，根据实施例中沙丘场的经纬度信息进行研究区裁剪，经纬度信息如下表：

表2实施例沙丘场经纬度

步骤2：获取单体沙丘的质心坐标。

a.对SAR影像进行过采样处理，提高影像精度。

由于沙丘通常具有较小的尺度和较高的空间变异性，常规的雷达遥感影像中往往难以准确地捕捉和表征沙丘区域的细节信息。因此，对影像进行过采样处理，利用最近邻插值法对低分辨率影像进行重构，在不增加感知信息量的情况下提高图像的分辨率，该方法考虑距离待求位置最近的四个领域位置上的像素点，从这四个像素点中选择与待求位置距离最近的像素点，并将该像素点的灰度值赋给待求位置的像素点。

b.单体沙丘质心坐标计算。

公式(1)对SAR影像强度信息的像素值在红(R)、绿(G)、蓝(B)三个通道进行归一化处理，获得像素(i，j)的归一化灰度值(img(i，j))。无植被覆盖的沙丘产生较强的后向散射，在雷达强度图像中显示出较小的灰度值。沙丘边界区域的后向散射强度逐渐递减，灰度值则逐渐变大。对影像的像素值进行归一化处理，基于像素灰度值的加权平均方法确定沙丘的质心坐标：

img(i，j)＝(0.30R(i，j)+0.59G(i，j)+0.11B(i，j))/255——(2)

将归一化系数用作加权平均的权重系数，上述公式(1)～(2)基于像素灰度值的加权平均法确定沙丘的质心坐标(x，y)，以避免对具体沙丘轮廓的计算造成误差。公式(2)中的0.30，0.59和0.11均为固定系数，它们是根据人眼对绿色敏感度最高、对蓝色敏感度最低确定的；R(i，j)、G(i，j)、B(i.j)分别表示影像中点(i，j)的红、绿、蓝三个通道的像素值。

步骤3：单体沙丘迁移速度求解。

计算沙丘在两幅影像时间间隔内的质心坐标沿雷达视线方向(LOS)的偏移距离和方向，根据雷达的入射角转换为沙丘场迁移速度，生成单体沙丘迁移速度矢量图，迁移速度包括迁移速率和迁移方向：

上述公式(3)～(5)计算出沙丘的迁移距离和迁移方向。其中，参数θ表示影像的入射角。在本实施例中，θ取值为38.9°。

单体沙丘的遥感影像识别与迁移速度测量图如图2所示，图2中左、中、右三图分别表示大范围至小范围的SAR影像，单体沙丘在最右侧的放大影像中可被清晰识别。

单体沙丘迁移速度和方向图如图3A√3D所示，选取了20个单体沙丘，本实施例中的沙丘场迁移速度与方向测量值如表3所示，单体沙丘平均迁移速度为4.9m/y，单体沙丘平均迁移方向211.12°，测量值与历史监测资料吻合度较高，说明本发明的方法测量结果准确。

表3实施例单体沙丘迁移距离、迁移速率与迁移方向

步骤4：研究区的地面控制点网格建立。

研究区的地面控制点网格建立可通过以下分步骤获得：

a.在配准后的SAR影像选定一个地面控制点网格；

结合研究区内沙丘场迁移情况的野外调查数据，选定一个地面控制点。

b.通过归一化互相关法确定地面控制点在主影像和从影像中的位置；

以两幅影像的控制点为中心设置一个搜寻窗口，在本实施例中搜索窗口为128×128，通过对窗口内像素的灰度值进行计算，建立控制点间一一对应的关系，该对应关系通过计算SAR影像内向量控制点的相关性表示。

以主影像(I

公式(6)～(8)中，I(x，y)、I(x

通过上述公式(6)～(8)计算影像块间的归一化互相关值，互相关值最大的位置为最佳匹配位置。基于SAR强度影像的像素灰度特征，设定主从影像中的搜寻窗口，窗口范围内灰度特征相关性最大的控制点对被匹配为对应的地面控制点。

c.将地面控制点连接，建立研究区的地面控制点网格。

本实施例中，网格范围与方位角间距设置为400×400米，总控制点网格点数为3648个。

步骤5：沙丘场迁移速度矢量场计算。

沙丘场迁移速度矢量场计算可通过以下分步骤获得：

a.基于SAR影像强度信息的像素偏移量计算：

控制点迁移速度矢量场包括地面控制点的坐标、迁移速度、迁移方向、偏移量、方位向偏移量和距离向偏移量。以步骤4计算得到的地面控制点网格为基础，根据互相关估计的偏移量计算选定地面控制点上的迁移速度。计算沙丘场地面控制点网格的像素偏移量信息，包括2020～2021年期间每个月沙丘场地面控制点的坐标、迁移速度、迁移方向、偏移量、方位向偏移量和距离向偏移量。通过对网格上计算的速度进行插值，生成研究区的迁移速度场图。绘制沙丘场每个月的迁移场图，以箭头表示迁移方向，箭头颜色表示迁移距离大小。

b.通过局部平均法填补异常地面控制点的偏移量：

调研研究区沙丘场迁移速度的研究资料，设置阈值为5m/day为地面控制点的最大偏移速度。匹配窗口的半径取决于地面控制点的最大偏移量，超出偏移量上限的地面控制点定义为无效点，从而提高主从影像中对同名控制点的配对成功率。对于超出偏移量上限的地面控制点，根据其周围的有效地面控制点的偏移量插值平均，以得到相应的偏移信息。

步骤6：获取时间序列的单体沙丘迁移速度和沙丘场迁移速度矢量场。

根据影像获取的时间顺序，逐一配对两幅影像，重复上述步骤1～5，获得时间序列的单体沙丘迁移速度图和沙丘场迁移速度矢量场。沙丘场时间序列迁移速度矢量场图如图4A～4C所示，图中的箭头方向表示研究周期内的迁移方向，箭头颜色则代表迁移距离，其中a～l分别为2020年1月至2021年1月的沙丘场时间序列迁移速度场图。本实施例中的沙丘场迁移速度与方向测量值如表4所示。12个月的沙丘场每个月的平均迁移距离为2.02米，平均迁移方向为162.23°，测量值与历史监测资料吻合度较高，说明本发明的方法测量结果准确。

表4实施例沙丘场迁移距离、迁移速率与迁移方向

虽然在上面论述中包含了若干具体实现细节，但是这些不应当被解释为对本发明的范围的限制。在单独的实施例的上下文中描述的某些特征还可以组合地实现在单个实现中。相反地，在单个实现的上下文中描述的各种特征也可以单独地或以任何合适的子组合的方式实现在多个实现中。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。