譜混合模型方法優化及其在海洋水團分析與水交換研究中的應用

宋　軍，劉玉龍，李　靜，2，郭俊如，牟　林，姚志剛，李希彬

（1.國家海洋信息中心，天津　300171；2.上海海洋大學　海洋科學學院，上海　201306；3.國家海洋局海洋減災中心，北京　100194；　4.中國海洋大學　教育部物理海洋學重點實驗室，山東　青島　266100；5.國家海洋局天津海洋環境監測中心站，天津　300457；6.國家海洋局空間遙感與應用研究重點實驗室，北京　100081）

?

譜混合模型方法優化及其在海洋水團分析與水交換研究中的應用

宋軍1，4，劉玉龍1，李靜1，2，郭俊如3，6，牟林1，姚志剛4，李希彬5

（1.國家海洋信息中心，天津300171；2.上海海洋大學海洋科學學院，上海201306；3.國家海洋局海洋減災中心，北京100194；4.中國海洋大學教育部物理海洋學重點實驗室，山東青島266100；5.國家海洋局天津海洋環境監測中心站，天津300457；6.國家海洋局空間遙感與應用研究重點實驗室，北京100081）

摘要：優化了譜混合模型（Spectral Mixture Model，SMM）分析方法，提出了譜混合模型方法中兩個關鍵參數的一般優化方案。優化后的方法能夠對大量的數據樣本進行快速聚類分析，并通過求解概率密度函數確定不同聚類之間的混合區域。以該方法在海洋水團以及水交換中的應用為例，詳細闡明了譜混合模型方法的工作原理及過程。在譜聚類方法基礎上建立的譜混合模型分析法，避免了傳統模糊聚類分析方法的不足，即使在物理量的散點數據分布呈現廣泛連續性時，仍然能夠抓住數據時空分布的主要變化方向，其在水團的辨別、水團邊界以及水交換混合區的分布及其變化規律的研究中具有廣泛的應用。關鍵詞：譜聚類；譜混合模型；海洋水團分析；水交換；模糊數學；溫度-鹽度圖

水團是源地與形成機制相近、具有大體相同的物理、化學和生物特征及變化趨勢，并與周圍水體有明顯差異的水體（張緒東等，2004）。不同的水團，其溫度、鹽度、密度等要素也不同，聲學、光學性質也有一定的差異，這些對于海軍潛艇的活動、水雷布設、水下通訊及監視，都有巨大的影響。在不同水團接壤、交匯的邊界水域，大多是有名的漁場。因此，對水團以及水團之間的交換混合區域給出準確的鑒別，并對其變化給出迅速的判斷和預測，能夠對軍事、漁業和水產事業與決策提供重要的信息保障。舉例來說，中國的近海，大部分地處中緯度溫帶季風區，四季交替明顯，季節變化顯著；深度不足200m的淺海，區域寬闊，島嶼棋布，岸線復雜；東部海域有強大的黑潮及其分支，西部有眾多的江河徑流入海（浦泳修，1981）。因而中國沿海水域水團分布復雜多變，尋找一種快速有效的水團鑒別方法具有重要的科學意義與實際意義。

Sverdrup等（1942）首次使用溫度-鹽度圖（T-S圖）進行水團分析，Miller等（1950）則首次應用T-S圖研究了水團的劃分問題。近年來，基于T-S圖的隸屬函數方法（Lietal，2004）與聚類分析方法（Su etal，1989）被廣泛用于研究水團的劃界與混合問題。但是隸屬函數方法只對特定的T-S曲線形狀有效，并且當使用同一方法應用到不同海域的時候，經常需要大量且復雜的調整工作才能使其適合研究的需要。經典的聚類分析方法有著更一般的應用方法，但當T-S分布呈現廣泛的連續性時，該方法很難抓住T-S曲線的主要變化方向。此外，以上兩種方法在處理大量數據點集的情況時，效率都非常低下。近年來，隨著觀測數據與數值模型結果數據的日益龐大，且可用于水團分析的要素從溫度、鹽度轉為多種要素，作為NP-hard問題的水團分析過程使經典的水團分析方法已不能滿足實際研究的需要。

