楊玉,吳云龍,單維鋒,石江浩
(1.防災科技學院,河北省地震動力學重點實驗室,河北 三河 065201;2.中國地震局地震研究所,中國地震局地震大地測量重點實驗室,湖北 武漢 430071; 3.北京航空航天大學,北京 100191)
歐洲空間局(ESA)于2009年3月17日發(fā)射了地球重力場和海洋環(huán)流探測衛(wèi)星[1](GOCE,Gravity field and steady-state Ocean Circulation Explorer),該衛(wèi)星以100km空間分辨率、1mGal的精度測定全球地球重力場,并以1-2cm的精度測定全球大地水準面[2]。GOCE衛(wèi)星搭載的高靈敏度重力梯度儀可對地球重力場每分鐘的變化進行三維測量[3],獲得的數據可更深入了解地球內部結構,觀察海洋和氣候變化。由于衛(wèi)星重力梯度在測量時受到衛(wèi)星系統(tǒng)周圍環(huán)境、檢校、定姿等問題的強烈干擾,造成衛(wèi)星重力梯度數據中不可避免地出現粗差[4]。因此在進行重力梯度數據的預處理時,必須進行相應的粗差探測與剔除來獲得數據保障。
粗差探測一直是國內外大地測量數據處理的熱點問題之一,許多專家學者提出不同角度的粗差探測方法。常見的一類方法是研究粗差數據套入模型中,如暴景陽等[5]基于逆?zhèn)鞑ド窠浘W絡對多波束測深數據中的粗差進行剔除,該方法用實驗模型對比分析擬合圖,驗證了逆?zhèn)鞑ド窠浘W絡方法的有效性;基于數學統(tǒng)計原理,閆廣峰等[6]人利用L1范數良好的抗差特性,在矩陣初等變換的理論上,求證了L1范數具有粗差定位的能力;為了豐富粗差探測理論,崔太岷等[7]基于半參數的方差膨脹模型,推導了Score檢驗公式,從而提高參數估計值接近真實值的精度。另一類方法是分析觀測值之間的相關性,如張建等[8]人利用基于卡方檢驗的粗差抗差算法減少觀測值由于相關性引起的粗差誤判,保證粗差定位的準確性;劉根友[9]模擬觀測值,采用最小二乘原理確定粗差;李玉芝等[10]人在不等權的條件下,以模擬觀測值和含粗差觀測值之間出現分群現象為目標,對粗差進行定位。
作為機器學習算法中的常用算法,聚類算法通過對衛(wèi)星數據進行數據挖掘和可視化分析,探究海量數據中有價值的特征信息。作為當前最流行、使用簡單廣泛、易于實現且計算效率較高的聚類算法,K-Means均值聚類算法對于大數據集,具備較小的時間和空間復雜度,該方法在地球科學領域中也得到驗證并使用。余德清[11]應用K-Means均值聚類算法對水面信息進行提取,通過算法本身的收縮運算,結合觀測數據進行對比,研究湖水資源量的變化。
本文基于K-Means均值聚類算法對數據進行聚類,分析了重力梯度數據之間的相似度,提出了基于K-Means均值聚類算法對衛(wèi)星重力梯度變量這樣的大數據集有高效的聚類效果,并設計了實驗,對比原始數據與聚類結果,驗證算法性能。K-Means均值聚類算法的研究未來還有更廣闊的空間,更多優(yōu)化算法會在粗差探測領域中發(fā)揮作用。
K-Means均值聚類算法屬于機器學習中的數據挖掘技術,它等價于數據集的分組,給定一個數據集,聚類算法可將數據集中的每個數據劃分成不同類型的組。理論上,同一組中的數據具有相似的屬性或特征,而不同組中的數據具有高度不同的屬性或特征。
典型的聚類算法內部計算過程,如圖1所示:

圖1 典型的聚類算法內部計算過程
K-Means均值聚類算法是一種基于迭代求解、距離的無監(jiān)督聚類分析方法,其運算速度快,操作簡單,常被應用在連續(xù)型的序列數據領域。K-Means均值聚類算法最關鍵的技術就是選取參數k,一般有兩種方法選擇最優(yōu)的k值,分別為:手肘法,輪廓系數法。
方法一:手肘法
手肘法的核心指標是誤差平方和(SSE,Sum of the Square errors)越小,即聚合程度越高。
(1)
當k小于最佳聚類中心個數時,由于k的增大會提高每個簇的聚合程度,故SSE的值會迅速下降;而當k接近最佳聚類中心個數時,隨著增加k值,聚合程度的效果不大,SSE下降的趨勢減小,k越增大而SSE越趨于平緩,即SSE和k的關系圖是一個手肘的形狀,隨著分類的類別數增加,這個肘部對應的k值就是數據的估計聚類數。實驗開始假設k>0,通過畫出K-SSE曲線,找出下降途中的拐點,即可較好地確定K值,曲線的第一個(或最顯著的)拐點表示估計的簇數。
方法二:輪廓系數法
a(i)是某樣本xi與同簇的其他樣本的平均距離,稱為凝聚度,b(i)是某樣本xi與最近簇中所有樣本的平均距離,稱為分離度。
首先計算樣本i到同簇其他樣本的平均距離a(i),a(i)越小,說明樣本i越應該被聚類到該簇。
然后計算樣本i到其他某簇的所有樣本的平均距離b(i),b(i)越大,說明樣本i越不屬于其他簇。
最后根據樣本i的凝聚度a(i)和分離度b(i),定義樣本i的輪廓系數,輪廓系數S(i)越大,聚類效果越好[12]。表達式(2)定義如下:
(2)
按照最鄰近原則把待分類樣本點分到各個簇。然后按平均法重新計算各個簇的質心,從而確定新的簇心。一直迭代,直到簇心的移動距離小于某個給定的值。
采用300階的EMG96模型,基于正常橢球GRS80模擬了徑向重力梯度觀測值Vzz,時長為1天,共 17 280個數據,將其視為“純凈”數據。表1是模擬衛(wèi)星重力梯度數據采用的相關參數,表2是數據中常見的統(tǒng)計值[13]。

