Landsat 5 TM沙丘迁移场的误差消除与时序反演：以毛乌素沙地西北部沙丘为例

杨定江，冯光财，冯志雄，李国帅，张洁

1.中南大学 地球科学与信息物理学院 雷达遥感实验室,长沙 410083;

2.中国科学院青藏高原研究所 环境变化与地表过程重点实验室&国家青藏高原科学数据中心,北京 100101;

3.中国科学院青藏高原研究所 青藏高原地球系统科学国家重点实验室,北京 100101;

4.青海省基础测绘院,西宁 810001

1 引言

沙丘迁移是在风力作用下，迎风坡的沙粒不断被吹蚀，背风坡不断堆积的方式来实现的（Ding等，2020b；Wang 等，2008）。沙丘迁移使得沙漠向外扩张，可能会导致掩埋耕地、污染水源、形成沙尘暴甚至破坏基础设施等社会经济与生态环境问题（Bruno等，2018）。量化沙丘的迁移模式以有效地服务于防沙治沙工程是非常必要的。从技术发展上看，量化沙丘迁移的手段基本上可以分为4类：（1）使用传统仪器进行野外调查监测（Dong 等，2000；Ould Ahmedou等，2007）；（2）通过数字摄影测量和InSAR 等技术获取时序DEM 变化（Brown和Arbogast，1999；Hugenholtz 和Barchyn，2010；Dabboor等，2013；郭华东，1986；金亚秋，1997）；（3）使用ArcGIS 等软件对高分辨率的影像进行目视解译获取沙丘的脊线移动和范围变化（Tsoar和Blumberg，2002；Yao等，2007）；（4）通过光学影像匹配技术获取主从影像之间的沙丘迁移场（Bridges等，2012；邵梅 等，2020）。由于传统的调查监测费时费力、时序DEM 对数据质量要求高、目视解译矢量化内业繁重，而光学影像匹配技术有着精度高、自动化、大范围等特点已成为沙丘迁移的时序监测的优先选择（Ding等，2020a；Ali等，2020）。

COSI-Corr（Co-Registration of Optically Sensed Images and Correlation）是常用的一款光学影像匹配与互相关计算软件包（Leprince 等，2007a，2007b）。2008 年，首次应用到非洲乍得Bodélé 沙丘的监测后（Vermeesch 和Drake，2008），就广泛地应用于地球（Ding 等，2020a，2020b；邵梅 等，2020；Ali 等，2020）和火星（Bridges 等，2012，2017；Ayoub等，2014；Silvestro等，2016）的沙丘活动监测。Necsoiu等（2009）以亚北极大科伯克沙丘（Great Kobuk Sand Dunes）为例，测试出ASTER与SPOT 波段匹配的最佳线性组合。Bridges 等（2012）利用COSI-Corr 处理高分辨率的HiRISE 影像，获取了火星的沙丘迁移场，结果表明火星沙丘的移动规律与地球沙丘相似，但风成活动不如地球剧烈，沙丘移动也比较缓慢，这和火星的大气稀薄有关。Ding等（2016，2020a，2020b，2021）对Landsat 8 与Sentinel-2 迁移场误差进行分析和讨论后，提出了“短时空基线和辐射基线”构网的时序反演方法来提升精度，并对敦煌和民勤等地进行了监测和风沙灾害评估。Landsat 5 TM 数据集有着覆盖范围广、存档时间长、质量稳定、可免费获取等特点，且L1TP 级数据可满足COSI-Corr软件自动化处理要求（Baird等，2019），有利于进行沙丘的长时间序列的定量监测，但目前关于Landsat 5 TM 获取沙丘迁移场的误差源分析和时序反演研究仍然较少。

这里，本文获取了1991 年—2000 年覆盖毛乌素沙地西北部的37 景Landsat 5 TM 影像。首先，我们对影像进行了时空构网并设计实验探讨了Landsat 5 TM 迁移场中存在的误差类型和可能的产生原因。然后，根据误差分布特点，叙述了常用的相应误差消除方法和时序反演方法，并对改正前后的结果进行精度评估。最后，通过矢量场和时间序列分析了毛乌素沙地西北部沙丘迁移场的时空特征，并与以往的研究做比较，展望了光学影像匹配技术在沙丘迁移监测中的应用前景。