鑒于以上現狀，宋軍（2010）將圖論中的最新研究成果——譜聚類方法（Spectral Clustering，SC）引入到海洋科學的研究當中，并基于此構建了可解決上述傳統水團分析方法問題的譜混合模型方法。本文將對該方法的構建過程與應用方案進行詳細的描述，并對模型中的兩個關鍵參數給出了具體的優化方案。

1算法綜述

為了更清楚地解釋這一方法，我們首先定義集合P = {P1，…，Pn}為所有目標研究區域的數據點，其中n為所有數據點的數量。集合M = {m1，…，mk}則被定義為所有聚類的集合，以滿足我們將上述n個數據點（P）劃分為k個聚類（M）的描述需要（Song，2011）。應該說明的是，為使本研究具有實際的物理意義，聚類的數目k應該小于數據點的數目n，并且每個點有且只屬于一個聚類。譜混合模型的計算過程主要包括三大步：第一步，使用譜聚類方法將所有目標數據點劃分為任意需要數目的聚類。所有的聚類以及其重心位置將在這一步得到確定；第二步，基于我們構造的概率密度函數，將得到每一個數據點Pi∈P映射到所有參與相互混合作用的聚類mi?M的隸屬概率；第三步也是最后一步，任意兩個屬于上一步考慮范圍的聚類之間的交換區將得到定義，并且以此為基礎，任意多個聚類的混合區域也將得到定義。不同于經典的混合模型，由于譜混合模型基于譜聚類方法，其基于數據點的連續程度考慮了數據點的主要分布方向，無論這一方向是直線的還是任意曲線的。此外，在計算概率密度函數時，譜混合模型還考慮了每一個聚類的勢力范圍，而并非僅使用各聚類重心的位置作為計算混合區域分布和位置的標準。在分析水團混合的研究中，各聚類的勢力范圍是指其擁有的數據點的數目以及這些數據點的平均分布密度。在本研究當中，由于整體數據點的分布密度差異不大，取各聚類擁有的數據點的數目作為其勢力范圍的指標。

1.1譜聚類

譜聚類（Spectral Clustering，SC）是一個21世紀剛剛發展起來的最新的聚類方法，但其明顯的優越性使其本身成為最熱門的研究學科之一（Von，2007）。譜聚類基于圖像分析理論，在計算機科學，特別是人工視覺系統（Jitendra et al，2001）和機器學習（Shietal，2000）等人工智能方面的研究中，甚至在與計算機結合緊密的最前沿的其他學科，比如生物化學（Menschaert et al，2009）的研究中，都有著廣泛的應用。由于譜聚類方法不僅考慮了所有目標數據點之間的差異性，而且考慮了整體數據的分布結構及其連通性，使其更容易抓住主要矛盾，該方法表現出較高的計算效率，這一點在數據量（包括數據點的數目和參數空間亦即區分不同數據點所使用的性質參數的數目）較大的時候更加明顯（Andrew etal，2002）。更重要的是，譜聚類可以更容易地得到全局最優解，從而避免了經典聚類分析方法使結果在凸型數據區域的計算中陷入局部最優的問題（Luo etal，2003）。簡單來說，譜聚類方法將所有的目標數據點Pi∈P考慮為一個多維無向連通圖G的頂點（如圖1所示）。在本研究中，頂點來自模型結果中對應的每個計算網格的水質點數據。G是一個無方向但各邊具有不同權重的圖（加權圖），我們使用一個半正定的，對稱權重的矩陣A（Aij=Aji≥0）來表示任意兩點Pi和Pj之間的相似性，這個相似性也即是圖G中各邊的權重。這里我們使用Aij來表示兩個數據點Pi和Pj在其參數空間下的距離。舉例來說，邊A23表示頂點P2與P3之間的相似度。在圖G （P，A）中，求解聚類分析的問題可以等價為求圖G的最小‘割權’問題。

