胡 潔，楊勤科，郭蘭勤

(1.西北農(nóng)林科技大學(xué)資源環(huán)境學(xué)院，陜西 楊凌 712100;2.西北大學(xué)城市與環(huán)境學(xué)院，西安 710069)

1 研究背景

根據(jù)傳統(tǒng)地貌學(xué)的研究，地形剖面線是沿固定位置和水平走向采集高程(或者其它信息)信息，并以距離為自變量(x軸)、以高程為應(yīng)變量(y軸)的曲線。剖面分析既是傳統(tǒng)地貌學(xué)研究的方法［1］，也是數(shù)字地形分析的基本方法之一［2－3］。所謂的信號是隨時間(空間)變化的物理量［4］。以時間特征量作為自變量來表示信號稱為時域信號，以空間特征量作為自變量來表示信號稱為空間域信號［5］?？梢姡匦纹拭婢€就是一個以距離特征量為自變量的一維空間域信號。在時間域信號中，時間頻率被定義為單位時間內(nèi)周期重復(fù)的次數(shù)，單位為Hz(即周/秒)［6］;類似地在空間域信號中，空間頻率被定義為單位距離內(nèi)周期重復(fù)的次數(shù)(單位為Hz)，即單位距離內(nèi)變化或出現(xiàn)次數(shù)較多的為高頻信息，反之就是低頻信息。與信號由不同頻率成分構(gòu)成一樣，地形剖面線也包含了不同空間域頻率的地形信息。

在地形分析中，地形剖面分析能夠直觀反映地形表面形態(tài)，結(jié)合構(gòu)造地質(zhì)，沉積地質(zhì)等研究方法，可以探討活動斷裂特征，盆山系統(tǒng)的耦合關(guān)系［7］，區(qū)域地表隆升、剝蝕演化及其坡度分布等之間的定量或半定量關(guān)系［8－9］。在土壤侵蝕地貌研究中，研究者用剖面線來分析侵蝕地貌的縱向、橫向和垂直結(jié)構(gòu)特征［10］。已有研究對剖面線的走向和位置關(guān)注不夠，也沒有進(jìn)行空間頻率的分析，限制了對地貌和侵蝕特征的定量分析。

本文基于較高分辨率DEM，布設(shè)具有固定位置的地形剖面;利用快速傅立葉變換(FFT)實現(xiàn)空間域和頻率域的轉(zhuǎn)換，并分析地形剖面所包含的地形頻率信息，結(jié)合對多種較低分辨率DEM所表達(dá)的侵蝕溝道即對應(yīng)的梁峁地形分析，將綜合地形剖面分解為具有不同空間頻率的成分，探索不同頻率成分侵蝕地形的數(shù)學(xué)表達(dá)。本研究對于促進(jìn)數(shù)字地形分析與土壤侵蝕分析的有機(jī)結(jié)合，對從DEM上辯識更多侵蝕地貌形態(tài)信息，將是一個嘗試。

2 材料與方法

2.1 基礎(chǔ)數(shù)據(jù)與剖面布設(shè)

研究區(qū)域:本研究以位于典型黃土丘陵溝壑區(qū)的縣南溝為研究區(qū)?？h南溝流域位于陜北黃土高原中部，總面積約44.85km2，屬延河流域子流域。該流域土壤侵蝕強烈，侵蝕地形發(fā)育，是研究侵蝕地貌形態(tài)的理想研究場所。

數(shù)據(jù)基礎(chǔ):本研究以縣南溝1∶10 000比例尺DLG(包括等高線、高程點和河流等專題層)為基礎(chǔ)，根據(jù)我們前期研究率定的參數(shù)［11］，利用 ANUDEM 生成一系列分辨率的DEM，如:2.5，5，10，25，30，40，50，60，70，80，90，100，110，120，130，140，150m;并且在ARCMAP中利用resample中的雙線性算法，將這些DEM轉(zhuǎn)為同一柵格尺寸(即 cellsize為2.5)。

