基于快速傅立叶变换(FFT)的地形剖面线的数学表达研究

胡 洁，杨勤科，郭兰勤

(1.西北农林科技大学资源环境学院，陕西 杨凌 712100;2.西北大学城市与环境学院，西安 710069)

1 研究背景

根据传统地貌学的研究，地形剖面线是沿固定位置和水平走向采集高程(或者其它信息)信息，并以距离为自变量(x轴)、以高程为应变量(y轴)的曲线。剖面分析既是传统地貌学研究的方法［1］，也是数字地形分析的基本方法之一［2－3］。所谓的信号是随时间(空间)变化的物理量［4］。以时间特征量作为自变量来表示信号称为时域信号，以空间特征量作为自变量来表示信号称为空间域信号［5］。可见，地形剖面线就是一个以距离特征量为自变量的一维空间域信号。在时间域信号中，时间频率被定义为单位时间内周期重复的次数，单位为Hz(即周/秒)［6］;类似地在空间域信号中，空间频率被定义为单位距离内周期重复的次数(单位为Hz)，即单位距离内变化或出现次数较多的为高频信息，反之就是低频信息。与信号由不同频率成分构成一样，地形剖面线也包含了不同空间域频率的地形信息。

在地形分析中，地形剖面分析能够直观反映地形表面形态，结合构造地质，沉积地质等研究方法，可以探讨活动断裂特征，盆山系统的耦合关系［7］，区域地表隆升、剥蚀演化及其坡度分布等之间的定量或半定量关系［8－9］。在土壤侵蚀地貌研究中，研究者用剖面线来分析侵蚀地貌的纵向、横向和垂直结构特征［10］。已有研究对剖面线的走向和位置关注不够，也没有进行空间频率的分析，限制了对地貌和侵蚀特征的定量分析。

本文基于较高分辨率DEM，布设具有固定位置的地形剖面;利用快速傅立叶变换(FFT)实现空间域和频率域的转换，并分析地形剖面所包含的地形频率信息，结合对多种较低分辨率DEM所表达的侵蚀沟道即对应的梁峁地形分析，将综合地形剖面分解为具有不同空间频率的成分，探索不同频率成分侵蚀地形的数学表达。本研究对于促进数字地形分析与土壤侵蚀分析的有机结合，对从DEM上辩识更多侵蚀地貌形态信息，将是一个尝试。

2 材料与方法

2.1 基础数据与剖面布设

研究区域:本研究以位于典型黄土丘陵沟壑区的县南沟为研究区。县南沟流域位于陕北黄土高原中部，总面积约44.85km2，属延河流域子流域。该流域土壤侵蚀强烈，侵蚀地形发育，是研究侵蚀地貌形态的理想研究场所。

数据基础:本研究以县南沟1∶10 000比例尺DLG(包括等高线、高程点和河流等专题层)为基础，根据我们前期研究率定的参数［11］，利用 ANUDEM 生成一系列分辨率的DEM，如:2.5，5，10，25，30，40，50，60，70，80，90，100，110，120，130，140，150m;并且在ARCMAP中利用resample中的双线性算法，将这些DEM转为同一栅格尺寸(即 cellsize为2.5)。

剖面布设:沿着较低分辨率DEM的某段等高线在高分辨率DEM上布设剖面，并让剖面尽可能布设在沟壑中部(见图1)，然后将这段等高线转换成1m间距的点文件，并将转换后的点文件在2.5m高分辨率DEM上采集高程值，进而生成研究所需的综合地形高程剖面记作剖面1(见图2)。

图1 2.5m分辨率DEM中剖面布设Fig.1 Location of profile in DEM of 2.5m resolution

图2 地形剖面1Fig.2 The terrain profile 1

2.2 研究方法与研究原理

2.2.1 频段确定与频率分解