2 研究区域与实验设计

2.1 研究区域与Landsat 5数据介绍

毛乌素沙地位处中国的鄂尔多斯高原与黄土高原过渡区，范围约在37°27.5′N—39°22.5′N 和107°20′E—110°30′E 之间，面积约42200 km2。如图1所示，毛乌素沙地高程从西北（1400—1500 m）至东南（950—1100 m）逐渐降低，个别区域高程达1600 m。年均气温为6.0 ℃—8.5 ℃；年降雨在250—440 mm 之间，自东南至西北递减，主要集中在7—9月。冬春季节盛行西北风，夏秋季节盛行东南风（庞营军 等，2019），年均风速3 m/s。起沙风集中在冬末春季，其中西北方向的起沙风占54.5%，频率最高为3 月份；东南方向的起沙风占28.9%，频率最高为5月（哈斯 等，2006）。从沉积物来看，毛乌素沙地早在第四纪时期就已存在，几经变迁、扩张、固定和缩小才发育成现在的样子，其变迁演化和鄂尔多斯高原古地理环境的变化息息相关，而从人类时期又开始受到人为活动的剧烈影响（董光荣 等，1983）。由于人为活动频繁，加上生态环境脆弱，20 世纪以来该地沙漠化十分严重，风沙灾害频发，给周边居民造成了极大的生命威胁和财产损失。为此，它被列为1978 年中国“三北防护林计划”的重点治理区域之一，经过几十年的治理保护，近年来我们终于见证了“人进沙退”的治理过程。目前作为中国第一个将要成功治理的沙漠，量化其沙丘迁移的时空规律，对于其他沙漠的治理工作有着重要的参考价值和指导意义。

图1 毛乌素沙地的地理空间分布Fig.1 Geographical and spatial distribution of Mu Us Sandy Land

Landsat 5 卫星是美国航空航天局1984-03-01发射的Landsat 系列卫星之一，卫星轨道高度为705 km，倾斜角度为98.2°，重复周期为16 d。它配有Multispectral Scanner（MSS）传感器和Thematic Mapper（TM）传感器。MSS 传感器有四个光谱波段，地面采样间隔57×79 m；TM 传感器有7 个光谱波段，波段中心分别为：0.485 μm、0.56 μm、0.66 μm、0.83 μm、1.60 μm、11.45 μm、2.215 μm；反射波段1—5、7 的空间分辨率为30 m，热波段6的为120 m。研究区域的Landsat 5 基本信息如表1所示，共获取了1991 年—2000 年的37 景Landsat 5 TM Band 3影像。

表1 研究区域Landsat 5 TM的基本信息Table 1 Basic information of Landsat 5 TM in the study area

2.2 影像时空构网

获取沙丘的时序迁移首先需要进行影像配对。最为简单的配对策略是相邻时间影像的两两配对，能以最少观测数构建适用于沙丘监测的时间序列。但该方法由于没有多余观测，因此受异常值影响较大。假设有N+1景影像，时间顺序为[t0，t1，…，tN]，则主从影像之间的配对数最大为若全部使用这些可能的影像对组合，对误差处理和时序反演来说是一个巨大的挑战，不排除粗差值和异常值对结果精度可能会有较大的影响。为提高时序反演精度，我们提出了一种新的影像时空构网策略，以此平衡时序反演精度和计算效率。

光学影像匹配技术是通过检索主从影像间的纹理差异来获取地表的水平迁移场，而Landsat 5主从影像之间的太阳高度角差值SED（Sun Elevation Difference）、太阳方位角差值SAD（Sun Azimuth Difference）、时间基线TB（Temporal Baseline）等变量会影响配准后的迁移场精度（Ding等，2020a，2020b；Ali 等，2020），所以，为了削弱系统误差的大小和提高时序反演精度，需要对主从影像的这3个变量进行约束，获取一个稳健的主从影像构网。本文提出的构网策略如下：

（1）假设有N+1 景影像，将其按时间顺序排列为[t0，t1，…，tN]；影像序号为i(i=0，1，…，N)；

