基于模糊C均值聚類-支持向量機的海底沉積物分類識別

尤加春, 毛慧慧 , 段文豪, 李紅星

(1.中國地質大學 地球物理與信息技術學院, 北京 100083; 2.中國科學院地質與地球物理研究所, 北京100029; 3.中國科學院大學, 北京100049; 4.東華理工大學 核工程與地球物理學院, 江西 撫州344000)

海底底質類型是一種重要的海洋環境參數, 底質類型的分布對海洋科學研究、海洋工程以及國防建設等具有重要的科學與實際意義。對于海底沉積物的探測主要有直接采樣法和間接采樣法, 由于直接采樣法的成本高, 采樣效率低(只能零星采樣)等缺點, 間接采樣方法是目前海底沉積物類型探測的主要方法。海底沉積物的間接采樣分類方法主要是利用聲學、光學、地震學、生物化學等方法接受不同海底沉積物的響應, 并根據這些響應的特征進行分類。

國內外對于海底底質的分類進行了大量的研究。國外早在1947~1948年瑞典科學家Arrhenius就開展了深海勘探計劃。經過數十年的發展[1-4], 國外已經形成了比較成熟的勘探技術, 海底底質分類的儀器和軟件也得到了極大的發展, 例如挪威 Simrad公司的 Triton軟件等。而國內對海底沉積物分類研究起步較晚。孟金生[5]、王正垠[6]利用不同的探測手段[7]研究了海底底質的響應差異并進行分類。對于實驗室基礎理論性的研究, 目前主要是以物理模擬為主, 例如卜英勇等[8]、鄧躍紅等[9]在水槽實驗平臺上建立沉積物模型, 利用超聲回波探測不同的沉積物,并基于回波信號對沉積物進行分類研究。

實驗室中的物理模擬主要是人為地按照一定比例配制海底沉積物樣本, 然后再利用裝置激發接受回波信號。物理模擬方法較為費時費力且實驗可重復性差, 而利用計算機模擬聲波探測海底沉積物并進行分類的研究鮮見于刊。本文嘗試采用計算機數值正演技術模擬實際的地震勘探數據采集過程, 數值模擬方法具有快速、高效、高重復性、易于操作、經濟等特點。基于不同的海底沉積物在聲學上將產生不同強度的回波信號, 本文分別采用一種無監督分類方式: 模糊C均值聚類(Fuzzy C Means, FCM)和一種監督分類方式: 支持向量機(Support VectorMachine, SVM)對提取的回波的特征向量進行分類識別, 并對上述兩種方法作了融合, 提出了一種新型、實用、快速高效的分類方法。此外, 數值模擬方法更有利于研究不同的分類識別方法、屬性提取技術對各種海底沉積物的識別效果并為實際應用提供理論依據。

1 雙相-隨機介質模型

實際的海底沉積介質為典型的雙相介質, 本文認為海底沉積物由巖石骨架(例如礫石、黏土等)和流體填充的孔隙構成, 利用基于雙相-隨機介質的彈性波動方程模擬實際地震勘探數據采集過程。基于Biot理論的雙相介質彈性波動方程[10-13]為:

式中:R為流相彈性系數;Q為固流相耦合彈性系數;e為固相體應變;ε為流相體應變;u為固相位移分量;U為流相位移分量;ρ11為單位體積中固體相對流體運動時固體部分總的等效質量;ρ22為單位體積中流體相對固體運動時流體部分總的等效質量;ρ12為單位體積中流體和固體之間的質量耦合系數;A, N相當于單相各向同性彈性理論中的拉梅常數;b為耗散系數。

為能更加準確地描述實際海底復雜的沉積環境,在雙相介質中引入隨機擾動構成雙相-隨機介質, 這種隨機介質[14-17]中的隨機擾動可以理解為雙相介質中巖石骨架彈性參數的隨機變化。為了簡化計算, 本文只考慮在排空情況下孔隙介質中巖石骨架彈性系數引起的非均一性, 并假設彈性系數各量的相對擾動是相同的, 從而可以只用一個相對擾動量來描述隨機介質在小尺度上的非均勻性[18]。考慮空間隨機介質擾動的孔隙介質彈性系數可以表示為N=N0(1+δ),A=A0(1+δ)。其中,A0,N0是背景介質的彈性參數,δ是空間隨機介質通過一定的自相關函數產生的擾動。本文在研究隨機介質的自相關函數時選擇指數型函數來產生隨機擾動:

