楊定江，馮光財，馮志雄，李國帥，張潔

1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院 雷達遙感實驗室,長沙 410083;

2.中國科學(xué)院青藏高原研究所 環(huán)境變化與地表過程重點實驗室&國家青藏高原科學(xué)數(shù)據(jù)中心,北京 100101;

3.中國科學(xué)院青藏高原研究所 青藏高原地球系統(tǒng)科學(xué)國家重點實驗室,北京 100101;

4.青海省基礎(chǔ)測繪院,西寧 810001

1 引言

沙丘遷移是在風(fēng)力作用下，迎風(fēng)坡的沙粒不斷被吹蝕，背風(fēng)坡不斷堆積的方式來實現(xiàn)的（Ding等，2020b；Wang 等，2008）。沙丘遷移使得沙漠向外擴張，可能會導(dǎo)致掩埋耕地、污染水源、形成沙塵暴甚至破壞基礎(chǔ)設(shè)施等社會經(jīng)濟與生態(tài)環(huán)境問題（Bruno等，2018）。量化沙丘的遷移模式以有效地服務(wù)于防沙治沙工程是非常必要的。從技術(shù)發(fā)展上看，量化沙丘遷移的手段基本上可以分為4類：（1）使用傳統(tǒng)儀器進行野外調(diào)查監(jiān)測（Dong 等，2000；Ould Ahmedou等，2007）；（2）通過數(shù)字?jǐn)z影測量和InSAR 等技術(shù)獲取時序DEM 變化（Brown和Arbogast，1999；Hugenholtz 和Barchyn，2010；Dabboor等，2013；郭華東，1986；金亞秋，1997）；（3）使用ArcGIS 等軟件對高分辨率的影像進行目視解譯獲取沙丘的脊線移動和范圍變化（Tsoar和Blumberg，2002；Yao等，2007）；（4）通過光學(xué)影像匹配技術(shù)獲取主從影像之間的沙丘遷移場（Bridges等，2012；邵梅 等，2020）。由于傳統(tǒng)的調(diào)查監(jiān)測費時費力、時序DEM 對數(shù)據(jù)質(zhì)量要求高、目視解譯矢量化內(nèi)業(yè)繁重，而光學(xué)影像匹配技術(shù)有著精度高、自動化、大范圍等特點已成為沙丘遷移的時序監(jiān)測的優(yōu)先選擇（Ding等，2020a；Ali等，2020）。

COSI-Corr（Co-Registration of Optically Sensed Images and Correlation）是常用的一款光學(xué)影像匹配與互相關(guān)計算軟件包（Leprince 等，2007a，2007b）。2008 年，首次應(yīng)用到非洲乍得Bodélé 沙丘的監(jiān)測后（Vermeesch 和Drake，2008），就廣泛地應(yīng)用于地球（Ding 等，2020a，2020b；邵梅 等，2020；Ali 等，2020）和火星（Bridges 等，2012，2017；Ayoub等，2014；Silvestro等，2016）的沙丘活動監(jiān)測。Necsoiu等（2009）以亞北極大科伯克沙丘（Great Kobuk Sand Dunes）為例，測試出ASTER與SPOT 波段匹配的最佳線性組合。Bridges 等（2012）利用COSI-Corr 處理高分辨率的HiRISE 影像，獲取了火星的沙丘遷移場，結(jié)果表明火星沙丘的移動規(guī)律與地球沙丘相似，但風(fēng)成活動不如地球劇烈，沙丘移動也比較緩慢，這和火星的大氣稀薄有關(guān)。Ding等（2016，2020a，2020b，2021）對Landsat 8 與Sentinel-2 遷移場誤差進行分析和討論后，提出了“短時空基線和輻射基線”構(gòu)網(wǎng)的時序反演方法來提升精度，并對敦煌和民勤等地進行了監(jiān)測和風(fēng)沙災(zāi)害評估。Landsat 5 TM 數(shù)據(jù)集有著覆蓋范圍廣、存檔時間長、質(zhì)量穩(wěn)定、可免費獲取等特點，且L1TP 級數(shù)據(jù)可滿足COSI-Corr軟件自動化處理要求（Baird等，2019），有利于進行沙丘的長時間序列的定量監(jiān)測，但目前關(guān)于Landsat 5 TM 獲取沙丘遷移場的誤差源分析和時序反演研究仍然較少。