（2）对于第i景影像，和后面所有影像配对，计算影像对的TB、SED和SAD；

（3）设置影像对TB 为32—832 天、SED 小于15°、SAD 小于15°，搜索所有的满足条件的影像对；

（4）在步骤（3）获得的所有满足条件的影像对中，通过设置每个主影像最多与时间最接近的3 个从影像做相关性计算的方式进一步来提高计算效率；

（5）重复步骤（2）—（4），直到i=N，结束配对。绘制时空构网图。

研究区域构网的情况和参数见图2、表2 和表3，共构建了88对影像对。我们提出的影像对构网策略以3个变量约束迁移场中的系统误差，并将影像对数量限制在影像数量的3倍以下，平衡了系统误差和计算效率（以影像对数量为代用指标）。

表2 影像对的TB、SED和SADTable 2 TB、SED and SAD of image pair

表3 最终影像对构网的特征参数Table 3 Characteristic parameters of final image pair network

2.3 参数配置与数据处理

因为COSI-Corr 软件包对推扫式和摆扫式卫星数据的适用性较强，且已经广泛地应用于行星表面的沙丘监测研究（Ding等，2020a，2020b；Bridges等，2012，2017），其配准精度达到1/20—1/50 像素（Leprince等，2007a，2007b），又可以批量处理，所以本文将COSI-Corr 软件包作为数据处理平台。但不同卫星由于成像模式和影像空间分辨率差异导致配准和互相关的参数配置也不同（Ayoub 等，2009），如Ding等（2020b）对Landsat 8和Sentinel-2参数配置（最初窗口、最终窗口、移动步长、迭代次数和信噪比SNR（Signal to Noise Ratio）阈值）测试后发现，处理窗口与信噪比阈值会影响匹配精度，移动步长和迭代次数则对处理时长影响较大。由于后续的时序反演需求，本文对影像进行时空构网，获取的影像对数量仍然比较大，所以设置参数时需要平衡匹配精度和计算效率。这里，我们使用的是COSI-Corr 软件包的频率域互相关算法，经过测试获得的Landsat 5 影像对的最优参数配置如下：X和Y方向的最初搜索窗口为128像素，最终窗口为64 像素，两方向的处理步长为3 像素，稳健性迭代次数为2，信噪比阈值为0.9。此外，同一影像的不同波段对匹配结果也存在影响，一般会优先选择高空间分辨率的波段，比如Landsat 8 Pan 波段（全色波段）。若不同波段的空间分辨率一样时就需要测试最优的匹配波段，比如贺礼家等（2019）对Sentinel-2的4个10 m波段试验后发现Band 8 的匹配结果标准差最小，是提取大量级地震同震形变场的最佳波段。Landsat 5 波段进行测试后发现，Band 3（0.63—0.69 μm）可以有效记录沙丘表层的纹理信息，使得匹配结果稳健精确，所以本文选择Landsat 5 TM Band 3来获取沙丘迁移场。

本文的数据处理分为预处理、影像匹配、误差分析与处理、时序反演4 个模块，处理流程如图3所示。第1个模块是预处理，将研究区域所有影像裁剪至同一范围，然后据2.2 节策略进行影像时空构网设计，并输出主从影像列表和参数配置列表。第2 个模块是影像匹配，通过COSI-Corr 软件批量处理，输出每一影像对的东西向迁移场、南北向迁移场和信噪比（SNR）文件。第3个模块是误差分析与处理，分析沙丘迁移场的误差分布和大小，确定误差消除顺序和和消除方法。第4是时序反演，将选取稳定区域计算权阵，通过最小二乘法反演沙丘迁移的时间序列，进而拟合得到年平均速度场。

图3 数据处理流程Fig.3 Data processing flow

3 误差消除与时序反演

本文通过Landsat 5 TM 影像对的东西向和南北向的迁移场的分析统计发现，迁移场主要存在失相关噪声、轨道误差和姿态角误差，如图4 所示。这些误差遮掩着真正的沙丘迁移，时序反演之前应该进行消除或削弱。