圖1　用于表示各點之間相似度的無向連通加權圖G（P，A）。其中A中各邊的權重（圖中由各邊的粗度表示）表示兩個頂點之間的相似度大小

為了避免任意孤立點從圖中被分離出去，如大多數的相關研究一樣，這里也采用一個升級的割權定義來代替上面的割權定義，稱這一升級的割權為標準化割權（Ncut），其被定義為如下形式：

其中，vol（mi）=∑i∈miAij。為求得上述標準化割權的最小值，譜聚類算法是目前為止被證明的最直接和有效的方法（Von，2007）。對于屬于集合Rl（字符R通常被用作表示實數集，其上標則表示此實數集的維數）的目標數據點集合P，將P劃分為k個聚類，其完整的譜聚類算法被簡單的描述如下：

1）構造親和度矩陣A∈Rn×n，在本研究中，我們使用高斯相似函數來定義A：

其中，σ為尺度參數，當計算pi與pj兩點間的相似度Aij時，它決定了隨兩點間距離的增加而造成的衰減率。雖然這里有很多方法可以自動的優化參數σ（Von Luxburg，2007），但在本研究中，由于所有的數據都已進行了標準化處理，所以該參數被簡單地設為1。

2）定義一個對角矩陣D，其中Dii等于矩陣A中第i行所有數值的總和：

3）定義一個拉普拉斯矩陣L，這是一個半正定的矩陣：

4）求解矩陣L，得到最小的k個特征值（λ1，…，λk）和這k個特征值對應的特征向量（，…，）；

5）構造目標矩陣O=（（x）1，…，（x）k）∈Rn×k，其每一列分別等于上一步求得的各個特征向量；

6）對矩陣O的每一行分別做標準化，使其擁有統一的單位標準，并由此構造出矩陣Q∈Rn×k如下：

7）將矩陣Q的每一行看作每個對應的數據點，這就意味著矩陣Q是一個有N個數據點，并且每個數據點有k個參數屬性的矩陣。最后使用K-means或其他經典的聚類分析方法對這一新構造的數據點集進行聚類分析，進而得到聚類集合M以及其相應的每個聚類和聚類的重心。

在本研究中，使用經典的K-means聚類分析方法來完成最后的一步運算。簡單來說，K-means聚類分析的核心算法可以描述為以下5步。如果希望提高其在計算機中的計算速度，可以參考Lai等（2009）發展的快速K-means聚類計算方法。

1）隨機地給出任意k個聚類的重心坐標；

2）分別計算每個數據點到這些聚類重心的距離。本研究采用歐式相空間距離公式來計算這一距離；

3）將每個點歸入距離其最近的聚類重心所代表的聚類中，這意味著第一次將所有的點分為了k個聚類；

4）通過計算每個聚類中所有數據點坐標的平均值，重新計算每個聚類的重心坐標；

5）重復做以上第2步到第4步的運算，直到所有的聚類集合都不再發生變化。至此，就得到了最終的分類方案和每個聚類的重心。

為了更直觀地展現傳統的K-means聚類分析法與譜聚類分析法在實際研究中的不同，本文分別用兩種方法對黑潮與東中國海地區的水團進行了聚類分析（如2所示）。圖2清楚地顯示出直接采用K-means聚類分析方法（圖2a）和譜聚類方法（圖2b）的一個結果對比圖。很顯然后者更能表現出黑潮水與東中國海陸架水的水團分布情況，從而能夠定義一個兩者之間的水交換帶。該結果主要是由于在東中國海，垂直于黑潮流向的溫度鹽度變化趨勢要遠大于沿著黑潮流向方向的變化趨勢造成的。更多類似于K-means方法與譜聚類方法在處理不同數據集的結果對比，可以參照Cai等（2005）的工作，其研究結果表現了譜聚類方法相對于傳統聚類方法的絕對優勢。

