時(shí)間序列小波分析(更新后)

1、時(shí)間序列的小波分析時(shí)間序列（Time Series）是地學(xué)研究中經(jīng)常遇到的問(wèn)題。在時(shí)間序列研究中，時(shí)域和頻域是常用的兩種基本形式。其中，時(shí)域分析具有時(shí)間定位能力，但無(wú)法得到關(guān)于時(shí)間序列變化的更多信息；頻域分析（如Fourier變換）雖具有準(zhǔn)確的頻率定位功能，但僅適合平穩(wěn)時(shí)間序列分析。然而，地學(xué)中許多現(xiàn)象（如河川徑流、地震波、暴雨、洪水等）隨時(shí)間的變化往往受到多種因素的綜合影響，大都屬于非平穩(wěn)序列，它們不但具有趨勢(shì)性、周期性等特征，還存在隨機(jī)性、突變性以及“多時(shí)間尺度”結(jié)構(gòu)，具有多層次演變規(guī)律。對(duì)于這類(lèi)非平穩(wěn)時(shí)間序列的研究，通常需要某一頻段對(duì)應(yīng)的時(shí)間信息，或某一時(shí)段的頻域信息。顯然，時(shí)域分析和頻

2、域分析對(duì)此均無(wú)能為力。20世紀(jì)80年代初，由Morlet提出的一種具有時(shí)-頻多分辨功能的小波分析（Wavelet Analysis）為更好的研究時(shí)間序列問(wèn)題提供了可能，它能清晰的揭示出隱藏在時(shí)間序列中的多種變化周期，充分反映系統(tǒng)在不同時(shí)間尺度中的變化趨勢(shì)，并能對(duì)系統(tǒng)未來(lái)發(fā)展趨勢(shì)進(jìn)行定性估計(jì)。目前，小波分析理論已在信號(hào)處理、圖像壓縮、模式識(shí)別、數(shù)值分析和大氣科學(xué)等眾多的非線性科學(xué)領(lǐng)域內(nèi)得到了廣泛的應(yīng)。在時(shí)間序列研究中，小波分析主要用于時(shí)間序列的消噪和濾波，信息量系數(shù)和分形維數(shù)的計(jì)算，突變點(diǎn)的監(jiān)測(cè)和周期成分的識(shí)別以及多時(shí)間尺度的分析等。一、小波分析基本原理1. 小波函數(shù)小波分析的基本思想是用一簇小

3、波函數(shù)系來(lái)表示或逼近某一信號(hào)或函數(shù)。因此，小波函數(shù)是小波分析的關(guān)鍵，它是指具有震蕩性、能夠迅速衰減到零的一類(lèi)函數(shù)，即小波函數(shù)且滿足： （1）式中，為基小波函數(shù)，它可通過(guò)尺度的伸縮和時(shí)間軸上的平移構(gòu)成一簇函數(shù)系： 其中， （2）式中，為子小波；a為尺度因子，反映小波的周期長(zhǎng)度；b為平移因子，反應(yīng)時(shí)間上的平移。需要說(shuō)明的是，選擇合適的基小波函數(shù)是進(jìn)行小波分析的前提。在實(shí)際應(yīng)用研究中，應(yīng)針對(duì)具體情況選擇所需的基小波函數(shù)；同一信號(hào)或時(shí)間序列，若選擇不同的基小波函數(shù)，所得的結(jié)果往往會(huì)有所差異，有時(shí)甚至差異很大。目前，主要是通過(guò)對(duì)比不同小波分析處理信號(hào)時(shí)所得的結(jié)果與理論結(jié)果的誤差來(lái)判定基小波函數(shù)的好壞，并

4、由此選定該類(lèi)研究所需的基小波函數(shù)。2. 小波變換若是由（2）式給出的子小波，對(duì)于給定的能量有限信號(hào)，其連續(xù)小波變換（Continue Wavelet Transform，簡(jiǎn)寫(xiě)為CWT）為： （3）式中，為小波變換系數(shù)；f(t)為一個(gè)信號(hào)或平方可積函數(shù)；a為伸縮尺度；b平移參數(shù)；為的復(fù)共軛函數(shù)。地學(xué)中觀測(cè)到的時(shí)間序列數(shù)據(jù)大多是離散的，設(shè)函數(shù)，（k=1,2,N; 為取樣間隔），則式（3）的離散小波變換形式為： （4）由式（3）或（4）可知小波分析的基本原理，即通過(guò)增加或減小伸縮尺度a來(lái)得到信號(hào)的低頻或高頻信息，然后分析信號(hào)的概貌或細(xì)節(jié)，實(shí)現(xiàn)對(duì)信號(hào)不同時(shí)間尺度和空間局部特征的分析。實(shí)際研究中，最主要

