基于SAR影像的干旱区冲/洪积扇地貌面定量分期研究-以河西走廊西部沙漠区的疏勒河洪积扇为例

杨勇忠， 任俊杰， 李东臣

1.中国科学院大学应急管理科学与工程学院，北京 100049；

2.应急管理部国家自然灾害防治研究院，北京 100085

0 引言

河流地貌是陆表系统最广泛的地貌类型之一，废弃的河流地貌为过去构造活动、气候变化和景观演变过程提供了重要的记录载体（Bull，1977，1991；Lubetkin and Clark，1988；徐锡伟等，2003；Matmon et al.，2005；Guralnik et al.，2010；黄飞鹏等，2021）。洪积扇是由山地河流带来的大量砾石和泥沙在山前堆积形成的扇形地貌，是河流系统中短暂存储沉积物的基本地貌单元（Huggett et al.，2016），洪积扇地貌单元的期次划分对于理解地貌演化、构造隆升和地表过程研究至关重要。

常用的洪积扇地貌面分期方法主要有3类：定性分析方法（郑荣章，2005；Regmi et al.，2014）、高分辨率地形数据法（Frankel and Dolan，2007；韩龙飞等，2019）和 遥 感 影 像 法（D’Arcy et al.，2018；苏 强 等，2020）。定性方法易受研究者经验和主观认识的影响，为地貌面的划分带来了诸多不确定性；高分辨率的数字高程模型（DEM）对数据的质量要求高，大范围数据获取的费用昂贵（Frankel and Dolan，2007；Fan and Atkinson，2018；Wu et al.，2018）；光学影像的参数易受到植被、大气条件等外界因素的影响，在定量描述地貌面时存在诸多不确定性。在大尺度干旱—半干旱地区洪积扇地貌面期次划分中，微波遥感方法由于其全天候、穿透性强、不易受环境变化影响等特点，成为近年来最理想的分期工具（Hetz et al.，2016；Sadeh et al.，2018）。

早在20世纪70年代，合成孔径雷达（SAR）技术就被用于地表参数（地表含水量和地表粗糙度）的测量（Ulaby et al.，1979）。近年来，遥感技术的蓬勃发展为干旱沙漠地形中的地貌地层冲积单元的测量提供了更好的数据源。利用RADARSAT-2影像数据，林国青等（2013）提出了更加适合干旱区地表粗糙度反演的新方法。Zhang and Guo（2013）通过研究SAR特征参数与地貌面年龄的直接相关性，表明C波段SAR数据（RADARSAT-2）在干旱地区地貌面分期中具有潜力。最新研究表明，L波段SAR数据对地表粗糙度最敏感，更适用于洪积扇地貌面分期（苏强等，2020；Ayari et al.，2021）。相关学者利用地貌面年龄、地表粗糙度和后向散射系数之间的相关关系，使用ALOS-SAR卫星数据完成了对地貌面的定量分期，并建立了后向散射系数值（R）和原位年龄（T）之间的稳定的幂律关系（Hetz et al.，2016；Su et al.，2023）。然而，由于大气条件的变化，同一地区不同时相的SAR数据反演的后向散射系数存在明显差异，使得反演结果容易受到数据质量的影响。以往的研究往往根据定性分析选择较好的数据源，并直接使用后向散射系数作为地表粗糙度的替代参数，缺乏对数据选取的评价方法。

文章以地貌面发育较好的疏勒河洪积扇为研究区，该洪积扇地形起伏小，各级洪积扇单元尺度大，全年气候干燥，植被稀少，适用于使用遥感方法开展洪积扇分期的研究。郑荣章（2005）通过卫星影像解译和野外调查认为该区域主要存在2期地貌面；Zhang and Guo（2013）利 用C波段RADARSAT-2雷达影像对洪积扇的部分区域进行了定量划分，共划分为4期地貌面。在已有研究成果的基础上，笔者尝试提出一种评估SAR数据源质量的方法，建立一套使用SAR数据进行洪积扇地貌单元定量分期的方案。该研究为更好地利用微波遥感数据开展地貌面分期提供了新思路。

1 研究区及影像概况

1.1 区域概况