图4 Landsat 5 TM初始迁移场的误差示意图（以1991-07-06—1991-08-23影像对为例；本文中所有迁移场向东和向北数值为正）Fig.4 Error diagram of Landsat 5 TM initial migration field（Take 1991-07-06—1991-08-23 image pair as an example；In this paper，the values of all migration fields eastward and northward are positive）

3.1 失相关噪声

失相关噪声与地物类型和地表辐射量相关。当地物变化过大或地表辐射属性异常，如辐射量过高、过低或者强烈的变化都会导致主从影像纹理不明显或缺失，在COSI-Corr 软件进行互相关匹配时无法得到最优结果或者出现匹配错误，产生失相关噪声。失相关噪声在SNR 文件中表现为SNR 值过低，在东西向和南北向的迁移场中表现为异常值，对于局部区域结果的可靠性和精确性影响较大。所以可以对其进行掩膜，掩膜区域有：（1）SNR文件中低于0.9的区域和NAN值区域，SNR文件和东西向南北向的迁移场是逐像元一一对应的，所以可以根据SNR 文件剔除迁移场中低SNR值区域；（2）异常值区域，对于东西向和南北向的迁移场，分别求其均值和中误差，超过均值±3 倍中误差的区域进行掩膜。

3.2 轨道误差

从美国地质勘探局USGS（United States Geological Survey）网站上下载的Landsat 5 TM 影像只做过简单的校正，再经过COSI-Corr 软件处理后迁移场中会存在一个系统性的偏移，这就是轨道误差。由于轨道误差量级较大，在后处理中，必须经过严格的几何模型进行消除。轨道误差为迁移场的主要误差源，本文用多项式拟合消除（Ding等，2016）。多项式拟合去除轨道误差在D-InSAR 处理中也极为常用，它是先对形变区域进行掩膜，接着使用多项式进行拟合得到误差场，然后再从原始迁移场中减去误差场进行消除。对于不同的卫星数据得到的迁移场的轨道误差大小和分布也有区别，Ding等（2016，2020b）用一次多项式消除Landsat 8 和Sentinel-2 的迁移场的轨道误差，但对于Landsat 5 TM 影像的迁移场，一次多项式一般是不能完全消除轨道误差的，有时需要用到二次或者三次才行，这可能是因为Landsat 5 的影像质量不如Landsat 8和Sentinel-2 造成的，毕竟Landsat 5 的服役时间最长，电子元件老化较为严重。

3.3 姿态角误差

宇宙环境并不稳定，太阳光压、温度变化和万有引力等或使卫星平台出现微小的抖动，导致卫星的姿态变化（俯仰、侧倾和偏航），这些不被记录的姿态变化再加上卫星影像不精确的稀疏采样，导致迁移场中沿卫星飞行方向分布着周期性的误差条纹，这就是卫星姿态角误差（Leprince 等，2007a，2007b）。与Landsat 8 相比，Landsat 5 TM不存在条带误差，但姿态角误差十分严重（Ding等，2016，2020b），对结果精度影响较大，必须进行消除。消除卫星姿态角误差常用的是传统的“均值消去法”（Scherler等，2008；Ding等，2016），即根据这些误差条纹的分布特点，沿卫星飞行方向的垂直向求取稳定区域的均值，然后用原始迁移场减去即可消除。但贺礼家等（2019）发现，Sentinel-2 影像对获取的地震形变场的卫星姿态角误差与传感器的12 个电荷耦合元件CCD（Chargecoupled Device）线阵列的排列有关，据此他们提出了改进的“均值消去法”消除姿态角误差，比之传统方法效果更好，因此本文也使用该方法。

3.4 时序反演

时序反演可以通过冗余观测数提高迁移场的精度和鲁棒性，并获取沙丘的累积迁移序列。本文在此使用的是较为常用的最小二乘法（Ali等，2020）。

假设有N+1 景影像，时间顺序为[t0，t1，…，tN]，经过构网和COSI-Corr 软件包匹配后获取了M对迁移场，记相邻时间的平均速度为VT=[v1，v2，…，vN]，则每一影像对的迁移场Δdi(i=1，2，…，M)的函数模型可以表示为