5、的就是要由小波變換方程得到小波系數(shù)，然后通過(guò)這些系數(shù)來(lái)分析時(shí)間序列的時(shí)頻變化特征。3. 小波方差將小波系數(shù)的平方值在b域上積分，就可得到小波方差，即 (5)小波方差隨尺度a的變化過(guò)程，稱為小波方差圖。由式(5)可知，它能反映信號(hào)波動(dòng)的能量隨尺度a的分布。因此，小波方差圖可用來(lái)確定信號(hào)中不同種尺度擾動(dòng)的相對(duì)強(qiáng)度和存在的主要時(shí)間尺度，即主周期。二、小波分析實(shí)例-時(shí)間序列的多時(shí)間尺度分析(Multi-time scale analysis)例題河川徑流是地理水文學(xué)研究中的一個(gè)重要變量，而多時(shí)間尺度是徑流演化過(guò)程中存在的重要特征。所謂徑流時(shí)間序列的多時(shí)間尺度是指：河川徑流在演化過(guò)程中，并不存在真正意義

6、上的變化周期，而是其變化周期隨著研究尺度的不同而發(fā)生相應(yīng)的變化，這種變化一般表現(xiàn)為小時(shí)間尺度的變化周期往往嵌套在大尺度的變化周期之中。也就是說(shuō)，徑流變化在時(shí)間域中存在多層次的時(shí)間尺度結(jié)構(gòu)和局部變化特征。表1給出了某流域某水文觀測(cè)站1966-2004年的實(shí)測(cè)徑流數(shù)據(jù)。試運(yùn)用小波分析理論，借助Matlab6.5、suffer8.0和相關(guān)軟件（Excel等），完成下述任務(wù)：計(jì)算小波系數(shù)；繪制小波系數(shù)圖（實(shí)部、模和模方）、小波方差圖和主周期變化趨勢(shì)圖，并分別說(shuō)明各圖在分析徑流多時(shí)間尺度變化特征中的作用。表1 某流域某水文觀測(cè)站1966-2004年實(shí)測(cè)徑流數(shù)據(jù)（×108m3）年份徑流量年份徑流

7、量年份徑流量年份徑流量年份徑流量19661.43819742.23519820.77419901.80619981.70919671.15119754.37419830.36719910.44919990.00019680.53619764.21919840.56219920.12020000.00019691.47019772.59019853.04019930.62720012.10419703.47619783.35019860.30419941.65820020.00919714.06819792.54019870.72819951.02520033.17719722.14719800

8、.80719880.49219960.95520040.92119733.93119810.57319890.00719971.341分析1. 選擇合適的基小波函數(shù)是前提在運(yùn)用小波分析理論解決實(shí)際問(wèn)題時(shí)，選擇合適的基小波函數(shù)是前提。只有選擇了適合具體問(wèn)題的基小波函數(shù)，才能得到較為理想的結(jié)果。目前，可選用的小波函數(shù)很多，如Mexican hat小波、Haar小波、Morlet小波和Meyer小波等。在本例中，我們選用Morlet連續(xù)復(fù)小波變換來(lái)分析徑流時(shí)間序列的多時(shí)間尺度特征。原因如下：1.1 徑流演變過(guò)程中包含“多時(shí)間尺度”變化特征且這種變化是連續(xù)的，所以應(yīng)采用連續(xù)小波變換來(lái)進(jìn)行此項(xiàng)分析。1.

9、2實(shí)小波變換只能給出時(shí)間序列變化的振幅和正負(fù)，而復(fù)小波變換可同時(shí)給出時(shí)間序列變化的位相和振幅兩方面的信息，有利于對(duì)問(wèn)題的進(jìn)一步分析。1.3 復(fù)小波函數(shù)的實(shí)部和虛部位相差為/2，能夠消除用實(shí)小波變換系數(shù)作為判據(jù)而產(chǎn)生的虛假振蕩，使分析結(jié)果更為準(zhǔn)確。2. 繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢(shì)圖是關(guān)鍵當(dāng)選擇好合適的基小波函數(shù)后，下一步的關(guān)鍵就是如何通過(guò)小波變換獲得小波系數(shù)，然后利用相關(guān)軟件繪制小波系數(shù)圖、小波方差圖和主周期變化趨勢(shì)圖，進(jìn)而根據(jù)上述三種圖形的變化識(shí)別徑流時(shí)間序列中存在的多時(shí)間尺度。具體步驟1. 數(shù)據(jù)格式的轉(zhuǎn)化2. 邊界效應(yīng)的消除或減小3. 計(jì)算小波系數(shù)4. 計(jì)算復(fù)小波系數(shù)的實(shí)部5