疏勒河洪积扇位于中国西北部河西走廊中西部沙漠地区，南邻祁连山，北接北山，扇体总面积约为2100 km2（图1；王萍，2003；Guo et al.，2017）。疏勒河洪积扇的地貌面主要由第四纪洪积物组成，以晚更新世的戈壁砾石地层为主，全新世砂砾沉积也广泛分布。疏勒河发源于祁连山，向北流经研究区，区域内第四纪活动断层十分发育，南面山前为规模巨大的左旋走滑的阿尔金断裂（王萍，2003；王萍等，2004；Zhang and Guo，2013；云龙等，2021）。研究区的年平均气温为6.95～9.42 ℃，1月份温度最低，7月温度最高，年均降水量为39.6～63.4 mm，雨季发生在6—8月；年均蒸发量为2469～2869 mm，蒸发强度最高的时段在5—9月（毛洪亮等，2007；Guo et al.，2017）。

图1 疏勒河洪积扇地区地质构造图（引自《1︰50万中国地质图》公开版（http://www.ngac.org.cn）；断裂分布据Zelenin et al.，2022修改）Fig.1 Tectonic map of the diluvial fan area of Shule River (The geological map is quoted from the public version of the 1:500000 Geological Map of China at http://www.ngac.org.cn; Fracture distribution is modified after Zelenin et al.,2022)

1.2 影像数据

SAR系统存在入射波波长不同的工作波段——X波段（3.1 cm）、C波段（5.6 cm）、L波段（23.5 cm），以及4种常用的极化方式（HH、HV、VV、VH）。波段越长，穿透性越好；而不同的极化方式对不同目标的敏感性存在差异（Freeman，1992）。以往研究表明，当入射波为L波段、HH极化并且入射角大于30°时，雷达的后向散射系数强度值对地表粗糙度的差 异 性 最 为 敏 感（Beaudoin et al.，1990；Coppo et al.，1995；Ayari et al.，2021）。依据微波的物理特性和研究需求，研究选用ALOS PALSAR影像的Level 1.1 HH极化数据作为河流地貌分期研究影像数据（表1）。数据在NASA下属的数据开发平台EarthData（https://disc.gsfc.nasa.gov）上获取，可获取研究区域附近共计10个时相的相同分辨率数据，时间跨度为2007年7月—2009年10月。每个时相选取位于同轨道的两景影像，覆盖了疏勒河洪积扇地区的大部分区域。

表1 选取的SAR影像数据信息Table 1 Information of selected SAR images

2 地貌面分期原理

2.1 洪积扇地貌面分期原理

洪积扇是在干旱—半干旱地区暂时性山地水流出山口堆积形成的扇形地貌，广泛发育于河流出山口地带。由于复杂的气候变化和构造作用，河流经历多期堆积和下切过程，在山口地带表现为海拔高度不同的洪积扇地貌面（Blair and McPherson，1994）。地貌面形成后，会受到后期的侵蚀和风化作用，地貌面形成时代越老，遭受侵蚀的时间越长，地貌面上的砾石粒径逐渐变小，地貌面的粗糙度随之 减 小（图2；Ulaby et al.，1979；Evans et al.，1992；Sadeh et al.，2018）。基于这一自然过程，可以通过地表粗糙度来区分不同期次的洪积扇地貌单元。

图2 洪积扇演化示意图（据Blair and McPherson，1994修改）Fig.2 Schematic diagram of the evolution of alluvial fans (modified after Blair and McPherson,1994)

雷达后向散射系数是反映辐射能量与入射能量比值的参数，而粗糙的地表会导致雷达入射波的多次反射，形成漫反射场（图3）。随着地表粗糙度的增加，SAR接收到的回波能量增强。由此，可通过后向散射系数反演地表粗糙度（Greeley et al.，1988；Mattia et al.，1997；Zribi et al.，2005；李 森，2007）。基于地貌面年龄、粗糙度和后向散射系数三者之间的相关关系，可近似地将雷达后向散射系数作为地貌面粗糙度的替代参数，进而对地貌面进行分期（Hetz et al.，2016）。

图3 雷达波后向散射模式图Fig.3 Radar wave backscatter pattern diagram

2.2 SAR后向散射系数反演