式中，tm(m=0，1，…，N-1)和ts(s=1，2，…，N)为迁移场Δdi的主从影像时刻；vj为tj和tj+1两时刻之间的平均速度。式（1）可以写为矩阵形式：

平差的随机模型为

式中，D为Δd的方差矩阵；为单位权方差，在此定为1；P为每一迁移场中稳定区域计算的权阵，在此我们使用的是方差—协方差函数定权（冯光财 等，2015）。当M＞N时，即可通过最小二乘法平差解算：

式中，B为设计矩阵，B(j，i)=ti-ti-1(tm+1 ≤ti-1，ti≤ts；j=1，2，…，M；i=1，2，…，N+1)，其他位置元素为0。有时设计矩阵B会出现秩亏的情况，导致无法平差，这时需要引入奇异值分解（Ding 等，2020b）。但在本文的2.2 节的基线构网方法获取的影像对时间上是有重叠冗余的，不存在秩亏情况，用最小二乘方法即可解算。

获取了相邻时间的平均速度后，设第一个时间形变量为0，进而获取累积形变序列L，最后可用最小二乘法拟合序列的年平均速度：

式中，W=[t1-t0，…，tN-tN-1]为相邻影像的时间间隔。

4 精度与结果分析

4.1 精度分析

根据Landsat 5 TM 影像对的误差特点和分布，我们确定了迁移场的精度改善流程和方法。如图5所示，原始迁移场的精度提升共分为4步：去失相关噪声（T1）、去轨道误差（T2）、去姿态角误差（T3）和时序反演（T4）。

图5 精度提升步骤与方法（以1991-07-06—1991-08-23结果为例，黑色多边形为精度评定区域）Fig.5 Steps and methods of improving accuracy（aking the results of 1991-07-06—1991-08-23 as an example，the black polygon is the accuracy evaluation area）

表4 和图6 为图5 的各步骤精度评估区域的结果统计情况。本文以标准差和均值作为误差消除和时序反演精度指标，结果如表4所示。对于东西向迁移场来说，经过T1—T3 处理后标准差相继降低2%、34%、20%，说明东西向的轨道误差＞姿态角误差＞失相关噪声。南北向迁移场经过T1—T3处理后标准差相继降低5%、23%、4%，说明南北向的轨道误差＞失相关噪声＞姿态角误差。而通过时序反演后，标准差又降低了13%和14%，表明时序反演能进一步提高迁移场精度。

表4 各处理步骤的精度评估Table 4 Accuracy evaluation of each processing step

图6 测量值概率分布（T0和T1结果差异不大，所以T0被T1结果遮掩了）Fig.6 Probability distribution of measured values（There is little difference between T0 and T1 results，so T0 is masked by T1 results）

通过表4和图6，能发现东西向的T0—T4标准差明显大于南北向，南北向的T0 和T1 均值的绝对值又远大于东西向，说明南北向误差量级较大但较为集中，东西向误差较小但分布较为稀疏。东西向经过各步骤处理后，均值逐渐趋于0，分布也更陡峭，其中轨道误差的消除最为明显。南北向各步骤后分布形状并无太大变化，但标准差逐渐减小，其中对均值改变最为明显的是轨道误差消除。

4.2 沙丘迁移的空间特征分析

结合研究区域的东西向和南北向的年平均速度可以得到沙丘迁移矢量，而沙丘迁移的玫瑰图则能显示沙丘迁移的空间特征。结果显示（图7），研究区域沙丘迁移和主风向（西北—东南向）基本一致，沙丘的主要移动方向为SE（约140°），次为WNN（约340°）。在研究区域中，63.9%的沙丘迁移速度为0—1 m/a，15.6%的迁移速度为1—1.5 m/a，8.2%的迁移速度为1.5—2 m/a，8.1%的迁移速度为2—3 m/a，4.2%的迁移速度则≥3 m/a。

图7 沙丘迁移的年平均速度玫瑰图（玫瑰手臂指向沙丘迁移的方向）Fig.7 Rose chart of annual average velocity of sand dune migration（Rose′s arm points in the direction of sand dune migration）