這里，本文獲取了1991 年—2000 年覆蓋毛烏素沙地西北部的37 景Landsat 5 TM 影像。首先，我們對影像進行了時空構(gòu)網(wǎng)并設(shè)計實驗探討了Landsat 5 TM 遷移場中存在的誤差類型和可能的產(chǎn)生原因。然后，根據(jù)誤差分布特點，敘述了常用的相應(yīng)誤差消除方法和時序反演方法，并對改正前后的結(jié)果進行精度評估。最后，通過矢量場和時間序列分析了毛烏素沙地西北部沙丘遷移場的時空特征，并與以往的研究做比較，展望了光學(xué)影像匹配技術(shù)在沙丘遷移監(jiān)測中的應(yīng)用前景。

2 研究區(qū)域與實驗設(shè)計

2.1 研究區(qū)域與Landsat 5數(shù)據(jù)介紹

毛烏素沙地位處中國的鄂爾多斯高原與黃土高原過渡區(qū)，范圍約在37°27.5′N—39°22.5′N 和107°20′E—110°30′E 之間，面積約42200 km2。如圖1所示，毛烏素沙地高程從西北（1400—1500 m）至東南（950—1100 m）逐漸降低，個別區(qū)域高程達1600 m。年均氣溫為6.0 ℃—8.5 ℃；年降雨在250—440 mm 之間，自東南至西北遞減，主要集中在7—9月。冬春季節(jié)盛行西北風(fēng)，夏秋季節(jié)盛行東南風(fēng)（龐營軍 等，2019），年均風(fēng)速3 m/s。起沙風(fēng)集中在冬末春季，其中西北方向的起沙風(fēng)占54.5%，頻率最高為3 月份；東南方向的起沙風(fēng)占28.9%，頻率最高為5月（哈斯 等，2006）。從沉積物來看，毛烏素沙地早在第四紀(jì)時期就已存在，幾經(jīng)變遷、擴張、固定和縮小才發(fā)育成現(xiàn)在的樣子，其變遷演化和鄂爾多斯高原古地理環(huán)境的變化息息相關(guān)，而從人類時期又開始受到人為活動的劇烈影響（董光榮 等，1983）。由于人為活動頻繁，加上生態(tài)環(huán)境脆弱，20 世紀(jì)以來該地沙漠化十分嚴(yán)重，風(fēng)沙災(zāi)害頻發(fā)，給周邊居民造成了極大的生命威脅和財產(chǎn)損失。為此，它被列為1978 年中國“三北防護林計劃”的重點治理區(qū)域之一，經(jīng)過幾十年的治理保護，近年來我們終于見證了“人進沙退”的治理過程。目前作為中國第一個將要成功治理的沙漠，量化其沙丘遷移的時空規(guī)律，對于其他沙漠的治理工作有著重要的參考價值和指導(dǎo)意義。

圖1 毛烏素沙地的地理空間分布Fig.1 Geographical and spatial distribution of Mu Us Sandy Land

Landsat 5 衛(wèi)星是美國航空航天局1984-03-01發(fā)射的Landsat 系列衛(wèi)星之一，衛(wèi)星軌道高度為705 km，傾斜角度為98.2°，重復(fù)周期為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。研究區(qū)域的Landsat 5 基本信息如表1所示，共獲取了1991 年—2000 年的37 景Landsat 5 TM Band 3影像。

表1 研究區(qū)域Landsat 5 TM的基本信息Table 1 Basic information of Landsat 5 TM in the study area

2.2 影像時空構(gòu)網(wǎng)

獲取沙丘的時序遷移首先需要進行影像配對。最為簡單的配對策略是相鄰時間影像的兩兩配對，能以最少觀測數(shù)構(gòu)建適用于沙丘監(jiān)測的時間序列。但該方法由于沒有多余觀測，因此受異常值影響較大。假設(shè)有N+1景影像，時間順序為[t0，t1，…，tN]，則主從影像之間的配對數(shù)最大為若全部使用這些可能的影像對組合，對誤差處理和時序反演來說是一個巨大的挑戰(zhàn)，不排除粗差值和異常值對結(jié)果精度可能會有較大的影響。為提高時序反演精度，我們提出了一種新的影像時空構(gòu)網(wǎng)策略，以此平衡時序反演精度和計算效率。