其中,a、b分別為介質在x方向和z方向上的自相關長度。

本文在已知泥巖、砂巖、礫巖的縱波、橫波和密度的情況下, 將這些參數轉換為泥質礫巖、泥、泥質砂巖的雙相介質彈性參數, 相關等效介質理論可參考文獻[19-21]。雙相介質彈性波動方程的高階交錯網絡有限差分離散形式參考文獻[10-13, 22]。

2 分類識別方法

2.1 支持向量機

傳統的統計學研究方法都是建立在大數定理這一基礎上的漸進理論, 要求學習樣本數目足夠多。然而在實際應用中, 由于各個方面的原因, 這一前提往往難以得到保證。因此在小樣本情況下, 建立在傳統統計學基礎上的學習方法也就很難取得理想的學習效果和泛化性能。

基于Vapnik的統計學習理論而提出的支持向量機[23-25]大大減小了算法設計的隨意性, 并很好地解決了在上述問題中提到的如何在有限樣本情況下實現機器學習的強泛化能力。利用支持向量機分類的基本思想是通過核函數變換的方法, 將低維空間非線性分類問題轉換為高維空間線性可分的問題, 然后在新空間中求解最優分類面。

在實際的海底沉積物底質調查中, 勘測區域的沉積物樣本是較少的, 研究人員需要對大量的反射波信號進行分類研究, 進而繪制勘測區域的海底沉積物底質分布情況。本文正是在這種實際需求下, 考慮有限樣本情況, 利用支持向量機對海底沉積物進行分類識別研究。

2.1.1 支持向量機基本原理

支持向量機最初用來對線性可分數據進行二值分類處理的, 對線性可分問題, 給定訓練樣本{(x1,y1) ,(x2,y2),… ,(xN,yN)}, 期 望 輸 出y∈ {-1,1},分別代表兩類的類別標識。用于分類的最優超平面方程為:ωTx+b=0, 其中ω為權重向量,x為輸入向量,b為偏置。支持向量機的目的是尋找一個分離邊緣最大的超平面(即最優超平面)使兩類數據最大可能地分離, 離最優超平面最近的特殊樣本稱為支持向量。則對于一確定的最優超平面, 所有的樣本都滿足:

對于非線性可分模式的分類問題, 會有一些樣本不滿足(2)式的約束條件, 而出現分類誤差。因此需要適當放寬該約束條件, 將其變為:

式中引入了松弛變量ξi,i= 1 ,2,L ,N, 它用于量度一個數據點對線性可分理想條件的偏離程度。當0≤ξ≤1時, 數據點落入分離區域的內部, 且在分類超平面的正確一側; 當ξ＞1時, 數據點進入分類超平面的錯誤一側; 當ξ=0時, 退化為線性可分問題。

對于非線性問題, 尋找ω和b的最優值, 使其在(3)式的約束下, 最小化關于ω和ξi的目標函數, 即

其中,Ck為懲罰系數。

利用 Lagrange乘數法求解上述最優化問題, 上述方程可變如下對偶問題:

其中,e為單位矩陣,Q為半正定矩陣,Qij=yiyjK(xi,xj),為核函數, 本文采用的 RBF核函數為, 其中2σ為給定參數。

在實際應用中,Ck和σ的選取對于支持向量機分類的效果影響極其重要, 而目前對于上述參數的選取并無統一的準則, 往往是采用大量試驗的方法來獲得較優的參數值, 但這種方法比較繁瑣、費時,而且獲得的參數也不一定能使分類效果最優。其實,選擇合適的Ck和σ值使支持向量機的分類效果達到最優的問題屬于優化問題。本文提出采用差分進化(Differential Evolution, DE)算法[26]實現對支持向量機參數的自動最優化搜索。

2.1.2 差分進化算法

差分進化算法[28-29]是一種隨機的并行直接搜索算法, 整個算法包含3個過程: 變異、交叉、選擇, 類似于遺傳算法的變異、交叉和選擇操作。3種操作描述如下:

(1) 變異

設有N個個體xi(t),i=1, 2, …,N, 對于第i個體xi(t), 根據下面公式生成下一代變異個體:

其中,xr3(t),xr2(t),xr1(t)為從進化群體中隨機選取的互不相同的 3個個體, 其中i和r1、r2、r3之間必須是不同的。F為縮放比例因子, 用于控制差向量的影響大小。

(2) 交叉

為了增加群體的多樣性, 交叉操作被引入差分進化算法。將個體xi(t)和變異個體進行二項分布雜交 生 成 雜 交 個 體:

具體操作如下:

其中,R∈ [ 0,1], 為雜交參數,P∈ [0,1]之間的隨機數,D為解空間維數。

(3) 選擇

在基本差分進化算法中, 選擇操作采取貪婪策略, 即只有當產生的子代個體優于父代個體時(對應目標函數值f(xi(t+ 1 ))≤f(xi(t)))才被保留, 否則父代個體被保留至下一代。

2.1.3 基于差分進化算法的支持向量機參數優化

本文采用交叉驗證(Cross Validation, CV)的思想利用差分進化算法求取支持向量機的最優化參數。CV是用來驗證分類器性能的一種統計分析方法, 基本思想是在某種意義下將原始數據進行分組, 一部分用于訓練集, 另一部分用于驗證集。其方法是首先用訓練集對分類器進行訓練, 再利用驗證集來測試訓練得到的模型, 以得到的分類準確率作為評價分類性能指標。本文選擇 K-折交叉驗證(K-fold Cross Validation, K-fold CV)來評估分類模型的泛化性能。

為驗證模型的測試效果, 本文以準確率作為差分進化算法的目標函數。準確率定義為: 準確率＝測試樣本中分類正確的樣本數/測試樣本總數。利用差分進化算法對支持向量機中Ck和σ的尋優計算可轉化為一個二維最優化問題。具體操作步驟為:

步驟1: 初始化參數, 給定差分進化算法的種群規模,F值, 最大迭代次數, 精度要求,R、Ck和σ的初始值。

步驟 2: 根據種群參數利用支持向量機對訓練樣本集進行訓練, 并用訓練好的模型對測試集數據進行預測, 計算每個個體的目標函數值(即準確率)。

步驟3: 根據差分進化算法中的變異、交叉由父代種群產生子代種群, 利用子代種群參數對支持向量機進行訓練和測試, 計算子代種群每個個體的目標函數值; 執行差分進化算法中的選擇操作。

步驟 4: 判斷是否滿足計算精度要求或是否達到最大迭代次數, 不滿足時, 返回步驟 3; 否則, 執行步驟5。

步驟5: 迭代結束, 輸出支持向量機模型最優的Ck和σ的值。

2.2 模糊C均值聚類

模糊C均值聚類[27-28]是Bezkek于1981年提出的, 它是目前廣泛采用的一種聚類算法, 其主要思想是將經典劃分的定義模糊化, 用隸屬度來確定屬于某個聚類程度的一種聚類方法。模糊 C均值聚類是模糊聚類算法中非常有效的一種, 即使對于很難明顯分類的變量, 模糊 C均值聚類也能得到較為滿意的效果。

考慮一個樣本集合X={x1,x2,… ,xn}, 將其分為c個模糊組, 并求每組的聚類中心cj(j= 1,2,… ,C0),使目標函數達到最小。目標函數定義如下:

其中:uij∈[0,1]間;ci為模糊組i的聚類中心,C0是期望聚類的數目,為第i個聚類中心與第j個數據點間的歐幾里德距離; 且m∈ [ 1 , ∞)是一個加權指數。式(8)需要滿足

模糊聚類就是通過迭代最優化目標函數Jc實現的, 這是一個進行優化的過程。其中模糊隸屬度uij和聚類中心cij分別為:

這個過程從一個隨機的聚類中心開始, 通過搜索目標函數的最小點, 不斷調整聚類中心和每一個樣本的模糊隸屬度, 達到確定樣本類別的過程。

加權指數m控制著聚類的模糊性。m越接近于1, 聚類越趨向于突變(crisp),m越大, 結果越模糊、相對更易于反映空間的漸變性, 但過大的m值將導致類別間的重疊太多, 聚類結構不清晰, 因此對m的選取需要在模糊度與清晰度的聚類結構間進行權衡。參考文獻[29]根據實驗建議最佳m位于區間 [1.5,2.5]之間。