圖2　兩種聚類方法對經過標準化的100 m層的溫度-鹽度（T-S）數據進行聚類分析的結果圖

2.2概率密度函數

為了定義一個模糊混合區域，我們首先需要知道每個數據點多大程度的隸屬于每個需要考慮參與交換的聚類。基于上面譜聚類分析得到的每一個聚類的重心坐標，將構造一個概率密度函數以計算每個數據點投影到相應的所有聚類的概率。具體算法構造如下：

1）對任意數據點Pi∈P，定義disi，kk（kk∈[1 k]）為數據點Pi到聚類mkk重心Ckk的距離（或者說Pi與mkk重心Ckk的差別）。由于所有參與計算的數據都已經進行了標準化轉換，所以這里采用歐式距離來計算disi，kk：

2）使用disi，kk的倒數構造一個變量wi，kk，以指出數據點Pi隸屬于聚類mkk的隸屬度權重。

如果?disi，kk（kk∈[1 k]）≠0，則：

這里，numkk是聚類mkk所包含的數據點數。

如果?disi，j（j∈[1 k]）=0，則：

3）構造混合分布矩陣H。使數據點Pi相對于聚類mkk的隸屬度權重作為計算隸屬度的權重wi∈P，kk∈M，則Hi，kk最終被定義為：

2.3對任意數量指定聚類混合區的定義

設T= {t1，…，tv}（ti-mj=?，1< v≤k）為我們所關注的聚類的集合，其是全部聚類集合m的部分或全部。為了定義這些聚類的混合區，首先設集合Rm，n與集合Sm，n是T中任意兩個聚類tm與tn（其中m≠n）在相空間上一個交集范圍的最大限度：

這里，Zm，n包括了聚類tm和tn所有位于混合區的的數據點的集合，如圖3所示。

圖3　譜混合模型在v=2，α=0.5，β=0.8，num1 = num2情況下，變量和集合之間相互關系差異的分布結構圖

為了方便地將這一概念在任意多個聚類參與的情況下進行描述，現定義∩（ai）：=a1∩a2…∩a1。設Z為位于聚類集合T中所有聚類混合區域的數據點集，則最終Z可以表示為：

這里值得指出的是，譜混合模型在水團分析和鋒面分析中具有潛在的以及非常深遠的應用前景。在此基礎上，定義出兩個水團的信息過渡區，并且求解所有信息過渡區的交集。此外，該方法還可用來定義任意多個水團之間的信息過渡區。圖4即是一個將集合數目擴大到3個而求解其信息過渡區的實例，實例中，v=3，α=0.5，β=1.0。圖中黑色粗線包裹的區域限定了三個聚類a，b和c混合區的最大界限。其中每個聚類所包含的數據點的數量已在各圖中各聚類名稱后面的括號中給出。

綜上，譜混合模型是在譜聚類方法的基礎上，通過構建概率密度函數得到每個數據點隸屬于某個聚類的隸屬度權重，然后引入參數α和β，創造任意兩個聚類最大限度的交集集合R和S，從而最終確定混合區（在水交換研究中即水交換區）數據點的集合為。其中參數α和β決定了混合區的空間位置和范圍，其中α表示數據點隸屬于不同聚類權重的差值，越大說明該數據點越偏離某一聚類；β則用于排除兩水團的外邊緣區域帶來的誤差。

圖4　譜混合模型在v=3，α=0.5，β=1.0情況下，變量和集合之間相互關系差異的分布結構圖

2.4參數α和β的優化

譜混合模型中的參數α和β，能夠保證將不同聚類之間的混合區域限制在各聚類中心連線的中點附近，對混合區的確定起到了關鍵性的作用，所以找到一種可以優化兩參數的方法是十分有必要的。本文對此提出了一套能夠優化這兩個關鍵參數的一般方案，具體闡述如下：