使用ENVI SARscape软件对获取的Level 1.1单视复数数据进行处理后，得到SAR后向散射系数。首先，需要导入影像并进行多视处理和滤波；然后进行辐射定标和地理编码，以实现后向散射系数反演（Holecz et al.，2000）。在多视处理中，根据计算选择距离向视数为1、方位向视数为4，这样可以使方位向分辨率与距离向分辨率接近，以达到抑制噪声的目的。经过测试，使用Gamma滤波和Frost滤波的组合可以有效地抑制斑点噪声。

辐射定标是将SAR影像的亮度灰度值转换为绝对的辐射亮度，定标过程如公式（1）所示：

式中，DN为SAR影像预处理后像元亮度值，CF为SAR影像后向散射系数校正值；σ0为辐射亮度。

2.3 反演结果统计分析

首先，需要获取各期地貌面统计特征值作为原始样本。参考光学影像和反演结果，在同一期地貌面区域内，手动选取多边形范围来确定训练样本的统计范围。在选取过程中要避开人为建筑，并确保在后向散射系数影像和光学影像上都可以明显区分。通过统计各样本的均值和方差，获得各期地貌面后向散射系数值分布曲线（图4）。同一期地貌面的后向散射系数符合正态分布，而地貌面的后向散射系数值越低，则其形成年龄越老（Hetz et al.，2016；苏强等，2020）。

F1—F4对应不同期地貌面；AC为现今河床区域；µ为均值，σ为标准差；S为重叠面积图4 各期地貌面后向散射强度值正态分布概率密度曲线离散程度与区分程度Fig.4 Calculation method of dispersion degree and discrimination degree of normal distribution probability density curve of backscatter intensity values of various geomorphic surfacesF1 to F4 correspond to different stages of landforms; AC represents the current riverbed area; µ denotes the mean, σ is the standard deviation; S stands for the overlapping area.

已有研究将后向散射系数作为粗糙度替代参数（Hetz et al.，2016），但实际上后向散射系数受到多种因素的影响，包括雷达系统参数（入射角、极化方式、天线增益和脉冲频率等）、地表参数（粗糙度、土壤含水量和土壤表面温度等）以及表面植被的覆盖 情况（Ulaby et al.，1979；Fung et al.，1992；Mattia et al.，1997；邵芸等，2002；Escorihuela et al.，2007）。尽管文中选取的研究区域地表干旱裸露，但由于雷达成像时的气候条件不同，仍然会导致影响后向散射系数的环境因子存在较大差异。苏强等（2020）通过经验认为，蒸发量最大的高温夏季时期进行成像反演可以获得最佳效果，但这种方式依赖于主观判断、缺乏定量的验证分析。因此，需要一种更准确的方法来选取以粗糙度为主并且受环境因素影响更小的数据源。

后向散射系数与各个因子存在耦合关系（Mattia et al.，1997；Aubert et al.，2011.），如果无法精确量化各个因素，就无法精确地确定粗糙度。因此，需要尽可能使得其他因子相对于粗糙度对后向散射系数的影响最小化。文中拟通过样本统计分析，初步判断所选影像统计指标与地貌分期结果的关联性，并结合相关气象数据产品协助选出最佳数据源。

样本统计分析主要包含各期地貌分布离散程度和重叠度分析。同一期地貌面对应后向散射系数分布越集中，表明该数据越能反应真实的地貌面年龄，而相邻分布间的高重叠程度会使得重叠范围的归属类错分概率增大，不利于分类（Zhou，2021；Su et al.，2023）。文中离散程度采用分布曲线标准差表征，重叠度采用相邻分布概率密度函数曲线重叠面积表征，曲线重叠面积使用黎曼求和方法计算。

2.4 反演结果监督分类

基于选取的最佳数据，文中采用最大似然分类法对影像反演结果进行分类（杨鑫，2008；Mather and Tso，2016）。最大似然分类是常用的遥感图像监督分类方法，其基本原理为计算属于预定义类别集合的每一个给定像素的归属概率，然后将该像素归为概率最大的类别。该准则计算过程为：

式中P(ωi|X)表示属于第i类的后验概率，X为观测样本，S表示类别总数。

当样本分布为正态分布时，判别函数为：

式中P(ωi)表示属于第i类的先验概率，∑i表 示样本方差，Ei表示样本均值。此时分类准则为：

3 地貌分期结果