3 實驗數據分析

根據前文建立的雙相-隨機介質模型, 利用高階有限差分技術正演計算海底沉積物的一次反射波信號。為簡化計算, 本文設計了一個雙層地質模型, 第一層(0~400 m)為海水, 第二層為海底沉積介質, 該部分介質為雙相-隨機介質, 彈性參數的隨機擾動由指數自相關函數給出, 見圖1。正演模型尺度為1 000 m×750 m,空間網格步長為 5 m, 時間步長為 0.1 ms, 所用Ricker子波主頻為30 Hz, 震源深度保持H= 300 m不變, 接受器置于炮點相同位置, 使之水平移動得到海底反射波的自激自收剖面, 每個模型通過正演模擬計算得到 100道地震記錄。海底沉積物底質顆粒的彈性參數見表1。利用等效介質理論計算的雙相介質的彈性參數見表2。根據Folk沉積物分類[30]方法, 在本文中, 泥質礫巖中礫巖: 泥巖=7∶3, 泥質砂巖中砂巖∶泥巖=1∶1。

圖1 由隨機函數產生的隨機擾動Fig.1 Random disturbance generated by the random function

表1 固體顆粒彈性參數Tab.1 Elastic property of solid particles

在得到海底反射波之后, 利用小波變換[31]提取該反射波的特征向量。提取步驟: (1)首先對信號進行Hilbert變換得到原信號的解析表達式; (2)對信號的解析表達式的實部和虛部分別作小波變換, 之后對某一級小波分解信號求模, 即為反射波的包絡特征向量。

表2 沉積物底質彈性參數Tab.2 Elastic property of marine sediments

綜合考量了輸入神經網絡的特征向量個數應盡可能少及特征向量要充分體現反射波特性時應盡可能多的特點, 經本文試驗計算后選擇對小波分解的第五層系數作求模運算。 本文選擇的小波基為harr小波。

有限差分正演計算的 3個模型的第 100道地震記錄見圖2。利用小波變換計算其包絡作為特征向量見圖3。本文在每個模型中選取20個一次反射波的特征向量作為支持向量機的訓練數據, 將其余數據作為預測數據, 同時利用算法對支持向量機要預測的數據進行分類以便對比這兩種方法的分類效果。從圖4中可見利用差分進化算法優化過的支持向量機的預測分類正確率達到了100%。圖5為模糊C均值聚類算法的聚類效果, 正確率也為100%。從圖4、圖5中可見, 支持向量機分類方法和模糊 C均值聚類分類方法均取得了極好的分類結果。

為了驗證本文所用方法的穩定性, 本文對反射波信號中分別加入10%, 30%, 50%的高斯白噪音, 某一道地震記錄及加噪后的波形見圖6。從圖中可見高斯白噪音幾乎全部掩蓋了原始信號的特征, 單純地依靠人的經驗已經無法區分海底沉積物底質類型了。本文仍然采用支持向量機和模糊 C均值聚類進行分類識別, 訓練樣本和預測樣本同上。此時, 本文以訓練樣本的最大分類正確率為目標函數, 利用差分進化算法對支持向量機分類參數進行優化選擇計算。最優化的支持向量機預測分類的結果與模糊 C均值聚類分類結果見表3。

圖2 三個模型第100道地震記錄對比Fig.2 Seismic record of trace No.100 in three models

圖3 三個模型反射波第100道的特征向量Fig.3 Characteristic vector of trace No.100 in three models

圖4 支持向量機預測分類結果Fig.4 The classification predicted by SVM

從表3中可見, 在信號中加一定的噪音之后,模糊C均值聚類相較于最優化支持向量機預測分類的正確率明顯較低。可見在低信噪比、小樣本情況下, 支持向量機仍然取得了比較好的分類效果。這充分說明了支持向量機對于分類數據具有很好的泛化能力和分類識別能力且具有較好的抗噪能力和魯棒性。

圖5 模糊C均值聚類3D聚類分類圖Fig.5 The clustering result produced by FCM

表3 加噪信號預測分類正確率Tab.3 The accuracy of using Optimized-SVM and FCM to classify the data with noise