若要研究的是混合區的時空變化情況，則在滿足混合區域被控制在各聚類之間的前提下，其變化越顯著，越有利于研究的進行。為了方便敘述，這里以某一事件的發生，能夠引起混合區的變化為例，對優化參數的方法步驟進行闡明。首先定義4個指數IRB、IRA、RA以及TR，其中IRB為事件發生時相對于事件發生前混合區的變化率；IRA為事件發生時相對于事件發生后混合區的變化率；RA為平均氣候狀態下混合區占總區域的比例；TR則是以上3個指數的總和。

圖5顯示了在不同參數的配置下，在一個案例中，4個指數的變化分布情況。

首先判斷該事件的發生對混合區的影響是加強還是減弱。若有加強作用，則IRB、IRA、RA以及TR 4個指數必須均為正數；反之，指數IRB和IRA為負值；然后，考慮到參數α表示數據點隸屬于不同聚類權重的差值，其值越小，所確定的混合區越趨向各聚類中心連線的中點；最后，考慮到混合區的顯著性體現，即在滿足混合區被控制在各聚類之間的前提下，其占據總區域的比例越大，且在事件前后的變化越明顯，越有利于研究的進行。

基于以上3個原則，結合圖5中a、b、c、d4張圖，最終可以確定該案例中的最佳參數配置為α=0.15，β=0.25。

圖5　指數IRB、IRA、RA以及TR在不同參數配置下的分布圖

此外，若只研究各聚類之間混合區的空間分布情況，而不考慮其隨時間的變化，則上述方法中的4個指數只有RA有效，而RA反映的是參數α和β的取值越大，混合區占總區域的面積越大，這顯然不能以此確定兩個參數的取值，因為該混合區還需要滿足位于各聚類之間這一前提條件。故針對這類案例，該參數優化方法并不適用，而是要根據案例的實際情況來對這兩個參數進行優化。

3討論

完善后的譜混合模型能夠快速有效地分析水團以及確定水團之間的交換區，但模型中有幾點還有必要做進一步得改進。本研究中簡單的使用聚類mkk所包含的數據點的數量numkk來表示mkk的勢力范圍，這個定義是基于所有聚類中數據點的分布密度都是一致的這一前提假設。因此，在對各聚類的數據點分布密度差異較大的情況研究中，建議使用unmkk/ddkk來代替公式（8）中的numkk，這里ddkk表示聚類mkk數據點的平均分布密度。此外，與大多數聚類分析研究一樣，聚類個數k在本研究中是給定的。實際上，在更一般的應用當中，最優化的k值可以通過分析求解矩陣L（方程5）得到的特征值序列而得到。一個基本的想法是將這些特征值從小到大排列起來，然后自動確定數目k，使前k個特征值的和相對于所有特征值的和的比例小于某一標準。

參考文獻

Andrew Y N，Michael I J，Yair W，2002.On spectral clustering：Analysis and an algorithm.Advances in Neural Information Processing Systems，14（1-2）：849-856.

JitendraM，SB，TL，etal，2001.Contourand textureanalysis for image segmentation.Int JComput，43（1）：7-27.

Lai J，Huang T J，Liaw Y C，2009.A fast k-means clustering algorithm using cluster center displacement.Pattern Recogn，42（11）：2551-2556.

Li F，Jun X，Yao L，2004.New methods of fitting the membership function of oceanic water masses.Journal of Ocean University of China（English Edition），3（1）：1-9.

Luo B，W ilson R C，Hancock E R，2003.Spectral clustering of graphs.Computer Analysis of Images and Patterns，Proceedings，2756：540-548.

Menschaert G，Vandekerckhove T T M，Landuyt B，et al.，2009.Spectral clustering in peptidomics studies helps to unravel modification profile of biologically active peptides and enhances peptide identification rate.Proteomics，9（18）：4381-4388.

Miller A R，1950.A study of mixing processes over the edge of the continentalshelf.J.Mar.Res.，9（2）：145-160.

Shi J，Malik J，2000.The 20th anniversary of the IEEE transactions on pattern analysis and machine intelligence.Ieee T Pattern Anal，22 （1）：1-3.