3.1 后向散射系数反演结果

在获得10个时相后向散射系数反演结果后，经过平衡对比度增强拉伸处理,可获得较为直观的、反映后向散射系数强度的图像（图5）。根据每一时相反演结果，并参考光学影像和相关研究(王萍，2003；郑荣章，2005；Zhang and Guo，2013)，在确定为同一期地貌面的范围内，手动选取多边形范围作为统计样本以表征地貌面的分期情况（图5a）。具有相同颜色区域统计结果在分布（均值方差）上十分接近，被认为是同一期地貌面（图5b）。

a—样本选取区域（彩色多边形为样本区域范围，不同颜色区域代表不同期地貌面样本）；b—蓝色多边形范围对应的后向散射系数统计结果图5 地貌面样本选取Fig.5 Selection of geomorphic surface samples(a) Sample selection area (Colored polygons represent the sample area range, different colored areas represent samples of different stages of landforms); (b) Statistical results of the backscatter coefficients corresponding to the blue polygon range

3.2 数据质量评定

研究基于MATLAB编程实现统计分析和监督分类过程。在统计分析中，选择的影像成像时间分布在5月、7月、8月和10月，不同时相的后散射系数分布差异明显（图6）。由于气候环境因素导致河流现今河床区域不断发生变化，固定的统计区域范围无法适用于所有时相数据。因此，在统计过程中，不考虑现今河床区域的后向散射系数。按照2.3节所述方法，对各时相分布进行统计分析。结果显示，成像时间为2008年5月22日的影像具有最高的区分度（重叠面积为0.5653）和较低的离散程度（方差和为2.9039）。而成像时间为2008年7月7日的影像具有最低的离散程度（方差和为2.8590）和较高的区分度（重叠面积为0.6144；表2，表3）。这两幅数据都是后续分期可选择的潜在理想数据源。

表2 不同时相各期地貌面后向散射系数分布方差Table 2 Variance of numerical distribution of backscatter coefficients on different levels of landforms at different time periods

表3 不同时相相邻地貌面后向散射系数分布重叠面积Table 3 Overlapping area of numerical distribution of backscatter coefficients on adjacent geomorphic surfaces at different time periods

F1—F4对应不同期地貌面；AC为现今河床区域图6 不同时相各期地貌面后向散射系数分布概率密度曲线Fig.6 PDF (Probability density function) of backscatter coefficient distribution on different stages of geomorphic surface at different time periodsF1 to F4 correspond to different stages of landforms; AC represents the current riverbed area

研究区域常年干燥、降水量稀少，文中仅考虑地表参数的变化对后向散射系数的影响。但是，从统计结果可以观察到不同时相后向散射系数的反演结果之间存在巨大差距。因此有必要进一步分析各成像时间内的具体气象数据，以探究造成分期差异的原因。

基于NASA官网提供的气象数据（https://disc.gsfc.nasa.gov）和历史天气记录网站（https://rp5.ru）提供的信息，文中对各个成像时期的相关气象条件进行了分析。气象信息包含降水量与地表土壤含水量数据（De Jeu and Owe，2012；Huffman et al.，2019），文中分别计算了目标区域内成像时间当天的平均地表土壤含水量和成像时间前5天内每日的累积降水量，分析成像近日是否有降水发生，作为土壤地表含水量结果的交叉验证。研究区域为干旱沙漠地表，蒸发强度高。研究表明，发生常规降水后4天内，地面表层水持续高速率蒸发，常规降水3～5天后，地面表层5 cm深度内的含水量恢复到降水前水平（孙朋，2017）。此外，文中还获取了成像时间前5天内的最高温度数据，因为气温是影响地表水分蒸发强度的主要因素。

对比分析统计数据（表征区分度的重叠面积和表征离散程度的方差之和）与气象数据后发现，土壤水分、方差和、重叠面积和三者之间存在一定相关性，特别是2008年5月22日之后的数据。文中选择2008年7月7日时相为最佳数据源，该时相具有较好的统计指标，同时地表土壤含水量低，降水量低，温度高。相比之下，2007年10月5日数据源指标最差，其重叠度最高（重叠面积为0.8512），离散程度也最高（方差和为3.2782），同时在该时相成像时地表含水量最高（图7）。