光學(xué)影像匹配技術(shù)是通過檢索主從影像間的紋理差異來獲取地表的水平遷移場，而Landsat 5主從影像之間的太陽高度角差值SED（Sun Elevation Difference）、太陽方位角差值SAD（Sun Azimuth Difference）、時間基線TB（Temporal Baseline）等變量會影響配準(zhǔn)后的遷移場精度（Ding等，2020a，2020b；Ali 等，2020），所以，為了削弱系統(tǒng)誤差的大小和提高時序反演精度，需要對主從影像的這3個變量進行約束，獲取一個穩(wěn)健的主從影像構(gòu)網(wǎng)。本文提出的構(gòu)網(wǎng)策略如下：

（1）假設(shè)有N+1 景影像，將其按時間順序排列為[t0，t1，…，tN]；影像序號為i(i=0，1，…，N)；

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

（3）設(shè)置影像對TB 為32—832 天、SED 小于15°、SAD 小于15°，搜索所有的滿足條件的影像對；

（4）在步驟（3）獲得的所有滿足條件的影像對中，通過設(shè)置每個主影像最多與時間最接近的3 個從影像做相關(guān)性計算的方式進一步來提高計算效率；

（5）重復(fù)步驟（2）—（4），直到i=N，結(jié)束配對。繪制時空構(gòu)網(wǎng)圖。

研究區(qū)域構(gòu)網(wǎng)的情況和參數(shù)見圖2、表2 和表3，共構(gòu)建了88對影像對。我們提出的影像對構(gòu)網(wǎng)策略以3個變量約束遷移場中的系統(tǒng)誤差，并將影像對數(shù)量限制在影像數(shù)量的3倍以下，平衡了系統(tǒng)誤差和計算效率（以影像對數(shù)量為代用指標(biāo)）。

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

表3 最終影像對構(gòu)網(wǎng)的特征參數(shù)Table 3 Characteristic parameters of final image pair network

2.3 參數(shù)配置與數(shù)據(jù)處理

因為COSI-Corr 軟件包對推掃式和擺掃式衛(wèi)星數(shù)據(jù)的適用性較強，且已經(jīng)廣泛地應(yīng)用于行星表面的沙丘監(jiān)測研究（Ding等，2020a，2020b；Bridges等，2012，2017），其配準(zhǔn)精度達到1/20—1/50 像素（Leprince等，2007a，2007b），又可以批量處理，所以本文將COSI-Corr 軟件包作為數(shù)據(jù)處理平臺。但不同衛(wèi)星由于成像模式和影像空間分辨率差異導(dǎo)致配準(zhǔn)和互相關(guān)的參數(shù)配置也不同（Ayoub 等，2009），如Ding等（2020b）對Landsat 8和Sentinel-2參數(shù)配置（最初窗口、最終窗口、移動步長、迭代次數(shù)和信噪比SNR（Signal to Noise Ratio）閾值）測試后發(fā)現(xiàn)，處理窗口與信噪比閾值會影響匹配精度，移動步長和迭代次數(shù)則對處理時長影響較大。由于后續(xù)的時序反演需求，本文對影像進行時空構(gòu)網(wǎng)，獲取的影像對數(shù)量仍然比較大，所以設(shè)置參數(shù)時需要平衡匹配精度和計算效率。這里，我們使用的是COSI-Corr 軟件包的頻率域互相關(guān)算法，經(jīng)過測試獲得的Landsat 5 影像對的最優(yōu)參數(shù)配置如下：X和Y方向的最初搜索窗口為128像素，最終窗口為64 像素，兩方向的處理步長為3 像素，穩(wěn)健性迭代次數(shù)為2，信噪比閾值為0.9。此外，同一影像的不同波段對匹配結(jié)果也存在影響，一般會優(yōu)先選擇高空間分辨率的波段，比如Landsat 8 Pan 波段（全色波段）。若不同波段的空間分辨率一樣時就需要測試最優(yōu)的匹配波段，比如賀禮家等（2019）對Sentinel-2的4個10 m波段試驗后發(fā)現(xiàn)Band 8 的匹配結(jié)果標(biāo)準(zhǔn)差最小，是提取大量級地震同震形變場的最佳波段。Landsat 5 波段進行測試后發(fā)現(xiàn)，Band 3（0.63—0.69 μm）可以有效記錄沙丘表層的紋理信息，使得匹配結(jié)果穩(wěn)健精確，所以本文選擇Landsat 5 TM Band 3來獲取沙丘遷移場。