剖面布設(shè):沿著較低分辨率DEM的某段等高線在高分辨率DEM上布設(shè)剖面，并讓剖面盡可能布設(shè)在溝壑中部(見圖1)，然后將這段等高線轉(zhuǎn)換成1m間距的點文件，并將轉(zhuǎn)換后的點文件在2.5m高分辨率DEM上采集高程值，進(jìn)而生成研究所需的綜合地形高程剖面記作剖面1(見圖2)。

圖1 2.5m分辨率DEM中剖面布設(shè)Fig.1 Location of profile in DEM of 2.5m resolution

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

2.2 研究方法與研究原理

2.2.1 頻段確定與頻率分解

本文將地形圖上表達(dá)的溝壑粗分為3類，暫且稱為較大溝、中等溝和較小溝。通過觀察不同分辨率DEM表面對侵蝕地形的表達(dá)，可以得到不同分辨率DEM所表達(dá)的地形信息(表1)。如圖3，2.5m高分辨率DEM上有本研究所需的3種溝，25m分辨率DEM則是只包含中等溝和較大溝地形，而150m分辨率DEM中只包含較大溝地形。然后，在2.5，25，150m分辨率DEM上分別生成剖面2、剖面3、剖面4(見圖4)。并使DEM相減，將所需的中頻成分分離出來，即用25m分辨率DEM減去150m分辨率DEM，使所得結(jié)果只表達(dá)中等溝信息，在相減的DEM上取剖面5(見圖5)。其次，通過對這4個剖面快速傅立葉變換后的頻譜分析，可以確定剖面線中中等溝的頻率即中頻地形成分的范圍，即也確定低頻的最大值和高頻的最小值。最后，利用中頻成分的范圍和快速傅立葉變換后的頻譜分析，將地形剖面所包含的地形頻率成分進(jìn)行分解，可得到低頻地形層、中頻地形層和高頻地形層。

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

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

2.2.2 傅里葉變換和傅立葉級數(shù)

傅里葉變換就是建立時域信號f(t)和與之對應(yīng)的頻域信號 F(jω)的關(guān)系［12］。信號的時域(空間域)描述簡單直觀，只能反應(yīng)信號的幅值隨時間(空間)變換的特性。而頻域信號則可用以分析信號的頻率構(gòu)成和各頻率成分的幅值大小和相位關(guān)系。快速傅里葉變換(FFT)由Cooley和Tukey于1965年提出，它巧妙地利用WN因子的周期性和對稱性，構(gòu)造了離散傅立葉變換(DFT)的快速算法［13］。本文就是采用FFT，將地形剖面線信號變換到頻域，通過對該信號進(jìn)行頻譜分析，來研究地形剖面線的地形頻率構(gòu)成。

根據(jù)傅立葉級數(shù)的理論，任何一個周期為T的周期信號f(t)，如果滿足狄里克雷條件，都可以展開成如下的三角級數(shù)。對于有限長的非周期信號，也可以用周期信號來實現(xiàn)［14］，因此，對于地形剖面線這樣一個有限長非周期信號也可以用傅立葉級數(shù)來表示，即地形剖面線可以由不同頻率的正弦和或余弦和的形式乘以不同的系數(shù)來表示。

若在MATLAB中對某一維信號y(t)進(jìn)行FFT，其采樣頻率為fs(此時 fs符合采樣定理)［15－16］，采樣點數(shù)為N(為了運算簡單，一般設(shè)fs與N值相等)。在該信號的幅值頻譜圖中，可得到主要的頻率成分為f1，f2，f3…fn，以及它們對應(yīng)的幅值 F1，F(xiàn)2，F(xiàn)3…FN，并且從相位頻譜圖中得到對應(yīng)的相位ω1，ω2，ω3…ωn，根據(jù)快速傅里葉變換和傅里葉級數(shù)原理，就可以推出該信號y(t)的表達(dá)式(如式1):