1991 年—2000 年，毛乌素沙地西北部的沙丘迁移速度显示出该研究区域沙丘的主要移动方向为SE，最大迁移速度超过6 m/a，移动方向与冬春盛行风向—西北风一致（图8）；虽然毛乌素沙地为双风向的主导区域，但是从图8不能推断出另一盛行风向—东南风。结合图7 和图9，这可能是由于西北风远强于东南风（庞营军 等，2019），长时间作用下向东南移动的沙丘运动掩盖了向西北移动的沙丘运动导致的。

图8 欧氏公式计算的年平均速度（（a）的底图为东西向与南北向合成速度大小，绿色圆点为图9中A点位置。（b）为（a）中蓝色边框的放大图）Fig.8 Annual average velocity calculated by European formula（The（a）bottom map of is the composite velocity of East-West and North-South directions，and the green dot is the position of point A in Fig.9.The（b）is an enlarged view of the blue border in（a））

4.3 沙丘迁移的时间特征分析

沙丘迁移的时间序列可以反映作用于沙丘表面的风力强度的时间变化。图9显示A 点的东西向年均速度为6.4 m/a，方向为正东；其南北向年均速度为-5.2 m/a，方向为正南；合成速度方向为东南向，随时间变化的迁移量呈现出较为明显的线性关系。由于双风向的作用，该研究区域的沙丘会沿西北—东南向往返运动，存在一个周期性波动，这样就会导致毛乌素沙地西北部沙丘横断面多是西北—东南向。5—6 月和10—11 月为西北风和东南风的交替月份，于是沙丘脊线会在西北—东南向摆动，迎风坡和背风坡也会相互转换（哈斯 等，2006）。在6月和11月左右沙丘迁移的方向会发生变化，但长时间尺度下沙丘整体上仍然向东南方向移动。

5 讨论

5.1 沙丘迁移的基本条件

强劲的风力和丰富的沙源物质是沙丘迁移的基本条件。沙丘的形成和发育本质上是沙粒在空间上的搬运，而只有达到有效风速（＞5 m/s）才能形成风沙流使得沙粒脱离地表（朱震达 等，1981）。所以沙丘的迁移模式由该区域的常年风况决定，而沙丘的形成和发育又依赖于有无充足的沙源补充。研究结果表明，毛乌素沙地在6月和11月左右西北风和东南风会相互转变（哈斯 等，2006），沙丘迁移方向也随之变化，所以该研究区域的沙丘会沿西北—东南方向往返运动（哈斯 等，2006；王静璞 等，2013），沙丘的常年排列方向也沿西北—东南方向排列。但西北向的起沙风频率和强度均大于东南向（庞营军 等，2019），所以该地沙丘向东南方向持续推进。风成沙乃是下伏沉积物经风力加工改造而成（朱震达 等，1981），而Liu和Yang（2018）的研究显示，毛乌素沙地西北部风成沙主要源于下伏砂岩和湖相沉积物，中沙和细沙含量较高，十分适合移动沙丘发育（哈斯 等，2006；李智佩 等，2006）。所以，毛乌素沙地西北部复杂的风况条件和丰富的沉积物分布，以及干旱少雨的气候、稀疏的植被和不合理的人类活动等因素决定了该区域常年普遍存在沙丘迁移的现象。

5.2 研究对比与新发现

本文获取的沙丘迁移结果大体上与其他研究的结果基本一致，但广阔的监测范围和精确的时序结果可以更详细的揭示风沙迁移特征，这是其他研究还不能做到的。图7 和图8 结果显示该研究区域普遍存在沙丘迁移活动，但占63.9%的区域迁移速度为0—1 m/a，≥3 m/a的区域只有4.2%，主要向东南向（约140°）移动，这说明该区域的风能为低能环境—输沙势在200 VU 以下（Fryberger，1979），且风况较为稳定。庞营军等（2019）根据毛乌素腹地2008 年—2011 年风况数据计算的年输沙势为66.75 VU，合成输沙势方向为146°，与本文结果较为一致。但许明静等（2020）监测毛乌素沙地西北部典型沙丘的平均沙丘迁移方向（122.5°）与我们的存在一定差异，这可能是由于该研究只矢量化了一部分沙丘造成的。而且本文所用技术的监测范围达几十公里至几百公里，时长可至几十年，结果连续且更为可靠，精度也完全满足监测需求，可以为防沙治沙工程的长期建设提供指导。