從兩種方法實現的原理上, 本文分析了模糊 C均值聚類和支持向量機對于分類識別問題各自的優缺點: 模糊 C均值聚類相較于支持向量機不需要先驗信息, 直接根據聚類方法中的距離函數即可進行分類, 屬于無監督分類方法, 但對噪音數據比較敏感, 分類正確率不高; 支持向量機需要一定數目的先驗信息(訓練樣本), 屬于監督分類方法, 但支持向量機對數據泛化能力強, 預測分類正確率較高。

鑒于此, 本文提出一種混合分類方法, 即模糊C均值聚類-支持向量機分類方法, 該方法充分利用上述兩種方法的優點而達到無需先驗信息的高精度、快速分類。這種混合分類方法, 從算法的實現上也可以稱為兩步分類法。兩步分類法計算步驟如下:

步驟1: 數據初步聚類分析, 此步驟主要是通過模糊C均值聚類實現。具體操作如下:

(a) 確定分類數C0、m、精度要求等參數;

(b) 初始化cj(j= 1,2,… ,C0);

(c) 根據式(9)計算u和cj(j= 1,2,… ,C0);

(d) 根據式(8)計算模糊聚類目標函數, 判斷是否滿足聚類精度要求, 滿足則算法終止, 否則返回(c)。

步驟2: 支持向量機訓練樣本篩選, 根據模糊聚類的結果選擇最靠近每類中心的樣本作為支持向量機的訓練樣本; 首先計算每類的類內均值, 然后計算每類的類中所有樣本到中心值的距離矩陣, 從每類的距離矩陣中篩選出距離最小的若干個樣本構成支持向量機的訓練集;

步驟3: 在用訓練數據訓練支持向量機時, 利用差分進化算法對支持向量機中關鍵參數進行最優化搜索;

步驟 4: 用上述訓練好的支持向量機模型對其余數據進行預測分類輸出。

基于上述對海底沉積物數值模擬分類方法的論證, 可見本文所用方法是全完可行的。為了便于進一步對不同海底沉積物底質產生的地震反射波進行深入研究, 也為研究不同特征提取技術或其他分類識別算法對海底沉積物類型的識別能力, 本文歸納總結了利用計算機數值模擬技術對海底沉積物進行分類識別的一般化研究流程:

(1) 模型構制, 先將固體顆粒的彈性參數按照一定的體積百分比利用等效介質理論轉變為海底沉積物的彈性參數;

(2) 地震反射波信號采集, 基于雙相-隨機介質彈性波動方程, 利用有限差分技術計算模型的地震反射波信號;

(3) 特征向量提取, 利用小波變換或其他特征提取技術提取海底沉積物底質反射波的特征向量;

(4) 分類效果評價, 利用模糊 C均值聚類-支持向量機分類方法預測分類并進行評價。

此外, 鑒于本文所提方法具備較好的函數泛化能力及模糊 C均值聚類的無監督模式識別的特點,本文所采用的方法也有利于應對實際復雜的海底沉積物底質的分類識別, 這也是本文下一步的研究目標。

圖6 原始地震記錄與加噪音地震記錄波形對比Fig.6 The comparison of the original data and the data with noise

4 總結

海底底質類型進行分類研究一直是一個經典、熱點的課題。本文在總結了前人關于海底沉積物分類研究的基礎上, 率先提出采用計算機數值模擬手段來開展海底沉積物的分類識別研究。首先利用計算機數值正演技術模擬實際地震勘探數據采集過程,然后分別利用模糊 C均值聚類和基于差分進化算法優化的支持向量機對地震反射波進行分類識別, 再分析了上述兩種方法的優缺點之后, 本文提出了一種對海底沉積物分類識別的兩步法操作, 該方法有機地結合了模糊 C均值聚類的無監督特性和支持向量機的強泛化能力。在論證了本方法的可行性之后,本文歸納總結了一套利用計算機數值模擬技術進行海底沉積物分類識別的一般化流程以便進一步開展更加廣泛深入的研究。

[1]Mackenzie K V.Reflection of sound from coastal bottoms [J].Acoust Soe Am, 1960, 32(2): 221-231.

[2]Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid, I.Low-frequency range[J].Acoust Soc Am, 1956, 28(2): 168-178.

[3]Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid, II.high-frequency range[J].Acoust Soc Am, 1956, 28(2): 168-179.