Song J，Xue H J，Bao X W，W U Dexing etal.，2011.A spectralmixture model analysis of the Kuroshiovariability andthe water exchange between the Kuroshioand the East China Sea.China Journal of Oceanology andLimnology，29（2）：446-459.

Su Y，Yu Z，Li F，1989.The application of the cluster analysis on the shallow water mass analysis.Periodical of Ocean University of China，15（3）：24-28.

Sverdrup H U，Johnson M W，Fleming R H，1942.The oceans，their physics，chemistry and general biology.New York：Pretince-Hall，1-1087.

Von L U，2007.A tutorial on spectral clustering.Stat.Comput，17（4）：395-416.

浦泳修，1981.關于東海北部表層的水系和環流.海洋通報，（5）：10-16

宋軍.水交換模型的理論方法及應用研究：（博士學位論文）.青島：中國海洋大學，2010.

張緒東，張國友，佟凱，等，2004.黑潮源區海域水團分析.海洋通報，23（1）：15-21.

（本文編輯：岳心陽）

Optim ization of the spectralm ixturem odelm ethod and its app lication in theanalysisofwaterm assand water exchange

SONG Jun1，4,LIU Yu-Long1,LIJing1，2,GUO Jun-Ru3,6,MU Lin1,YAO Zhi-Gang4,LIXi-Bin5
（1.NationalMarineDataand Information Service,Tianjin 300171,China;2.CollegeofMarine Science,ShanghaiOcean University,Shanghai201306,China;3.NationalMarineHazard Mitigation Service,SOA,Beijing100194,China;4.MinistryofEducation Key LaboratoryofPhysicalOceanography,Ocean UniversityofChina,Qingdao266100,China;5.Tianjin Marine EnvironmentalMonitoring CentralStation,SOA,TianjinMarine EnvironmentalMonitoringand ForecastingCenter,Tianjin 300451,China;6.Key LaboratoryofSpace Ocean Remote Sensingand Application,Beijing100081,China）

Abstract：The method of spectralmixture model is optimized by determining its two parameters in a better way.The optimized method can make rapidclustering analysis of a large number of data samples.By solving the probability density function，itcan determine themixing areabetween differentclusters.W ith the application of themethod in the ocean water mass and water exchange as an example，this paper illuminates the working principle and process of spectralmixture modelin detail.The analysismethod ofspectralmixturemodelwhich isbased on the spectral clusteringmethod，avoids the inadequacy of thetraditional fuzzy clustering analysis.Even if thescattered data distribution ofphysicalquantity presents a wide range of continuity，it is still able to seize the main changedirectionof the time and space distribution of the data.Hence，ithasawide rangeofapplications in termsof the discrimination ofwatermass，the distribution and change rulesof watermassboundary and themixed areaofwaterexchange.

Keywords：SpectralClustering;SpectralMixtureModel;watermassanalysis;waterexchange;FuzzyMathematics;T-S diagram

通訊作者：郭俊如，博士，電子郵箱：874623647@qq.com。

作者簡介：宋軍（1983-），男，博士，主要從事近海動力學、業務化海洋學方面的研究。

基金項目：國家自然科學基金（41206013；41376014；41430963；41206004）；教育部物理海洋重點實驗室開放基金；2011年度高等學校博士學科點專項科研基金（20110132130001）；海洋公益性行業科研專項（201205018；201005019）；國家科技支撐計劃（2014BAB12B02）；天津市科技支撐計劃（14ZCZDSF00012）；國家海洋局青年科學基金重點項目（2012202；2013203；2012223）；國家建設高水平大學公派研究生項目（留金出[2008]3019;[2012]3013）；國家海洋向空間遙感與應用研究重點實驗室開放基金重點課題（201601003）。

收稿日期：2015-06-10；

修訂日期：2015-07-03

Doi：10.11840/j.issn.1001-6392.2016.01.010

中圖分類號：P733

文獻標識碼：A

文章編號：1001-6932（2016）01-0074-07