本文的數(shù)據(jù)處理分為預(yù)處理、影像匹配、誤差分析與處理、時序反演4 個模塊，處理流程如圖3所示。第1個模塊是預(yù)處理，將研究區(qū)域所有影像裁剪至同一范圍，然后據(jù)2.2 節(jié)策略進行影像時空構(gòu)網(wǎng)設(shè)計，并輸出主從影像列表和參數(shù)配置列表。第2 個模塊是影像匹配，通過COSI-Corr 軟件批量處理，輸出每一影像對的東西向遷移場、南北向遷移場和信噪比（SNR）文件。第3個模塊是誤差分析與處理，分析沙丘遷移場的誤差分布和大小，確定誤差消除順序和和消除方法。第4是時序反演，將選取穩(wěn)定區(qū)域計算權(quán)陣，通過最小二乘法反演沙丘遷移的時間序列，進而擬合得到年平均速度場。

圖3 數(shù)據(jù)處理流程Fig.3 Data processing flow

3 誤差消除與時序反演

本文通過Landsat 5 TM 影像對的東西向和南北向的遷移場的分析統(tǒng)計發(fā)現(xiàn)，遷移場主要存在失相關(guān)噪聲、軌道誤差和姿態(tài)角誤差，如圖4 所示。這些誤差遮掩著真正的沙丘遷移，時序反演之前應(yīng)該進行消除或削弱。

圖4 Landsat 5 TM初始遷移場的誤差示意圖（以1991-07-06—1991-08-23影像對為例；本文中所有遷移場向東和向北數(shù)值為正）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 失相關(guān)噪聲

失相關(guān)噪聲與地物類型和地表輻射量相關(guān)。當(dāng)?shù)匚镒兓^大或地表輻射屬性異常，如輻射量過高、過低或者強烈的變化都會導(dǎo)致主從影像紋理不明顯或缺失，在COSI-Corr 軟件進行互相關(guān)匹配時無法得到最優(yōu)結(jié)果或者出現(xiàn)匹配錯誤，產(chǎn)生失相關(guān)噪聲。失相關(guān)噪聲在SNR 文件中表現(xiàn)為SNR 值過低，在東西向和南北向的遷移場中表現(xiàn)為異常值，對于局部區(qū)域結(jié)果的可靠性和精確性影響較大。所以可以對其進行掩膜，掩膜區(qū)域有：（1）SNR文件中低于0.9的區(qū)域和NAN值區(qū)域，SNR文件和東西向南北向的遷移場是逐像元一一對應(yīng)的，所以可以根據(jù)SNR 文件剔除遷移場中低SNR值區(qū)域；（2）異常值區(qū)域，對于東西向和南北向的遷移場，分別求其均值和中誤差，超過均值±3 倍中誤差的區(qū)域進行掩膜。

3.2 軌道誤差

從美國地質(zhì)勘探局USGS（United States Geological Survey）網(wǎng)站上下載的Landsat 5 TM 影像只做過簡單的校正，再經(jīng)過COSI-Corr 軟件處理后遷移場中會存在一個系統(tǒng)性的偏移，這就是軌道誤差。由于軌道誤差量級較大，在后處理中，必須經(jīng)過嚴(yán)格的幾何模型進行消除。軌道誤差為遷移場的主要誤差源，本文用多項式擬合消除（Ding等，2016）。多項式擬合去除軌道誤差在D-InSAR 處理中也極為常用，它是先對形變區(qū)域進行掩膜，接著使用多項式進行擬合得到誤差場，然后再從原始遷移場中減去誤差場進行消除。對于不同的衛(wèi)星數(shù)據(jù)得到的遷移場的軌道誤差大小和分布也有區(qū)別，Ding等（2016，2020b）用一次多項式消除Landsat 8 和Sentinel-2 的遷移場的軌道誤差，但對于Landsat 5 TM 影像的遷移場，一次多項式一般是不能完全消除軌道誤差的，有時需要用到二次或者三次才行，這可能是因為Landsat 5 的影像質(zhì)量不如Landsat 8和Sentinel-2 造成的，畢竟Landsat 5 的服役時間最長，電子元件老化較為嚴(yán)重。