本文将地形图上表达的沟壑粗分为3类，暂且称为较大沟、中等沟和较小沟。通过观察不同分辨率DEM表面对侵蚀地形的表达，可以得到不同分辨率DEM所表达的地形信息(表1)。如图3，2.5m高分辨率DEM上有本研究所需的3种沟，25m分辨率DEM则是只包含中等沟和较大沟地形，而150m分辨率DEM中只包含较大沟地形。然后，在2.5，25，150m分辨率DEM上分别生成剖面2、剖面3、剖面4(见图4)。并使DEM相减，将所需的中频成分分离出来，即用25m分辨率DEM减去150m分辨率DEM，使所得结果只表达中等沟信息，在相减的DEM上取剖面5(见图5)。其次，通过对这4个剖面快速傅立叶变换后的频谱分析，可以确定剖面线中中等沟的频率即中频地形成分的范围，即也确定低频的最大值和高频的最小值。最后，利用中频成分的范围和快速傅立叶变换后的频谱分析，将地形剖面所包含的地形频率成分进行分解，可得到低频地形层、中频地形层和高频地形层。

表1 不同分辨率DEM的地形分析Table 1 Terrain analysis of DEMs of different resolutions

图3 不同分辨率DEMFig.3 Contrast of DEMs of different resolutions

2.2.2 傅里叶变换和傅立叶级数

傅里叶变换就是建立时域信号f(t)和与之对应的频域信号 F(jω)的关系［12］。信号的时域(空间域)描述简单直观，只能反应信号的幅值随时间(空间)变换的特性。而频域信号则可用以分析信号的频率构成和各频率成分的幅值大小和相位关系。快速傅里叶变换(FFT)由Cooley和Tukey于1965年提出，它巧妙地利用WN因子的周期性和对称性，构造了离散傅立叶变换(DFT)的快速算法［13］。本文就是采用FFT，将地形剖面线信号变换到频域，通过对该信号进行频谱分析，来研究地形剖面线的地形频率构成。

根据傅立叶级数的理论，任何一个周期为T的周期信号f(t)，如果满足狄里克雷条件，都可以展开成如下的三角级数。对于有限长的非周期信号，也可以用周期信号来实现［14］，因此，对于地形剖面线这样一个有限长非周期信号也可以用傅立叶级数来表示，即地形剖面线可以由不同频率的正弦和或余弦和的形式乘以不同的系数来表示。

若在MATLAB中对某一维信号y(t)进行FFT，其采样频率为fs(此时 fs符合采样定理)［15－16］，采样点数为N(为了运算简单，一般设fs与N值相等)。在该信号的幅值频谱图中，可得到主要的频率成分为f1，f2，f3…fn，以及它们对应的幅值 F1，F2，F3…FN，并且从相位频谱图中得到对应的相位ω1，ω2，ω3…ωn，根据快速傅里叶变换和傅里叶级数原理，就可以推出该信号y(t)的表达式(如式1):

其中，频率ω为0，即该成分为直流分量，其幅度为幅值频谱图中对应的幅值F除以采样点数N。

图4 3种侵蚀沟的剖面Fig.4 The terrain profiles of three erosion ditches

3 结果与分析

3.1 地形剖面的地貌学分析

据调查和制图研究，研究区域侵蚀地貌可划分为3～4个等级的梁和相应的沟［10］。类似于土壤侵蚀研究中的浅沟、切沟和冲沟及其对应的分水。本文将研究区县南沟流域的侵蚀沟简单分为3级即较大沟、中等沟和较小沟，并分别在剖面2、剖面3、剖面4上选取较小沟、中等沟、较大沟的地形剖面(见图4)，并基于2.5m分辨率DEM生成的综合地形剖面1(图2)。这些剖面均具有固定的位置、走向，也包含了研究区的3种侵蚀地形即较大沟、中等沟、较小沟，可作为本研究对地形频率分析的基础。

3.2 频谱分析

在幅值频谱图中可分析该信号频率构成及各频率成分对应的幅值。在地形剖面的幅值频谱图中，幅值与其对应频率地形的绝对高程成正比关系，即幅值越大，对应频率地形的相对高差也越大。相位频谱图则反应了各频率成分对应的相位，而相位对地形剖面复杂的起伏度有一定的影响。这是因为相位影响其对应频率地形的高程极值点出现的位置，进而影响到对应频率地形的起伏特征。