10、. 繪制小波系數(shù)實(shí)部等值線圖6. 繪制小波系數(shù)模和模方等值線圖7. 繪制小波方差圖8. 繪制主周期趨勢(shì)圖下面，我們以上題為例，結(jié)合軟件Matlab 6.5、Suffer 8.0和Excel，詳細(xì)說(shuō)明小波系數(shù)的計(jì)算和各圖形的繪制過(guò)程，并分別說(shuō)明各圖在分析徑流多時(shí)間尺度變化特征中的作用。1. 數(shù)據(jù)格式的轉(zhuǎn)化和保存將存放在Excel表格里的徑流數(shù)據(jù)（以時(shí)間為序排為一列）轉(zhuǎn)化為Matlab 6.5識(shí)別的數(shù)據(jù)格式（.mat）并存盤(pán)。具體操作為：在Matlab 6.5 界面下，單擊“File-Import Data”，出現(xiàn)文件選擇對(duì)話框“Import”后，找到需要轉(zhuǎn)化的數(shù)據(jù)文件（本例的文件名為runoff

11、.xls），單擊“打開(kāi)”。等數(shù)據(jù)轉(zhuǎn)化完成后，單擊“Finish”，出現(xiàn)圖1顯示界面；然后雙擊圖1中的Runoff，彈出“Array Editor: runoff”對(duì)話框，選擇File文件夾下的“Save Workspace As”單擊，出現(xiàn)圖2所示的“Save to MAT-File:”窗口，選擇存放路徑并填寫(xiě)文件名（runoff.mat），單擊“保存”并關(guān)閉“Save to MAT-File”窗口。圖1 數(shù)據(jù)格式的轉(zhuǎn)化圖2數(shù)據(jù)的保存2. 邊界效應(yīng)的消除或減小因?yàn)楸纠械膶?shí)測(cè)徑流數(shù)據(jù)為有限時(shí)間數(shù)據(jù)序列，在時(shí)間序列的兩端可能會(huì)產(chǎn)生“邊界效用”。為消除或減小序列開(kāi)始點(diǎn)和結(jié)束點(diǎn)附近的邊界效應(yīng)，須對(duì)其

12、兩端數(shù)據(jù)進(jìn)行延伸。在進(jìn)行完小波變換后，去掉兩端延伸數(shù)據(jù)的小變換系數(shù)，保留原數(shù)據(jù)序列時(shí)段內(nèi)的小波系數(shù)。本例中，我們利用Matlab 6.5小波工具箱中的信號(hào)延伸（Signal Extension）功能，對(duì)徑流數(shù)據(jù)兩端進(jìn)行對(duì)稱性延伸。具體方法為：在Matlab 6.5界面的“Command Window”中輸入小波工具箱調(diào)用命令“Wavemenu”，按Enter鍵彈“Wavelet Toolbox Main Menu”（小波工具箱主菜單）界面（圖）；然后單擊“Signal Extension”，打開(kāi)Signal Extension / Truncation窗口，單擊“File”菜單下的“Load

13、Signal”，選擇runoff.mat文件單擊“打開(kāi)”，出現(xiàn)圖4信號(hào)延伸界面。Matlab 6.5的Extension Mode菜單下包含了6種基本的延伸方式（Symmetric、Periodic、Zero Padding、Continuous、Smooth and For SWT）和Direction to extend菜單下的3種延伸模式（Both、Left and Right），在這里我們選擇對(duì)稱性兩端延伸進(jìn)行計(jì)算。數(shù)據(jù)延伸的具體操作過(guò)程是：在Extension Mode下選擇“ Symmetric”，Dircetion to extend下選擇“Both”，單擊“Extend”按鈕進(jìn)