其中，頻率ω為0，即該成分為直流分量，其幅度為幅值頻譜圖中對應(yīng)的幅值F除以采樣點數(shù)N。

圖4 3種侵蝕溝的剖面Fig.4 The terrain profiles of three erosion ditches

3 結(jié)果與分析

3.1 地形剖面的地貌學(xué)分析

據(jù)調(diào)查和制圖研究，研究區(qū)域侵蝕地貌可劃分為3～4個等級的梁和相應(yīng)的溝［10］。類似于土壤侵蝕研究中的淺溝、切溝和沖溝及其對應(yīng)的分水。本文將研究區(qū)縣南溝流域的侵蝕溝簡單分為3級即較大溝、中等溝和較小溝，并分別在剖面2、剖面3、剖面4上選取較小溝、中等溝、較大溝的地形剖面(見圖4)，并基于2.5m分辨率DEM生成的綜合地形剖面1(圖2)。這些剖面均具有固定的位置、走向，也包含了研究區(qū)的3種侵蝕地形即較大溝、中等溝、較小溝，可作為本研究對地形頻率分析的基礎(chǔ)。

3.2 頻譜分析

在幅值頻譜圖中可分析該信號頻率構(gòu)成及各頻率成分對應(yīng)的幅值。在地形剖面的幅值頻譜圖中，幅值與其對應(yīng)頻率地形的絕對高程成正比關(guān)系，即幅值越大，對應(yīng)頻率地形的相對高差也越大。相位頻譜圖則反應(yīng)了各頻率成分對應(yīng)的相位，而相位對地形剖面復(fù)雜的起伏度有一定的影響。這是因為相位影響其對應(yīng)頻率地形的高程極值點出現(xiàn)的位置，進(jìn)而影響到對應(yīng)頻率地形的起伏特征。

對剖面5(圖5)的頻譜分析，可得到中等溝頻率的范圍，而對剖面2和剖面4的頻譜分析可分別得到較小溝和較大溝頻率的范圍，通過對布設(shè)不同剖面進(jìn)行分析及對所得頻率進(jìn)行合適的取舍，可得到中等溝的頻率范圍:2.25＜fmid＜62.5(周/km)。

圖5 布設(shè)地形剖面5Fig.5 Locating the terrain profile 5

通過對地形剖面1進(jìn)行FFT，得到該地形剖面信號的幅值頻譜圖和相位頻譜圖。在Matlab中，通過觀察幅值頻譜圖和相位頻譜圖，將得到地形成分的頻率、振幅和相位。結(jié)合式(1)和中頻剖面的頻率范圍，不僅可得到該地形剖面線的表達(dá)式，還可得到各頻率成分的表達(dá)式，即將綜合地形剖面線分解為低頻成分(flow)、中頻成分(fmid)和高頻成分(fhigh)，其和等于綜合剖面。

其中，低頻成分的表達(dá)式為

中頻成分的表達(dá)式為

高頻成分的表達(dá)式為

3.3 地形頻率劃分

結(jié)合綜合剖面的頻譜圖和已知中頻地形的頻率范圍可看出:該地形剖面包含了低、中、高3種地形頻率成分，即地形剖面中包含了較大溝、中等溝和較小溝侵蝕地形。在MATLAB中將3種地形頻率進(jìn)行分離，并將低頻地形成分與原綜合地形剖面、中頻地形成分與低頻地形成分、中頻地形成分與高頻地形成分分別進(jìn)行對比(見圖6)。

圖6 3種地形頻率成分對比Fig.6 Comparison of the three terrain frequency components

如圖6，剖面線的低頻成分對剖面線整體起伏趨勢和地形相對高差影響較大，低頻成分相對于中頻成分，其高程起伏較大且單位距離內(nèi)起伏變化小;同理，中頻相對于高頻成分其高程起伏較大、單位距離內(nèi)起伏變化小。高頻成分主要影響地形剖面起伏的復(fù)雜度即粗糙度，其相對起伏小、單位距離內(nèi)起伏變化大。