对剖面5(图5)的频谱分析，可得到中等沟频率的范围，而对剖面2和剖面4的频谱分析可分别得到较小沟和较大沟频率的范围，通过对布设不同剖面进行分析及对所得频率进行合适的取舍，可得到中等沟的频率范围:2.25＜fmid＜62.5(周/km)。

图5 布设地形剖面5Fig.5 Locating the terrain profile 5

通过对地形剖面1进行FFT，得到该地形剖面信号的幅值频谱图和相位频谱图。在Matlab中，通过观察幅值频谱图和相位频谱图，将得到地形成分的频率、振幅和相位。结合式(1)和中频剖面的频率范围，不仅可得到该地形剖面线的表达式，还可得到各频率成分的表达式，即将综合地形剖面线分解为低频成分(flow)、中频成分(fmid)和高频成分(fhigh)，其和等于综合剖面。

其中，低频成分的表达式为

中频成分的表达式为

高频成分的表达式为

3.3 地形频率划分

结合综合剖面的频谱图和已知中频地形的频率范围可看出:该地形剖面包含了低、中、高3种地形频率成分，即地形剖面中包含了较大沟、中等沟和较小沟侵蚀地形。在MATLAB中将3种地形频率进行分离，并将低频地形成分与原综合地形剖面、中频地形成分与低频地形成分、中频地形成分与高频地形成分分别进行对比(见图6)。

图6 3种地形频率成分对比Fig.6 Comparison of the three terrain frequency components

如图6，剖面线的低频成分对剖面线整体起伏趋势和地形相对高差影响较大，低频成分相对于中频成分，其高程起伏较大且单位距离内起伏变化小;同理，中频相对于高频成分其高程起伏较大、单位距离内起伏变化小。高频成分主要影响地形剖面起伏的复杂度即粗糙度，其相对起伏小、单位距离内起伏变化大。

4 结论与讨论

本文通过对不同分辨率的DEM进行分析，可以将DEM中所表现的侵蚀沟简单分为较小沟、中等沟、较大沟3种侵蚀地形。在不同分辨率DEM上布设地形剖面，并从信号的角度对地形剖面进行分析，可以得到中等沟频率即中频的范围:2.25＜fmid＜62.5(fx/dx)。利用该频率范围，本文对2.5m高分辨率DEM的某段综合地形剖面进行分析，可发现，该地形剖面中包含较大沟、中等沟、较小沟3种侵蚀地形。通过对该地形剖面的频谱分析及傅里叶变换原理，可得出该地形剖面线中3种侵蚀地形(较大沟、中等沟、较小沟)的数学表达，进而就得到地形剖面线的表达式。总之，利用高分辨率DEM进行高程剖面分析，实现地形剖面的数学表达，可对不同尺度的侵蚀地貌地形特征实现定量或半定量研究，也可推动DEM剖面分析在地学研究特别是在侵蚀区地形地貌研究中的应用。

本研究只是粗略地将复杂的侵蚀地形分成3类，如何根据实际侵蚀地形将剖面线中地形频率进行细分将是下一步研究的重点。另外，基于DEM的地形剖面的频率划分与该DEM图像本身的频率划分之间有什么必然的关系，也需要进一步探讨。

［1］(苏)邦达楚克.地貌学［M］.北京:地质出版社，1954.(B.Γ.ZHUK.Geomorphography［M］.Beijing:Geological Publishing Press，1954.(in Chinese))

［2］汤国安，闾国年.数字高程模型及地学分析的原理与方法［M］.北京:科学出版社，2005.(TANG Guo-an，LU Guo-nian.Digital Elevation Model and the Principles of Earth Science Analysis［M］.Beijing:Science Press，2005.(in Chinese))

［3］LONGLEY P A，GOODCHILD M F，MAGUIRE D J，etal.Geographical Information Systems:Principles，Techniques，Management and Applications，(2nd Edition，A-bridged)［M］.New York:John Wiley ＆ Sons，Inc.，2005.