14、行對(duì)稱性兩端延伸計(jì)算，然后單擊“File”菜單下的“Save Tranformed Signal”，將延伸后的數(shù)據(jù)結(jié)果存為erunoff.mat文件。從erunoff文件可知，系統(tǒng)自動(dòng)將原時(shí)間序列數(shù)據(jù)向前對(duì)稱延伸12個(gè)單位，向后延伸13個(gè)單位。圖3 小波工具箱主菜單圖4 徑流時(shí)間序列的延伸3. 計(jì)算小波系數(shù)圖5 小波變換菜單界面選擇Matlab 6.5小波工具箱中的Morlet復(fù)小波函數(shù)對(duì)延伸后的徑流數(shù)據(jù)序列（erunoff.mat）進(jìn)行小波變換，計(jì)算小波系數(shù)并存盤(pán)。小波工具箱主菜單界面見(jiàn)圖3，單擊“Wavelet 1-D”下的子菜單“Complex Continuous Wavelet 1-

15、D”，打開(kāi)一維復(fù)連續(xù)小波界面，單擊“File”菜單下的“Load Signal”按鈕，載入徑流時(shí)間序列erunoff.mat（圖5）。圖5的左側(cè)為信號(hào)顯示區(qū)域，右側(cè)區(qū)域給出了信號(hào)序列和復(fù)小波變換的有關(guān)信息和參數(shù)，主要包括數(shù)據(jù)長(zhǎng)度（Data Size）、小波函數(shù)類(lèi)型（Wavelet：cgau、shan、fbsp和cmor）、取樣周期（Sampling Period）、周期設(shè)置（Scale Setting）和運(yùn)行按鈕（Analyze），以及顯示區(qū)域的相關(guān)顯示設(shè)置按鈕。本例中，我們選擇cmor (1-1.5)、取樣周期為1、最大尺度為32，單擊“Analyze”運(yùn)行按鈕，計(jì)算小波系數(shù)。然后單擊“Fi

16、le”菜單下的“Save Coefficients”，保存小波系數(shù)為cerunoff.mat文件。注意：上面涉及到的數(shù)據(jù)保存，其格式均為.mat。4. 計(jì)算Morlet復(fù)小波系數(shù)的實(shí)部將復(fù)小波系數(shù)轉(zhuǎn)存到Excel表格，去掉兩端延伸數(shù)據(jù)的小波系數(shù)，并計(jì)算小波系數(shù)實(shí)部。在Matlab 6.5界面下的Workspace中將cerunoff.mat文件導(dǎo)入，然后雙擊打開(kāi)，全部復(fù)制到Excel后去掉延伸數(shù)據(jù)的小波變換系數(shù)（本例中去掉前12列和后13列），或只復(fù)制原時(shí)間序列的小波變換系數(shù)到Excel，最后使用Excel中的IMREAL函數(shù)計(jì)算原時(shí)間序列的小波系數(shù)實(shí)部（圖6）。圖6 復(fù)小波系數(shù)及實(shí)部計(jì)算示意

17、圖Excel 中IMREAL函數(shù)的調(diào)用方法為：?jiǎn)螕簟安迦搿辈藛蜗碌摹昂瘮?shù)（F）”按鈕，彈出圖7所示的“插入函數(shù)”窗口，在“搜索函數(shù)（S）：”框中輸入：“IMREAL”后單擊“轉(zhuǎn)到”，再單擊“確定”，出現(xiàn)函數(shù)參數(shù)窗口（圖8）。在“Inumber”一欄的空白處輸入所要計(jì)算的數(shù)據(jù)（圖6），單擊“確定”即可得到小波系數(shù)實(shí)部值。圖7 IMEAL 函數(shù)調(diào)用圖8 IMEAL 函數(shù)參數(shù)需要說(shuō)明的是，從cerunoff.mat文件中轉(zhuǎn)到Excel里的復(fù)小波系數(shù)，在其實(shí)部和虛部中間包含許多“空格”，在計(jì)算之前需要先將其去掉。5. 借助Suffer 8.0，繪制小波系數(shù)實(shí)部等值線圖5.1 小波系數(shù)實(shí)部等值線圖的繪制