我们对结果进行分析发现：（1）毛乌素西北部区域的西北至东南方向沙丘速率呈降低趋势（图8）；（2）毛乌素的东西方向风力和南北方向风力在各地的分布并不均匀（图9），图7 的玫瑰图也表明该地沙丘主要迁移方向在110°—150°，而不是完全向145°方向迁移；（3）两方向的迁移量结果（图9）表明，东南季风无法改变沙丘迁移的趋势，但减缓了沙丘向东南方向侵蚀。因此，在防沙治沙时需要从东南方向往西北方向推进，先治理东南部的小沙、慢沙，后治理西北部的大沙、快沙；固沙时草方格排列方向也需要根据沙丘迁移方向（一般为当地风向）做细致的调整，不能严格按西北至东南方向放置；在夏秋季节盛行东南季风时，沙丘需要时间改造形态向反方向迁移，可以趁这时向东南方向的沙丘速率减缓而又未向西北方向迁移时进行草方格固沙，等待雨季来临时种植耐旱植被，提升固沙防沙的治理效果。而且当地相关部门可以根据本文的沙丘迁移矢量图规划农田防护林和防风林的种植范围和种植顺序，依据各地沙丘迁移速率和降雨、风速、土壤情况等使用不同的治沙手段和种植相应品种的植被，以达到提升治沙效率、减少植被死亡和节约治沙成本等效果。

所以，与许明静等（2020）和王静璞等（2013）通过Google Earth 历史影像矢量化典型沙丘的迁移运动相比，本文获取的各像素点的沙丘迁移时间序列（图9）量化了沙粒随沙丘表面风流运动的时空特征，且监测范围更大，监测时间更长，可以为防沙治沙工程实施等提供有力的技术支持和有益建议。

6 结论

综上所述，本文以1991 年—2000 年毛乌素沙地西北部37 景Landsat 5 TM 影像为基础，使用COSI-Corr 软件包获取了沙丘迁移的时间序列和对应的年平均速度，并分析了研究区域沙丘迁移的时空特征。为了限制初始迁移场的系统误差，我们提出了一种新的Landsat 5 TM 影像时空构网策略，即对影像对的SED、SAD 和TB 等3 个变量进行约束得到了一个稳健的主从影像构网，提高了计算效率和精度。通过研究分析，本文得出的主要结论为：（1）Landsat 5 TM 沙丘迁移场的误差源主要有轨道误差、姿态角误差和失相关噪声，消除后精度相继提升23%—34%、4%—20%、2%—5%，再利用最小二乘法反演后可提升13%—14%。（2）沙丘迁移的时间序列和年平均速度场可以较好的量化沙丘的迁移特征。（3）风况条件、沉积物分布及其他因素决定了毛乌素沙地西北部的沙丘迁移模式。（4）本文的监测结果大体上和其他研究基本一致，但本文的时序结果更加详细地揭示出沙粒随沙丘表面风流结构运动的时空特征，监测范围也更广阔，可以为防沙治沙工程、沙漠化评估和风沙灾害的定量监测等提供可靠的技术支持。

此外，我们提出的影像对构网策略、设计的误差处理步骤和使用的时序反演方法除稳健精确的量化沙丘迁移场外，也同样适用于冰川流速提取、滑坡监测等大范围、大量级、长时序的地表形变监测。

志 谢本文所用Landsat 5 TM数据源于USGS（glovis.usgs.gov/），文中各种图件主要由GMT5.3.1绘制，数据处理平台为采用美国加州理工学院研发的软件包COSI-Corr（http：//www.tectonics.caltech.edu/）。十分感谢几位审稿人和编辑提出的宝贵意见，这让我们许多没有思考到的地方得到补充，对文章的细节部分和逻辑结构上得到进一步的处理，使得该篇文章的学术水平和我们的写作水平得到提高；在此，我们表示由衷的感谢。