3.3 姿態(tài)角誤差

宇宙環(huán)境并不穩(wěn)定，太陽光壓、溫度變化和萬有引力等或使衛(wèi)星平臺出現(xiàn)微小的抖動，導(dǎo)致衛(wèi)星的姿態(tài)變化（俯仰、側(cè)傾和偏航），這些不被記錄的姿態(tài)變化再加上衛(wèi)星影像不精確的稀疏采樣，導(dǎo)致遷移場中沿衛(wèi)星飛行方向分布著周期性的誤差條紋，這就是衛(wèi)星姿態(tài)角誤差（Leprince 等，2007a，2007b）。與Landsat 8 相比，Landsat 5 TM不存在條帶誤差，但姿態(tài)角誤差十分嚴(yán)重（Ding等，2016，2020b），對結(jié)果精度影響較大，必須進行消除。消除衛(wèi)星姿態(tài)角誤差常用的是傳統(tǒng)的“均值消去法”（Scherler等，2008；Ding等，2016），即根據(jù)這些誤差條紋的分布特點，沿衛(wèi)星飛行方向的垂直向求取穩(wěn)定區(qū)域的均值，然后用原始遷移場減去即可消除。但賀禮家等（2019）發(fā)現(xiàn)，Sentinel-2 影像對獲取的地震形變場的衛(wèi)星姿態(tài)角誤差與傳感器的12 個電荷耦合元件CCD（Chargecoupled Device）線陣列的排列有關(guān)，據(jù)此他們提出了改進的“均值消去法”消除姿態(tài)角誤差，比之傳統(tǒng)方法效果更好，因此本文也使用該方法。

3.4 時序反演

時序反演可以通過冗余觀測數(shù)提高遷移場的精度和魯棒性，并獲取沙丘的累積遷移序列。本文在此使用的是較為常用的最小二乘法（Ali等，2020）。

假設(shè)有N+1 景影像，時間順序為[t0，t1，…，tN]，經(jīng)過構(gòu)網(wǎng)和COSI-Corr 軟件包匹配后獲取了M對遷移場，記相鄰時間的平均速度為VT=[v1，v2，…，vN]，則每一影像對的遷移場Δdi(i=1，2，…，M)的函數(shù)模型可以表示為

式中，tm(m=0，1，…，N-1)和ts(s=1，2，…，N)為遷移場Δdi的主從影像時刻；vj為tj和tj+1兩時刻之間的平均速度。式（1）可以寫為矩陣形式：

平差的隨機模型為

式中，D為Δd的方差矩陣；為單位權(quán)方差，在此定為1；P為每一遷移場中穩(wěn)定區(qū)域計算的權(quán)陣，在此我們使用的是方差—協(xié)方差函數(shù)定權(quán)（馮光財 等，2015）。當(dāng)M＞N時，即可通過最小二乘法平差解算：

式中，B為設(shè)計矩陣，B(j，i)=ti-ti-1(tm+1 ≤ti-1，ti≤ts；j=1，2，…，M；i=1，2，…，N+1)，其他位置元素為0。有時設(shè)計矩陣B會出現(xiàn)秩虧的情況，導(dǎo)致無法平差，這時需要引入奇異值分解（Ding 等，2020b）。但在本文的2.2 節(jié)的基線構(gòu)網(wǎng)方法獲取的影像對時間上是有重疊冗余的，不存在秩虧情況，用最小二乘方法即可解算。

獲取了相鄰時間的平均速度后，設(shè)第一個時間形變量為0，進而獲取累積形變序列L，最后可用最小二乘法擬合序列的年平均速度：

式中，W=[t1-t0，…，tN-tN-1]為相鄰影像的時間間隔。