18、首先，將小波系數(shù)實(shí)部數(shù)據(jù)按照?qǐng)D9格式排列，其中列A為時(shí)間，列B為尺度，列C為不同時(shí)間和尺度下所對(duì)應(yīng)的小波系數(shù)實(shí)部值。圖9 小波系數(shù)實(shí)部數(shù)據(jù)格式其次，將圖9數(shù)據(jù)轉(zhuǎn)化成Suffer 8.0識(shí)別的數(shù)據(jù)格式。具體操作為：在Suffer 8.0界面下，單擊“網(wǎng)格”菜單下的“數(shù)據(jù)”按鈕，在“打開(kāi)”窗口選擇要打開(kāi)的文件（小波系數(shù)實(shí)部.xls），單擊“打開(kāi)”后彈出“網(wǎng)格化數(shù)據(jù)”對(duì)話框（圖10）。它給出了多種不同的網(wǎng)格化方法、文件輸出路徑及網(wǎng)格線索幾何學(xué)等信息。這里我們選擇“克里格“網(wǎng)格方法”，單擊“確定”，完成數(shù)據(jù)格式的轉(zhuǎn)化。圖10 小波系數(shù)實(shí)部數(shù)據(jù)格式轉(zhuǎn)化圖11 Suffer8.0中的小系數(shù)實(shí)部等值線圖最后

19、，繪制小波系數(shù)實(shí)部等值線圖。在Suffer 8.0界面下，單擊“地圖”菜單下的“等值線圖-新建等值線圖”按鈕，彈出“打開(kāi)網(wǎng)格”窗口后，選擇“小波系數(shù)實(shí)部.grd”文件，單擊“打開(kāi)”，完成等值線圖的繪制并存盤(pán)（圖11）。5.2 小波系數(shù)實(shí)部等值線圖在多時(shí)間尺度分析中的作用圖12 小系數(shù)實(shí)部等值線圖小波系數(shù)實(shí)部等值線圖能反映徑流序列不同時(shí)間尺度的周期變化及其在時(shí)間域中的分布，進(jìn)而能判斷在不同時(shí)間尺度上，徑流的未來(lái)變化趨勢(shì)。為能比較清楚的說(shuō)明小波系數(shù)實(shí)部等值線圖在徑流多時(shí)間尺度分析中的作用，我們利用Suffer 8.0對(duì)其進(jìn)一步處理和修飾，得到圖12顯示的小波系數(shù)實(shí)部等值線圖。其中，橫坐標(biāo)為時(shí)間（年

20、份），縱坐標(biāo)為時(shí)間尺度，圖中的等值曲線為小波系數(shù)實(shí)部值。當(dāng)小波系數(shù)實(shí)部值為正時(shí)，代表徑流豐水期，在圖中我們用實(shí)線繪出，“H”表示正值中心；為負(fù)時(shí)，表示徑流枯水期，用虛線繪出，“L”表示負(fù)值中心。由圖12可以清楚的看到徑流演化過(guò)程中存在的多時(shí)間尺度特征。總的來(lái)說(shuō)，在流域徑流演變過(guò)程中存在著1832年，817年以及37年的3類(lèi)尺度的周期變化規(guī)律。其中，在1832年尺度上出現(xiàn)了枯-豐交替的準(zhǔn)兩次震蕩；在817年時(shí)間尺度上存在準(zhǔn)5次震蕩。同時(shí)，還可以看出以上兩個(gè)尺度的周期變化在整個(gè)分析時(shí)段表現(xiàn)的非常穩(wěn)定，具有全域性；而310年尺度的周期變化，在1980s以后表現(xiàn)的較為穩(wěn)定。6. 繪制小波系數(shù)模和模方等

21、值線圖6.1 小波系數(shù)模和模方等值線圖的繪制參考4、5兩步，繪制小波系數(shù)模和模方等值線圖（圖13、14）。說(shuō)明：在Excel中，復(fù)數(shù)模的計(jì)算函數(shù)為“IMABS”。圖13 小波系數(shù)模等值線圖圖14 小波系數(shù)模方等值線圖6.2 小波系數(shù)模等值線圖在多時(shí)間尺度分析中的作用Morlet小波系數(shù)的模值是不同時(shí)間尺度變化周期所對(duì)應(yīng)的能量密度在時(shí)間域中分布的反映，系數(shù)模值愈大，表明其所對(duì)應(yīng)時(shí)段或尺度的周期性就愈強(qiáng)。從圖13可以看出，在流域徑流演化過(guò)程中，1832年時(shí)間尺度模值最大，說(shuō)明該時(shí)間尺度周期變化最明顯，1822年時(shí)間尺度的周期變化次之，其他時(shí)間尺度的周期性變化較?。?6.2 小波系數(shù)模方等值線圖在多