为便于对比趋势，对部分数据进行了平移和缩放，其中方差和×10；重叠面积和×10+30，地表含水量+30kg/m2图7 统计指标与大气条件分析Fig.7 Statistical indicators and atmospheric condition analysisFor ease of trend comparison, partial data has been shifted and scaled, where variance is multiplied by 10; overlapping area is multiplied by 10 plus 30, and surface soil moisture is increased by 30 kg/m².

3.3 数据分期结果

根据公式（4）中的分类准则，当d1(x)=d2(x)时，属于这两类中任一类的概率相等。通过代入具体数据，可计算分类决策面（一维情况下为直线x=a；杨 鑫，2008；Mather and Tso，2016）。对 于 地 貌 面F1概率密度曲线左边界和现今河床AC概率密度曲线右边界，分别使用µ-2σ、µ+2σ（µ为均值，σ为标准差）作为边界值（正态曲线（µ-2σ，µ+2σ）内的面积为95.4%）。按照上述分类规则，最终确定后向散射系数值分类区间（表4），并且得到疏勒河洪积扇的分期结果（图8）。

表4 最大似然法获得的分类区间Table 4 Classification interval obtained by maximum likelihood method

F1—F4对应不同期地貌面；AC为现今河床区域a—成像时间为2007年10月5日；b—成像时间为2008年7月7日；c—C波段RADARSAT（据Zhang and Guo，2013修改）图8 不同数据源划分的疏勒河洪积扇地貌面分期结果Fig.8 Classification of the geomorphic surface of the Shule River alluvial fan(a) Image acquisition date: October 5, 2007; (b) Image acquisition date: July 7, 2008; (c) C-band RADARSAT (modified after Zhang and Guo, 2013)F1 to F4 correspond to different stages of landforms; AC represents the current riverbed area

根据2008年7月7日时相获得的分期结果，疏勒河流域共发育4期河流地貌面。在这个时相下，侵蚀冲沟十分发育，地貌面分布不连续。F4为最老一期地貌面，主要分布在洪积扇东北部和西北部扇缘区域。西北部扇中地区存在一期保留较为完整的地貌面。F3地貌面主要分布在扇中地区，少量以条带状分布在扇缘区域。F2、F1地貌面主要分布在河道周缘地区，受到河流改造作用较强，AC则为现今河床所在区域。

4 讨论

4.1 数据质量影响因素

影响后向散射系数的因素主要包括地表参数和雷达系统参数，还可能包括地表植被覆盖引起的干扰。文中选择该区域临近成像时间的LANDSAT5遥感影像（2007年10月16日与2008年7月28日），计算了该区域的归一化植被指数（NDVI）值（郭铌，2003；Robinson et al.，2017），结果显示，无论夏季（图9a）还是秋季（图9b），除现代河床周边出现少量植被覆盖，该区域绝大部分NDVI值都在0.1以下（裸土或岩石地表），表明植被覆盖并非为造成后向散射系数明显差异的主要原因。

在分析气象条件对后向散射系数分期差异产生的影响时，主要分析了温度、降水量、土壤含水量3种因素，并进行了相互验证（图7）。地表含水量是影响地表复介电常数的主要因素（李森，2007），仅分析最优数据源与最差数据源的反演结果，发现土壤地表含水量与统计指标出现了最为接近的趋势。此外，数据源较好的时相都出现在夏季。尽管研究区域在夏季处于雨季，但高温条件导致地表水分蒸发速度加快，而区域年均蒸发量远大于年均降水量，说明在分析地表含水量时需要更多地考虑地表蒸发的影响。例如，2007年8月22日时相数据，尽管成像时间前一天内该区域出现较大幅度降雨，但高温天气使得土壤地表含水量并不高。这与主观经验的结果基本一致，但需要注意的是，并非所有处于夏季时相的数据都适用于地貌面分期。例如，2008年8月25日时相数据，根据统计指标分析，该时相较高的重叠度会增加相邻两期地貌面之间错误划分概率，这表明主观经验分析方法具有一定的局限性。