4 精度與結(jié)果分析

4.1 精度分析

根據(jù)Landsat 5 TM 影像對的誤差特點和分布，我們確定了遷移場的精度改善流程和方法。如圖5所示，原始遷移場的精度提升共分為4步：去失相關(guān)噪聲（T1）、去軌道誤差（T2）、去姿態(tài)角誤差（T3）和時序反演（T4）。

圖5 精度提升步驟與方法（以1991-07-06—1991-08-23結(jié)果為例，黑色多邊形為精度評定區(qū)域）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 的各步驟精度評估區(qū)域的結(jié)果統(tǒng)計情況。本文以標(biāo)準(zhǔn)差和均值作為誤差消除和時序反演精度指標(biāo)，結(jié)果如表4所示。對于東西向遷移場來說，經(jīng)過T1—T3 處理后標(biāo)準(zhǔn)差相繼降低2%、34%、20%，說明東西向的軌道誤差＞姿態(tài)角誤差＞失相關(guān)噪聲。南北向遷移場經(jīng)過T1—T3處理后標(biāo)準(zhǔn)差相繼降低5%、23%、4%，說明南北向的軌道誤差＞失相關(guān)噪聲＞姿態(tài)角誤差。而通過時序反演后，標(biāo)準(zhǔn)差又降低了13%和14%，表明時序反演能進一步提高遷移場精度。

表4 各處理步驟的精度評估Table 4 Accuracy evaluation of each processing step

圖6 測量值概率分布（T0和T1結(jié)果差異不大，所以T0被T1結(jié)果遮掩了）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，能發(fā)現(xiàn)東西向的T0—T4標(biāo)準(zhǔn)差明顯大于南北向，南北向的T0 和T1 均值的絕對值又遠大于東西向，說明南北向誤差量級較大但較為集中，東西向誤差較小但分布較為稀疏。東西向經(jīng)過各步驟處理后，均值逐漸趨于0，分布也更陡峭，其中軌道誤差的消除最為明顯。南北向各步驟后分布形狀并無太大變化，但標(biāo)準(zhǔn)差逐漸減小，其中對均值改變最為明顯的是軌道誤差消除。

4.2 沙丘遷移的空間特征分析

結(jié)合研究區(qū)域的東西向和南北向的年平均速度可以得到沙丘遷移矢量，而沙丘遷移的玫瑰圖則能顯示沙丘遷移的空間特征。結(jié)果顯示（圖7），研究區(qū)域沙丘遷移和主風(fēng)向（西北—東南向）基本一致，沙丘的主要移動方向為SE（約140°），次為WNN（約340°）。在研究區(qū)域中，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 年，毛烏素沙地西北部的沙丘遷移速度顯示出該研究區(qū)域沙丘的主要移動方向為SE，最大遷移速度超過6 m/a，移動方向與冬春盛行風(fēng)向—西北風(fēng)一致（圖8）；雖然毛烏素沙地為雙風(fēng)向的主導(dǎo)區(qū)域，但是從圖8不能推斷出另一盛行風(fēng)向—東南風(fēng)。結(jié)合圖7 和圖9，這可能是由于西北風(fēng)遠強于東南風(fēng)（龐營軍 等，2019），長時間作用下向東南移動的沙丘運動掩蓋了向西北移動的沙丘運動導(dǎo)致的。

圖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 沙丘遷移的時間特征分析

沙丘遷移的時間序列可以反映作用于沙丘表面的風(fēng)力強度的時間變化。圖9顯示A 點的東西向年均速度為6.4 m/a，方向為正東；其南北向年均速度為-5.2 m/a，方向為正南；合成速度方向為東南向，隨時間變化的遷移量呈現(xiàn)出較為明顯的線性關(guān)系。由于雙風(fēng)向的作用，該研究區(qū)域的沙丘會沿西北—東南向往返運動，存在一個周期性波動，這樣就會導(dǎo)致毛烏素沙地西北部沙丘橫斷面多是西北—東南向。5—6 月和10—11 月為西北風(fēng)和東南風(fēng)的交替月份，于是沙丘脊線會在西北—東南向擺動，迎風(fēng)坡和背風(fēng)坡也會相互轉(zhuǎn)換（哈斯 等，2006）。在6月和11月左右沙丘遷移的方向會發(fā)生變化，但長時間尺度下沙丘整體上仍然向東南方向移動。