4 結(jié)論與討論

本文通過對不同分辨率的DEM進(jìn)行分析，可以將DEM中所表現(xiàn)的侵蝕溝簡單分為較小溝、中等溝、較大溝3種侵蝕地形。在不同分辨率DEM上布設(shè)地形剖面，并從信號的角度對地形剖面進(jìn)行分析，可以得到中等溝頻率即中頻的范圍:2.25＜fmid＜62.5(fx/dx)。利用該頻率范圍，本文對2.5m高分辨率DEM的某段綜合地形剖面進(jìn)行分析，可發(fā)現(xiàn)，該地形剖面中包含較大溝、中等溝、較小溝3種侵蝕地形。通過對該地形剖面的頻譜分析及傅里葉變換原理，可得出該地形剖面線中3種侵蝕地形(較大溝、中等溝、較小溝)的數(shù)學(xué)表達(dá)，進(jìn)而就得到地形剖面線的表達(dá)式。總之，利用高分辨率DEM進(jìn)行高程剖面分析，實現(xiàn)地形剖面的數(shù)學(xué)表達(dá)，可對不同尺度的侵蝕地貌地形特征實現(xiàn)定量或半定量研究，也可推動DEM剖面分析在地學(xué)研究特別是在侵蝕區(qū)地形地貌研究中的應(yīng)用。

本研究只是粗略地將復(fù)雜的侵蝕地形分成3類，如何根據(jù)實際侵蝕地形將剖面線中地形頻率進(jìn)行細(xì)分將是下一步研究的重點。另外，基于DEM的地形剖面的頻率劃分與該DEM圖像本身的頻率劃分之間有什么必然的關(guān)系，也需要進(jìn)一步探討。

［1］(蘇)邦達(dá)楚克.地貌學(xué)［M］.北京:地質(zhì)出版社，1954.(B.Γ.ZHUK.Geomorphography［M］.Beijing:Geological Publishing Press，1954.(in Chinese))

［2］湯國安，閭國年.數(shù)字高程模型及地學(xué)分析的原理與方法［M］.北京:科學(xué)出版社，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］林 梓.信號與線性系統(tǒng)分析基礎(chǔ)［M］.北京:北京郵電大學(xué)出版社，2005.(LIN Zi.Signal and Linear System Analysis［M］.Beijing:Beijing University of Posts and Telecommunications Press，2005.(in Chinese))

［5］康宜華，楊叔子.空間域信號的采樣方法［J］.華中理工大學(xué)學(xué)報，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ìn)行地形高程剖面分析的新方法［J］.地學(xué)前緣，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，周榮軍.青藏高原東緣數(shù)字高程剖面及其對晚新生代河流下切深度和下切速率的約束［J］.第四紀(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］中國科學(xué)院西北水土保持研究所.黃土高原杏子河流域自然資源與水土保持［M］.西安:陜西科學(xué)技術(shù)出版社，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］張彩霞，楊勤科，段建軍.高分辨率數(shù)字高程模型的構(gòu)建方法［J］.水利學(xué)報，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］高西全，丁玉美，闊永紅.數(shù)字信號處理:原理，實現(xiàn)及應(yīng)用［M］:北京:電子工業(yè)出版社，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.數(shù)字圖像處理［M］.北京:電子工業(yè)出版社，2003.(GONZALEZ R C.Digital Image Processing［M］.Beijing:Electronic Industry Press，2003.(in Chinese))

［15］劉進(jìn)明，應(yīng)懷樵.FFT譜連續(xù)細(xì)化分析的富里葉變換法［J］.振動工程學(xué)報，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ì)躍波.基于解析信號和帶通濾波的頻率細(xì)化分析［J］.重慶大學(xué)學(xué)報 (自然科學(xué)版)，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 ))