单独分析统计指标时发现，2007年7月5日时相数据也具有较好的离散程度（方差和为2.9139）和区分度（重叠面积为0.6439），但是气象数据分析表明，该时相在成像前出现较大规模降雨，成像时地表含水量较高。研究认为，该时相气象因素使得后向散射系数反演结果整体偏高（图6），但并未影响反演结果的离散程度与区分度。而后向散射系数反演结果的整体偏差会影响不同区域分期结果的比对，因此在选取数据过程中要综合考虑统计指标与气象指标。

4.2 分期结果对比

综合分析2007年10月5日时相成像结果和气象因素，推测成像时前一天的小幅度降水以及当日较低的温度造成地面湿度较高，不均匀的地表含水量使得地表物质复介电常数偏高，导致后向散射强度反演结果存在较大误差，分期结果不佳。对比认为，SAR数据质量主要受到成像时间区域土壤地表含水量的影响，而与季节相关性不大。因此，在未经过后验定量分析的情况下，成像时间处于高温且降雨稀少天气的影像是最理想数据源。

经过相同的地貌面分期方法处理，得到最差数据源2007年10月5日与最佳数据源2008年7月7日时相的疏勒河洪积扇分期结果（图8a、8b），两次结果差异明显。前者分期结果中F4、F3面积较小，F1、F2地貌面较大，同一期的地貌面分布也更为散乱。东部区域差异最明显，在光学影像和2008年7月7日时相SAR影像上，被认为是同一期的大量地貌面在2007年10月5日时相分期结果中表现为杂乱分布的零散单元。这可能是由于用于统计的样本单个分布的离散程度较高，并且相邻分布之间存在较高的重叠度，导致后续样本的错分概率增大。而北部区域大部分地貌面出现了完全不同的分期结果，但是在扇面现代河床附近少部分区域，2007年10月5日时相得到了更清晰的分期结果。整体来看，2008年7月7日地貌面分期结果明显好于2007年10月5日的分期结果。

将文中地貌面分期结果与Zhang and Guo（2013）得出的结果对比（图8b、8c），其主要区别在于F4与F1地貌面的划分。Zhang and Guo（2013）研究认为在扇缘地区存在一期更古老的地貌面为F4，而F3地貌面区域与文中F4地貌面空间位置基本对应。文中分期结果可以明显地区分出比C波段分期F1更年轻一期的地貌面以及现今河床区域，并且对于较年轻的F2、F3地貌面（空间位置对应Zhang and Guo（2013）划分的F1、F2地貌面）分期的边界特征更明显，这与Ayari et al.（2021）研究中指出的L波段对于高粗糙度地貌面更加敏感相一致。由于研究区域、研究尺度并不完全相同，更加具体细微的差异无法进行对比。文中分期方法的另一优势在于获得数据更方便、分期过程更简单。Zhang and Guo（2013）给出的C波段分期结果是通过RADARSAT-2全极化数据计算SAR多种特征参数得出的，而文中仅采用的是HH单极化数据。总体而言，文中分期方案对于不同阶段年龄地貌面的划分更具优势，对于大范围河流地貌定量分期与对比更具有潜力和进一步的应用价值。

5 结论

（1）研究建立了一套完整地分析和处理L波段SAR数据并进行地貌面分期的系统方案。首先，使用ALOS PARSAR数据反演得到区域内后向散射系数值；然后，通过分析不同时相影像后验统计指标和大气条件来评估数据质量；最后使用统计学分类方法（最大似然分类法）完成河流地貌面分期。并以疏勒河地区洪积扇为例，选定2008年7月7日时相为最佳分期结果，将地貌面分为4期地貌面和现今河床区域。

（2）地貌面分期结果与后验统计分析表现出的明显差异表明了选取数据的重要性。总体而言，影像质量和分期结果与成像时大气条件紧密相关，夏季高蒸发强度的数据通常具有较好的质量。主观经验在选择最佳时相方面具有局限性，目前仍未有一种准确选取最佳时相的先验参考指标的方法。因此，最为严谨的方案仍是下载研究区域内所有时相数据进行评估。

（3）文中建立的地貌面划分的数据选取方案与定量分期方法可广泛用于干旱地区洪积扇地貌面的定量分期。未来的研究将进一步通过编程实现地貌面的自动化定量分期，从而实现大范围地区河流地貌面快速分期。这将为地质构造活动、地貌演化过程、气候环境变化等研究提供重要的地貌分期参考依据。