5 討論

5.1 沙丘遷移的基本條件

強勁的風(fēng)力和豐富的沙源物質(zhì)是沙丘遷移的基本條件。沙丘的形成和發(fā)育本質(zhì)上是沙粒在空間上的搬運，而只有達到有效風(fēng)速（＞5 m/s）才能形成風(fēng)沙流使得沙粒脫離地表（朱震達 等，1981）。所以沙丘的遷移模式由該區(qū)域的常年風(fēng)況決定，而沙丘的形成和發(fā)育又依賴于有無充足的沙源補充。研究結(jié)果表明，毛烏素沙地在6月和11月左右西北風(fēng)和東南風(fēng)會相互轉(zhuǎn)變（哈斯 等，2006），沙丘遷移方向也隨之變化，所以該研究區(qū)域的沙丘會沿西北—東南方向往返運動（哈斯 等，2006；王靜璞 等，2013），沙丘的常年排列方向也沿西北—東南方向排列。但西北向的起沙風(fēng)頻率和強度均大于東南向（龐營軍 等，2019），所以該地沙丘向東南方向持續(xù)推進。風(fēng)成沙乃是下伏沉積物經(jīng)風(fēng)力加工改造而成（朱震達 等，1981），而Liu和Yang（2018）的研究顯示，毛烏素沙地西北部風(fēng)成沙主要源于下伏砂巖和湖相沉積物，中沙和細(xì)沙含量較高，十分適合移動沙丘發(fā)育（哈斯 等，2006；李智佩 等，2006）。所以，毛烏素沙地西北部復(fù)雜的風(fēng)況條件和豐富的沉積物分布，以及干旱少雨的氣候、稀疏的植被和不合理的人類活動等因素決定了該區(qū)域常年普遍存在沙丘遷移的現(xiàn)象。

5.2 研究對比與新發(fā)現(xiàn)

本文獲取的沙丘遷移結(jié)果大體上與其他研究的結(jié)果基本一致，但廣闊的監(jiān)測范圍和精確的時序結(jié)果可以更詳細(xì)的揭示風(fēng)沙遷移特征，這是其他研究還不能做到的。圖7 和圖8 結(jié)果顯示該研究區(qū)域普遍存在沙丘遷移活動，但占63.9%的區(qū)域遷移速度為0—1 m/a，≥3 m/a的區(qū)域只有4.2%，主要向東南向（約140°）移動，這說明該區(qū)域的風(fēng)能為低能環(huán)境—輸沙勢在200 VU 以下（Fryberger，1979），且風(fēng)況較為穩(wěn)定。龐營軍等（2019）根據(jù)毛烏素腹地2008 年—2011 年風(fēng)況數(shù)據(jù)計算的年輸沙勢為66.75 VU，合成輸沙勢方向為146°，與本文結(jié)果較為一致。但許明靜等（2020）監(jiān)測毛烏素沙地西北部典型沙丘的平均沙丘遷移方向（122.5°）與我們的存在一定差異，這可能是由于該研究只矢量化了一部分沙丘造成的。而且本文所用技術(shù)的監(jiān)測范圍達幾十公里至幾百公里，時長可至幾十年，結(jié)果連續(xù)且更為可靠，精度也完全滿足監(jiān)測需求，可以為防沙治沙工程的長期建設(shè)提供指導(dǎo)。