［4］林 梓.信号与线性系统分析基础［M］.北京:北京邮电大学出版社，2005.(LIN Zi.Signal and Linear System Analysis［M］.Beijing:Beijing University of Posts and Telecommunications Press，2005.(in Chinese))

［5］康宜华，杨叔子.空间域信号的采样方法［J］.华中理工大学学报，1992，(增1):183－188.(KANG Yi-hua，YANG Shu-zi.A Sampling Method for Spatial Signals［J］.Journal of Huazhong University of Science and Technology，1992，(Sup.1):183－188.(in Chinese))

［6］李宗扬.时间频率计量［M］.北京:原子能出版社，2002.(LI Zong-yang.Time and Frequency Measurement［M］.Beijing:Atomic Energy Press，2002.(in Chinese))

［7］张会平，刘少峰.利用DEM进行地形高程剖面分析的新方法［J］.地学前缘，2004，11(3):226－226.(ZHANGHui-ping，LIU Shao-feng.A New Method of Analyzing Terrain Elevation Profile by Use of DEM［J］.Earth Science Frontiers，2004，11(3):226－226.(in Chinese))

［8］SUMMERFIELD M A.Geomorphology and Global Tectonics［M］.New York:John Wiley ＆ Sons，Inc.，2000.

［9］李 勇，DENSMORE A L，周荣军.青藏高原东缘数字高程剖面及其对晚新生代河流下切深度和下切速率的约束［J］.第四纪研究，2006，26(2):236－243.(LI Yong，DENSMORE A L，ZHOU Rong-jun.Profiles of Digital Elevation Models(DEMS)Crossing the Eastern Margin of the Tibetan Plateau and Their Constraints on Dissection Depths and Incision Rates of Late Cenozoic Rivers［J］.Quaternary Sciences，2006，26(2):236－243.(in Chinese))

［10］中国科学院西北水土保持研究所.黄土高原杏子河流域自然资源与水土保持［M］.西安:陕西科学技术出版社，1986.(Institute of Soil and Water Conservation of CAS＆ MWR.The Natural Resources and Soil and Water Conservation of the Xingzihe Watershed［M］.Xi’an:Shannxi Science and Technology Press，1986.(in Chinese))

［11］张彩霞，杨勤科，段建军.高分辨率数字高程模型的构建方法［J］.水利学报，2006，37(8):1009－1014.(ZHANG Cai-xia，YANGQin-ke，DUAN Jian-jun.Method for Establishing High Resolution Digital Elevation Model［J］.Journal of Hydraulic Engineering，2006，37(8):1009－1014.(in Chinese))

［12］高西全，丁玉美，阔永红.数字信号处理:原理，实现及应用［M］:北京:电子工业出版社，2006.(GAO Xiquan，DING Yu-mei，KUO Yong-hong.Digital Signal Processing:Principles，Implementation and Application［M］.Beijing:Electronic Industry Press，2006.(in Chinese))

［13］NAZEERUDDIN M K，DE ANGELIS F.FANTACCI S.Combined Experimental and DFT-TDDFT Computational Study of Photoelectrochemical Cell Ruthenium Sensitizers［J］.Journal of the American Chemical Society，2005，127(48):16835－16847.

［14］冈萨雷斯R C.数字图像处理［M］.北京:电子工业出版社，2003.(GONZALEZ R C.Digital Image Processing［M］.Beijing:Electronic Industry Press，2003.(in Chinese))

［15］刘进明，应怀樵.FFT谱连续细化分析的富里叶变换法［J］.振动工程学报，1995，8(2):162－166.(LIU Jinming，YING Huai-qiao.Zoom FFT Spectrum by Fourier Transform［J］.Journal of Vibration Engineering，1995，8(2):162－166.(in Chinese))

［16］郭 瑜，汤宝平，纪跃波.基于解析信号和带通滤波的频率细化分析［J］.重庆大学学报 (自然科学版)，2001，24(4):l7－22.(GUO Yu，TANG Bao-ping，JI Yue-bo.Zoom FFT Technology Based on Analytic Signal and Band Pass Filter［J］.Journal of Chongqing University(Natural Science Edition)，2001，24(4):l7－22.(in Chinese ))