數據模擬采用的相關參數 表1

模擬重力梯度數據的統(tǒng)計值(單位:10-9s-2) 表2
計算整個數據序列的標準差σ(如表2),以0為期望、0.01σ為標準差模擬生成含白噪聲的序列,再向白噪聲序列中加入169個的粗差,其中149個粗差以隨機分布的形式加入整個數據序列中;20個粗差以連續(xù)的形式加入數據序列,最終得到含有白噪聲和粗差的衛(wèi)星重力梯度數據。
在衛(wèi)星重力梯度數據粗差探測的模擬研究中,通常采用以下兩個指標來評價粗差探測的效果:
(1)成功率:
ORS=ns/no
(3)
式中:ns為探測到的粗差的個數,no為數據集中粗差的總個數。
(2)失敗率:
ORF=nf/n
(4)
式中:nf為有效數據被錯誤地探測為粗差的個數,n為總的數據個數。
需要指出的是,GOCE實測數據中的粗差個數是未知的。本文通過模擬計算和分析,利用上述兩個指標對粗差探測方法進行評價,可為GOCE實測數據的粗差探測提供參考依據。
K-Means均值聚類算法是將數據集按照樣本之間的距離大小,自動劃分為k組,選定k個初始的簇中心,按照以下步驟迭代細化[15](如圖2所示):

圖2 K-Means均值聚類算法實驗流程
訓練策略中,數據樣本處理步驟為:
(1)輸入數據集,選定好k個簇中心。隨機選擇數據集中的k個點分別作為k個簇的初始中心,然后計算其他點到這k個簇初始中心點的距離,從而決定這些其他點都屬于那個簇中心。比如某個點A到這k個點的距離產生k個數值,最小的那個值就是點A對應的簇;
(2)從數據集中隨機選擇k個分類對象作為初始化的k個聚類質心。簇中所有數據的均值通常被稱為這個簇的“質心”。在一個二維平面中,一簇數據點的質心的橫坐標就是這一簇數據點的橫坐標的均值,質心的縱坐標就是這一簇數據點的縱坐標的均值;
(3)將帶有聚類的數據放到一個聚類集合中,利用歐式距離計算相似度;
歐式距離表達式:

(5)
(4)根據聚類結果,不斷迭代,更新聚類質心,使得類內的相似度最大,類間的相似度最小。
(5)計算每個簇包含全部點的坐標平均值作為新的質心,進行2次迭代,直到誤差平方和最小,即結果趨于平穩(wěn)收斂。
K-Means均值聚類算法對于大數據集樣本處理具有較高的效率且是可伸縮性的,該算法的迭代優(yōu)化功能改善了初始監(jiān)督學習樣本識別不合理的地方。當樣本結果密集時,該算法的效果也較好,是其他聚類算法(如:譜聚法,高斯混合模型算法等)的基礎算法。
該方法的聚類階段如圖3,圖4所示(圖中x和y分別代表數據點的橫縱坐標值)。該聚類算法包含兩個部分:其一是k,k代表著類的數目,實驗中將數據聚為5類。其二是means,表示每次計算聚類中心的時候采取的是計算平均值。

圖3 原始數據

圖4 聚類結果(k=5)
從圖中可以看出,離群點是數據集中部分區(qū)域的粗差,在訓練時無法在聚類過程中統(tǒng)一到集合中。
如圖5是數據集在K-Means均值聚類模型下的粗差探測效果,圖中清晰地標注了粗差的位置,可以看出,聚類算法能夠準確地識別粗差。實驗中,將訓練所得數據與輸入的原始數據(含粗差數據)進行對比處理后,以此來初步判斷網絡模型的有效性。

圖5 在K-Means均值聚類模型下的粗差探測效果
K-Means均值聚類算法在重力梯度異常點檢測數據處理中具有顯著的優(yōu)勢。表3為應用K-Means均值聚類算法對重力梯度數據粗差探測的結果,圖6展示了應用K-Means均值聚類算法前后的粗差探測效果對比。結果顯示,經過模型訓練后的數據具有顯著的粗差識別能力。

重力梯度數據粗差探測結果 表3

圖6 重力梯度數據粗差探測前后的對比(片段)
實驗證明,基于K-Means均值聚類算法,在選定聚類中心個數(簇的個數)、質心、迭代次數后,數據與數據之間區(qū)別較明顯,且質心大小相近時,其聚類結果較理想。對于處理大數據集合,該算法非常高效,適應性較好。
本文基于K-Means均值聚類算法進行詳細的闡述,應用該方法對重力梯度變量的模擬數據進行粗差探測,結果顯示,K-Means均值聚類算法方法能夠高效地應用在衛(wèi)星重力梯度數據的預處理中。未來的研究工作可將改進的聚類算法與深度學習結合進行探測,以探索得到更適用于粗差探測的方法。