我們對結(jié)果進行分析發(fā)現(xiàn)：（1）毛烏素西北部區(qū)域的西北至東南方向沙丘速率呈降低趨勢（圖8）；（2）毛烏素的東西方向風(fēng)力和南北方向風(fēng)力在各地的分布并不均勻（圖9），圖7 的玫瑰圖也表明該地沙丘主要遷移方向在110°—150°，而不是完全向145°方向遷移；（3）兩方向的遷移量結(jié)果（圖9）表明，東南季風(fēng)無法改變沙丘遷移的趨勢，但減緩了沙丘向東南方向侵蝕。因此，在防沙治沙時需要從東南方向往西北方向推進，先治理東南部的小沙、慢沙，后治理西北部的大沙、快沙；固沙時草方格排列方向也需要根據(jù)沙丘遷移方向（一般為當(dāng)?shù)仫L(fēng)向）做細(xì)致的調(diào)整，不能嚴(yán)格按西北至東南方向放置；在夏秋季節(jié)盛行東南季風(fēng)時，沙丘需要時間改造形態(tài)向反方向遷移，可以趁這時向東南方向的沙丘速率減緩而又未向西北方向遷移時進行草方格固沙，等待雨季來臨時種植耐旱植被，提升固沙防沙的治理效果。而且當(dāng)?shù)叵嚓P(guān)部門可以根據(jù)本文的沙丘遷移矢量圖規(guī)劃農(nóng)田防護林和防風(fēng)林的種植范圍和種植順序，依據(jù)各地沙丘遷移速率和降雨、風(fēng)速、土壤情況等使用不同的治沙手段和種植相應(yīng)品種的植被，以達到提升治沙效率、減少植被死亡和節(jié)約治沙成本等效果。

所以，與許明靜等（2020）和王靜璞等（2013）通過Google Earth 歷史影像矢量化典型沙丘的遷移運動相比，本文獲取的各像素點的沙丘遷移時間序列（圖9）量化了沙粒隨沙丘表面風(fēng)流運動的時空特征，且監(jiān)測范圍更大，監(jiān)測時間更長，可以為防沙治沙工程實施等提供有力的技術(shù)支持和有益建議。

6 結(jié)論

綜上所述，本文以1991 年—2000 年毛烏素沙地西北部37 景Landsat 5 TM 影像為基礎(chǔ)，使用COSI-Corr 軟件包獲取了沙丘遷移的時間序列和對應(yīng)的年平均速度，并分析了研究區(qū)域沙丘遷移的時空特征。為了限制初始遷移場的系統(tǒng)誤差，我們提出了一種新的Landsat 5 TM 影像時空構(gòu)網(wǎng)策略，即對影像對的SED、SAD 和TB 等3 個變量進行約束得到了一個穩(wěn)健的主從影像構(gòu)網(wǎng)，提高了計算效率和精度。通過研究分析，本文得出的主要結(jié)論為：（1）Landsat 5 TM 沙丘遷移場的誤差源主要有軌道誤差、姿態(tài)角誤差和失相關(guān)噪聲，消除后精度相繼提升23%—34%、4%—20%、2%—5%，再利用最小二乘法反演后可提升13%—14%。（2）沙丘遷移的時間序列和年平均速度場可以較好的量化沙丘的遷移特征。（3）風(fēng)況條件、沉積物分布及其他因素決定了毛烏素沙地西北部的沙丘遷移模式。（4）本文的監(jiān)測結(jié)果大體上和其他研究基本一致，但本文的時序結(jié)果更加詳細(xì)地揭示出沙粒隨沙丘表面風(fēng)流結(jié)構(gòu)運動的時空特征，監(jiān)測范圍也更廣闊，可以為防沙治沙工程、沙漠化評估和風(fēng)沙災(zāi)害的定量監(jiān)測等提供可靠的技術(shù)支持。

此外，我們提出的影像對構(gòu)網(wǎng)策略、設(shè)計的誤差處理步驟和使用的時序反演方法除穩(wěn)健精確的量化沙丘遷移場外，也同樣適用于冰川流速提取、滑坡監(jiān)測等大范圍、大量級、長時序的地表形變監(jiān)測。

志 謝本文所用Landsat 5 TM數(shù)據(jù)源于USGS（glovis.usgs.gov/），文中各種圖件主要由GMT5.3.1繪制，數(shù)據(jù)處理平臺為采用美國加州理工學(xué)院研發(fā)的軟件包COSI-Corr（http：//www.tectonics.caltech.edu/）。十分感謝幾位審稿人和編輯提出的寶貴意見，這讓我們許多沒有思考到的地方得到補充，對文章的細(xì)節(jié)部分和邏輯結(jié)構(gòu)上得到進一步的處理，使得該篇文章的學(xué)術(shù)水平和我們的寫作水平得到提高；在此，我們表示由衷的感謝。