22、時(shí)間尺度分析中的作用小波系數(shù)的模方相當(dāng)于小波能量譜，它可以分析出不同周期的震蕩能量。由圖14知，2532年時(shí)間尺度的能量最強(qiáng)、周期最顯著，但它的周期變化具有局部性（1980s前）；1015年時(shí)間尺度能量雖然較弱，但周期分布比較明顯，幾乎占據(jù)整個(gè)研究時(shí)域(19742004年)。7. 繪制小波方差圖7.1小波方差圖的繪制圖15 小波方差圖將不同時(shí)間尺度下的小波系數(shù)代入式(5)可得徑流變化的小波方差，以小波方差為縱坐標(biāo)，時(shí)間尺度a為橫坐標(biāo)，可繪制小波方差圖（圖15）。7.2小波方差圖在多時(shí)間尺度分析中的作用小波方差圖能反映徑流時(shí)間序列的波動(dòng)能量隨尺度a的分布情況。可用來(lái)確定徑流演化過(guò)程中存在的主周期

23、。流域徑流的小波方差圖中（圖15）存在4個(gè)較為明顯的峰值，它們依次對(duì)應(yīng)著28年、14年、8年和4年的時(shí)間尺度。其中，最大峰值對(duì)應(yīng)著28年的時(shí)間尺度，說(shuō)明28年左右的周期震蕩最強(qiáng)，為流域年徑流變化的第一主周期；14年時(shí)間尺度對(duì)應(yīng)著第二峰值，為徑流變化的第二主周期，第三、第三峰值分別對(duì)應(yīng)著8年和4年的時(shí)間尺度，它們依次為流域徑流的第三和第四主周期。這說(shuō)明上述4個(gè)周期的波動(dòng)控制著流域徑流在整個(gè)時(shí)間域內(nèi)的變化特征。氣象家園提問(wèn)截圖，小波方差圖：8. 主周期趨勢(shì)圖的繪制及其在多時(shí)間尺度分析中的作用根據(jù)小波方差檢驗(yàn)的結(jié)果，我們繪制出了控制流域徑流演變的第一和第二主周期小波系數(shù)圖（圖16）。從主周期趨勢(shì)圖中

24、我們可以分析出在不同的時(shí)間尺度下，流域徑流存在的平均周期及豐-枯變化特征。圖16a顯示，在14年特征時(shí)間尺度上，流域徑流變化的平均周期為9.5年左右，大約經(jīng)歷了4個(gè)豐-枯轉(zhuǎn)換期；而在28年特征時(shí)間尺度上（圖16b），流域的平均變化周期為20年左右，大約2個(gè)周期的豐-枯變化。圖16 大沽夾河流域年徑流變化的13年和28年特征時(shí)間尺度小波實(shí)部過(guò)程線參考文獻(xiàn)王文圣，丁晶，李耀清. 2005. 水文小波分析M. 北京：化學(xué)工業(yè)出版社曹素華等. 1998. 實(shí)用醫(yī)學(xué)多因素統(tǒng)計(jì)方法M. 上海：上海醫(yī)科大學(xué)出版社方開(kāi)泰. 1989. 實(shí)用多元統(tǒng)計(jì)分析M. 上海：華東師范大學(xué)出版社何清波，蘇炳華，錢(qián)亢. 2002. 醫(yī)學(xué)統(tǒng)計(jì)學(xué)及其軟件包M. 上海：上海科學(xué)技術(shù)文獻(xiàn)出版社胡秉民. 1987. 微電腦在農(nóng)業(yè)科學(xué)中的應(yīng)用M. 北京：科學(xué)出版社孫尚拱. 1990. 實(shí)用多元變量統(tǒng)計(jì)方法與計(jì)算程序M. 北京：北京醫(yī)科大學(xué)、中國(guó)協(xié)和醫(yī)科大學(xué)聯(lián)合出版社唐守正. 1986. 多元統(tǒng)計(jì)分析方法M。北京：中國(guó)林業(yè)出版社王學(xué)仁. 1982. 地址數(shù)據(jù)的多變量統(tǒng)計(jì)分析. 北京：科學(xué)出版社徐振邦，金淳浩，婁元仁. 1986. 距離系數(shù)和距離系數(shù)尺度在聚類(lèi)分析中的應(yīng)用M. 趙旭東等主編，中國(guó)數(shù)學(xué)地質(zhì)（1）. 北京：地質(zhì)出版社於崇文. 1978. 數(shù)學(xué)地質(zhì)的方法與應(yīng)用M. 北京：冶金工業(yè)出版社