基于快速傅立葉變換(FFT)的地形剖面線的數(shù)學(xué)表達(dá)研究

胡 潔，楊勤科，郭蘭勤

(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］。所謂的信號(hào)是隨時(shí)間(空間)變化的物理量［4］。以時(shí)間特征量作為自變量來(lái)表示信號(hào)稱為時(shí)域信號(hào)，以空間特征量作為自變量來(lái)表示信號(hào)稱為空間域信號(hào)［5］?？梢?，地形剖面線就是一個(gè)以距離特征量為自變量的一維空間域信號(hào)。在時(shí)間域信號(hào)中，時(shí)間頻率被定義為單位時(shí)間內(nèi)周期重復(fù)的次數(shù)，單位為Hz(即周/秒)［6］;類似地在空間域信號(hào)中，空間頻率被定義為單位距離內(nèi)周期重復(fù)的次數(shù)(單位為Hz)，即單位距離內(nèi)變化或出現(xiàn)次數(shù)較多的為高頻信息，反之就是低頻信息。與信號(hào)由不同頻率成分構(gòu)成一樣，地形剖面線也包含了不同空間域頻率的地形信息。

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

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

2 材料與方法

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

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

數(shù)據(jù)基礎(chǔ):本研究以縣南溝1∶10 000比例尺DLG(包括等高線、高程點(diǎn)和河流等專題層)為基礎(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間距的點(diǎn)文件，并將轉(zhuǎn)換后的點(diǎ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表面對(duì)侵蝕地形的表達(dá)，可以得到不同分辨率DEM所表達(dá)的地形信息(表1)。如圖3，2.5m高分辨率DEM上有本研究所需的3種溝，25m分辨率DEM則是只包含中等溝和較大溝地形，而150m分辨率DEM中只包含較大溝地形。然后，在2.5，25，150m分辨率DEM上分別生成剖面2、剖面3、剖面4(見圖4)。并使DEM相減，將所需的中頻成分分離出來(lái)，即用25m分辨率DEM減去150m分辨率DEM，使所得結(jié)果只表達(dá)中等溝信息，在相減的DEM上取剖面5(見圖5)。其次，通過對(duì)這4個(gè)剖面快速傅立葉變換后的頻譜分析，可以確定剖面線中中等溝的頻率即中頻地形成分的范圍，即也確定低頻的最大值和高頻的最小值。最后，利用中頻成分的范圍和快速傅立葉變換后的頻譜分析，將地形剖面所包含的地形頻率成分進(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 傅里葉變換和傅立葉級(jí)數(shù)

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

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

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

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

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

3 結(jié)果與分析

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

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

3.2 頻譜分析

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

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

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

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

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

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

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

3.3 地形頻率劃分

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

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

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

4 結(jié)論與討論

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

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

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

［2］湯國(guó)安，閭國(guó)年.數(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］林 梓.信號(hào)與線性系統(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］康宜華，楊叔子.空間域信號(hào)的采樣方法［J］.華中理工大學(xué)學(xué)報(bào)，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］李宗揚(yáng).時(shí)間頻率計(jì)量［M］.北京:原子能出版社，2002.(LI Zong-yang.Time and Frequency Measurement［M］.Beijing:Atomic Energy Press，2002.(in Chinese))

［7］張會(huì)平，劉少峰.利用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ù)字高程剖面及其對(duì)晚新生代河流下切深度和下切速率的約束［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］中國(guó)科學(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é)報(bào)，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ù)字信號(hào)處理:原理，實(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］.振動(dòng)工程學(xué)報(bào)，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ì)躍波.基于解析信號(hào)和帶通濾波的頻率細(xì)化分析［J］.重慶大學(xué)學(xué)報(bào) (自然科學(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 ))