[4]Biot M A.Mechanics of deformation and acoustic propagation in porous Media [J].Journal of Applied Physics, 1964, 33(4): 1482-1498.

[5]孟金生, 關定華.海底沉積物的聲學方法分類[J].聲學學報, 1982, 7(6): 337-343.

[6]王正垠, 馬遠良.寬帶聲吶湖底沉積物分類研究[J].聲學學報, 1996, 4: 517-524.

[7]劉建國, 李志舜.基于連續小波變換的湖底回波特征提取[J].西北工業大學學報, 2006, 1: 111-114.

[8]卜英勇, 張超, 聶雙雙.基于離散小波變換的水下回波信號尾波包絡特征提取[J].鄭州大學學報(工學版),2007, 4: 80-83.

[9]鄧躍紅, 聶雙雙.基于小波變換的水下超聲波測距方法研究[J].鄭州大學學報(工學版), 2007, 28(4):75-79.

[10]裴正林.雙相各向異性介質彈性波傳播交錯網格高階有限差分法模擬[J].石油地球物理勘探, 2006,41(2): 137-143.

[11]楊頂輝.雙相各向異性介質中彈性波方程的有限元解法及波場模擬[J].地球物理學報, 2002, 45(4): 575-583.

[12]Zhu X, Mcmechan G A.Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory[J].Geophysics, 1991, 56(3): 328-339.

[13]王東, 張海瀾, 王秀明.部分飽和孔隙巖石中聲波傳播數值研究[J].地球物理學報, 2006, 49(2): 524-532.

[14]姚姚, 奚先.隨機介質模型正演模擬及其地震波場分析[J].石油物探, 2002, 41(1): 31-36.

[15]奚先, 姚姚.隨機介質模型的模擬與混合型隨機介質[J].地球科學-中國地質大學學報, 2002, 27(1): 67-71.

[16]Li C P, Liu X W.Study on the scales of heterogeneous geologic bodies in random media[J].Applied Geophysics, 2011, 4: 363-369.

[17]殷學鑫, 劉洋.二維隨機介質模型正演模擬及其波場分析[J].石油地球物理勘探, 2011, 6: 862-872

[18]李紅星, 陶春輝.雙相各向異性隨機介質偽譜法地震波場特征分析[J].物理學報, 2009, 4: 2836-2842.

[19]Mavko G, Mukerji T, Dvorkin J.the rock physics handbook: tools for seismic analysis in porous media[M].Cambridge: Cambridge university press, 2003.

[20]Gassmann F.Uber die Elastizitat porous media[J].Vier der Natur Gesellschaft in Zurich, 1951, 96: 1-23.

[21]Gessrtsna J, Smit D C.Some aspects of elastic wave propagation in a fluid-saturated porous solids[J].Geophysics, 1990, 95: 15643-15656.

[22]Zeng Y Q, He J Q, Liu Q H.The applieation of the perfeetly matehed layer in numerical modeling of wave propagation in poroelastic media[J].Geophsics, 2001,66(4): 1258-1266.

[23]Haykin S.神經網絡原理(第二版)[M].葉世偉譯.北京: 機械工業出版社, 2004: 229-237.

[24]Cortes C, Vapnik V N.Support vector networks[J].Machine Learning, 1995, 20: 273-297.

[25]Vapnik V N.Statistical learning theory[M].New York:Wiley, 1998.

[26]Rainer S, Kenneth P.Differential evolution: A simple and efficient adaptive scheme for global optimization over continuous spaces [J].Global Optimization, 1997,11: 341-359.

[27]Bezdek J C.Pattern Recognition with Fuzzy objective Function Algorithms[M].New York: Plenum Press,1981.

[28]Bezdek J C, Hathaway R J.Recent convergence results for the fuzzy c-means clustering algorithm[J].Classifieation, 1988, 5(2): 237-247.

[29]高新波, 李潔, 謝維信.FCM算法中參數m的優選[J].模式識別與人工智能, 2000, 13(1): 7-11.

[30]Folk R L, Andrews P B, Lewis D W.Detrital sedimentary rock classification and nomenclature for use in New Zealand[J].New Zealand Journal of Geology and Geophysics, 1970, 13: 937-968.

[31]Daubechies I.小波變換十講[M].李建平, 楊萬年,譯.北京: 國防工業出版社, 2004: 56-99.