基于無監(jiān)督機(jī)器學(xué)習(xí)的噪聲信號(hào)聚類分析
——以郯廬斷裂帶濰坊段短周期密集臺(tái)陣觀測(cè)為例

楊勇剛， 鈕鳳林,2*

1 中國石油大學(xué)(北京)非常規(guī)油氣科學(xué)技術(shù)研究院， 北京 102249 2 萊斯大學(xué)地球科學(xué)學(xué)院， 美國休斯敦 77005

0 引言

現(xiàn)代先進(jìn)的地震計(jì)具有較寬的動(dòng)態(tài)范圍，可記錄來自一系列震源的震動(dòng)信號(hào)，包括與板塊構(gòu)造運(yùn)動(dòng)相關(guān)的信號(hào)，如地震事件(Ross et al., 2017)、火山噴發(fā)和非火山活動(dòng)震顫信號(hào) (Eibl et al., 2015; Shelly et al., 2007)，以及非構(gòu)造活動(dòng)的震源，如海水波動(dòng)(Gerstoft and Tanimoto, 2007)、風(fēng)(Johnson et al., 2019; Tanimoto and Valovcin, 2015)、公路和鐵路交通產(chǎn)生的地面運(yùn)動(dòng)信號(hào)(Riahi and Gerstoft, 2015; 劉磊和蔣一然, 2019; 李麗等, 2004; 徐善輝等, 2017).本文將由非構(gòu)造活動(dòng)引起的地面震動(dòng)信號(hào)統(tǒng)稱為地震噪聲信號(hào).

近十多年來，許多先進(jìn)的算法被提出來用于微地震或者震顫事件的檢測(cè)(Aguiar and Beroza, 2014; Meng and Ben-Zion, 2018b; Reynen and Audet, 2017; Ross et al., 2018; Yoon et al., 2015; 趙明等, 2019)，但這些算法在實(shí)際應(yīng)用時(shí)也面臨諸多挑戰(zhàn)，如在某些觀測(cè)環(huán)境下，地震噪聲信號(hào)的波形和微地震事件非常相似，甚至噪聲信號(hào)的振幅可能大于2級(jí)以下的微地震事件信號(hào)(Johnson et al., 2019; Withers et al., 1996)，由地震噪聲引起的地震事件的誤檢可能會(huì)對(duì)我們理解地震活動(dòng)機(jī)理產(chǎn)生較大的影響(Inbal et al., 2015; Yang et al., 2021).進(jìn)入21世紀(jì)以來，由地震噪聲場獲取的面波經(jīng)驗(yàn)格林函數(shù)被廣泛應(yīng)用于地下結(jié)構(gòu)成像或者動(dòng)態(tài)監(jiān)測(cè)，而經(jīng)驗(yàn)格林函數(shù)的計(jì)算精度取決于噪聲源空間分布是否相對(duì)均勻(Bensen et al., 2007; Shapiro et al., 2005)，因此了解噪聲的特征與噪聲源的時(shí)空分布對(duì)計(jì)算經(jīng)驗(yàn)格林函數(shù)至關(guān)重要.現(xiàn)在廣泛使用的以微地震識(shí)別與去噪為目的開發(fā)的監(jiān)督類機(jī)器學(xué)習(xí)算法，需要正確標(biāo)注的地震事件和噪聲信號(hào)作為標(biāo)簽進(jìn)行學(xué)習(xí)訓(xùn)練獲取可靠分類器或回歸函數(shù)(Johnson et al., 2019).因此，為了正確處理地震記錄，有必要對(duì)噪聲信號(hào)的組成及其可能的物理來源開展詳細(xì)的研究.然而，傳統(tǒng)的噪聲信號(hào)分類主要聚焦于簡單的時(shí)間(Johnson et al., 2019)、頻率域特征分析(Díaz, 2016; Johnson et al., 2019; Meng and Ben-Zion, 2018a; Riahi and Gerstoft, 2015; Suriach et al., 2005)，或者統(tǒng)計(jì)散點(diǎn)圖分析(Meng et al., 2019).針對(duì)密集臺(tái)陣記錄的海量連續(xù)波形數(shù)據(jù)，近些年發(fā)展起來的機(jī)器學(xué)習(xí)算法為我們研究噪聲信號(hào)提供了非常有力的手段.Johnson等(2020)針對(duì)一個(gè)布設(shè)于郊區(qū)位置的單分量密集臺(tái)陣數(shù)據(jù)，利用無監(jiān)督聚類算法識(shí)別出了幾種與風(fēng)擾動(dòng)灌木叢或建筑物有關(guān)的噪聲信號(hào)；Snover等(2020)則針對(duì)布設(shè)于城市區(qū)域的密集臺(tái)陣數(shù)據(jù)，采用深度自動(dòng)編碼器來識(shí)別與城市活動(dòng)有關(guān)的噪聲信號(hào).

機(jī)器學(xué)習(xí)是一種數(shù)據(jù)驅(qū)動(dòng)的技術(shù)，大部分機(jī)器學(xué)習(xí)算法可被劃分為監(jiān)督學(xué)習(xí)(supervised learning)與無監(jiān)督學(xué)習(xí)(unsupervised learning)類算法.前者通過分析大量標(biāo)簽(Label)數(shù)據(jù)中的特征，建立數(shù)據(jù)與模型之間隱含的非線性關(guān)系(Bishop, 2006).而后者則無需標(biāo)注樣本，通過算法自身來挖掘數(shù)據(jù)與模型之間的潛在關(guān)系(Mousavi et al., 2019; Yoon et al., 2015).鑒于我們事先未知噪聲信號(hào)的組成或者標(biāo)簽，所以無監(jiān)督類方法更適合本論文的研究問題.

無監(jiān)督學(xué)習(xí)算法中聚類和降維是兩種比較常用的簡單算法，其原理是根據(jù)相似性分別將輸入樣本進(jìn)行聚類以及對(duì)輸入數(shù)據(jù)進(jìn)行降維處理，以獲取樣本的分布規(guī)律與主要特征.常見的聚類算法有K-means、DBSCAN(density-based spatial clustering of application with noise)等算法.其中，K-means算法按照聚類中心與樣本之間的距離對(duì)樣本進(jìn)行聚類，最終使得簇內(nèi)的樣本足夠緊密，簇間樣本足夠分離(Jin and Han, 2017).K-means算法的聚類效果依賴于聚類簇中心的初始化，同時(shí)需要在聚類前對(duì)聚類簇個(gè)數(shù)進(jìn)行假定.而DBSCAN聚類算法則通過在特征空間中不斷尋找樣本點(diǎn)分布密集的區(qū)域來對(duì)樣本進(jìn)行聚類成簇(Ester et al., 1996)，因此不需要先驗(yàn)地設(shè)定簇的個(gè)數(shù)，但缺點(diǎn)是相對(duì)于K-means算法而言運(yùn)行速度較慢.常用的降維算法包括主成分分析法(Principal Component Analysis, PCA)、t-SNE(t-distributed stochastic neighbor embedding)等算法.其中，PCA算法的原理是通過分析樣本的特征，去除噪音與不重要的特征，保留最重要的特征，以實(shí)現(xiàn)提升數(shù)據(jù)處理的速度.處理上相應(yīng)地將n維特征映射到k(

利用布設(shè)于華北東部地區(qū)郯廬斷裂帶濰坊段的密集臺(tái)陣觀測(cè)數(shù)據(jù)(圖1)，我們使用便捷高效且被廣泛應(yīng)用的K-means算法(Jain, 2010; Wu et al., 2008; Xu and Wunsch, 2005)來聚類分析不同類別的地震噪聲信號(hào).此臺(tái)陣布設(shè)于郯廬斷裂帶上的一個(gè)地震空區(qū)位置，目標(biāo)是獲取地震空區(qū)內(nèi)的微地震活動(dòng)性與郯廬斷裂帶下方的精細(xì)結(jié)構(gòu).因此，本論文關(guān)于此區(qū)域噪聲信號(hào)特征的研究對(duì)空區(qū)內(nèi)的微地震識(shí)別與檢測(cè)有著至關(guān)重要的作用.

圖1 郯廬斷裂帶濰坊段密集觀測(cè)臺(tái)陣共計(jì)300個(gè)短周期臺(tái)站，臺(tái)站平均間距約5 km. 臺(tái)陣北部區(qū)域城市化程度較高，南部區(qū)域以山區(qū)為主，圖中道路交通線密集的范圍代表城鎮(zhèn)區(qū)域. 本文定義列的索引為[0 1 2 3 4 5 6 7 8 9 A B C D E F G H I J]，行的索引為[0 1 2 3 4 5 6 7 8 9 A B C D E]，臺(tái)站的名稱由行和列的索引決定，如臺(tái)站SPMN表示位于第M行N列位置處的臺(tái)站.Fig.1 Dense array deployed at the Weifang segment of the Tanlu FaultIt includes 300 seismic stations with a station average spacing of ～5 km. The northern part of array covers more urbanized areas, while its southern part is mainly rural areas with hills and mountains, and areas with dense road lines correspond to major cities and towns. The column index is [0 1 2 3 4 5 6 7 8 9 A B C D E F G H I J], and line index is [0 1 2 3 4 5 6 7 8 9 A B C D E]. The station is named according to its line and column indexes, such as the station SPMN located at the line M and column N.

1 臺(tái)陣概況和數(shù)據(jù)選取

我們使用的數(shù)據(jù)來源于一個(gè)布設(shè)于郯廬斷裂帶濰坊段的短周期密集臺(tái)陣(以下簡稱密集臺(tái)陣)，臺(tái)陣的地理坐標(biāo)范圍為36.2°N—37.0°N，118.5°E—119.9°E，包含300個(gè)短周期臺(tái)站(EPS便攜式數(shù)字地震儀，頻帶范圍為5 s～200 Hz)，臺(tái)站平均間距約為5 km(圖1)，以每秒200個(gè)采樣點(diǎn)(200 Hz) ，從2017年8月至10月，連續(xù)采集了3個(gè)月的三分量波形記錄，我們只選取了運(yùn)行穩(wěn)定的9月份的垂向分量連續(xù)波形記錄進(jìn)行分析.密集臺(tái)陣以濰坊市為中心，同時(shí)覆蓋其它較大城市，包括昌邑市和壽光市.總體上，人口分布密集的城市區(qū)域主要位于臺(tái)陣的北部，而南部則以農(nóng)村或山區(qū)為主.此密集臺(tái)陣為我們?cè)诔鞘?農(nóng)村混合的復(fù)雜環(huán)境下，研究地震噪聲信號(hào)的識(shí)別提供了非常好的契機(jī).

圖2 數(shù)據(jù)處理流程Fig.2 Workflow of data processing

2 數(shù)據(jù)處理

數(shù)據(jù)處理過程如圖2所示，我們首先對(duì)連續(xù)波形進(jìn)行1 Hz的高通濾波處理，然后，從國家地震數(shù)據(jù)中心(http:∥10.5.160.18/)收集2017年9月1日至30日時(shí)間范圍內(nèi)的地震目錄，共包含187個(gè)近震事件，9個(gè)M4+的區(qū)域地震事件(震中距小于2000 km)，和99個(gè)M5+的遠(yuǎn)震事件(震中距大于2000 km).針對(duì)這些地震事件，我們利用Taup程序(Crotwell et al., 1999) 計(jì)算它們?cè)诟鱾€(gè)臺(tái)站上的直達(dá)P波到時(shí)，并以此到時(shí)為中心，將-60～+60 s時(shí)間范圍內(nèi)的記錄直接從連續(xù)波形中剔除.在這里我們假設(shè)區(qū)域地震和遠(yuǎn)震中的高頻橫波和面波成分經(jīng)過長距離的衰減后可以忽略不計(jì).然后，我們將剩余的連續(xù)波形裁剪分段，每個(gè)波形片段的長度為2 s，共計(jì)50,858,207個(gè)波形片段數(shù)據(jù)(波形樣本).

針對(duì)波形樣本，我們首先進(jìn)行了去除線性趨勢(shì)和均值的預(yù)處理，然后，借鑒以往的分析波形特征來識(shí)別地震信號(hào)的思路(Küperkoch et al., 2012)，對(duì)高頻噪聲信號(hào)(1～100 Hz)的樣本波形進(jìn)行相應(yīng)的特征提取并根據(jù)這些特征進(jìn)行分類.我們提取了基于信號(hào)能量(Hulbert et al., 2019; Rouet-Leduc et al., 2019)和頻譜(Hammer et al., 2012; Mousavi et al., 2016)的7個(gè)特征(Johnson et al., 2020)(表1)來處理噪聲信號(hào)的分類問題，這些特征也常被用于地震事件檢測(cè).其中，波形的平方積分代表的是波形片段的能量特征.平方積分和最大振幅特征可以被用來區(qū)分具有不同震動(dòng)強(qiáng)度的噪聲信號(hào).譜峰頻率和中心頻率特征描述了噪聲信號(hào)的主要頻率成分.過零率和譜峰率來自于語音信號(hào)處理領(lǐng)域.其中，過零率指每幀信號(hào)中，信號(hào)通過零點(diǎn)(本文指向上過零率，對(duì)應(yīng)于振幅從負(fù)變?yōu)檎?的次數(shù).而譜峰率指信號(hào)峰值出現(xiàn)的次數(shù).這2個(gè)特征在語音識(shí)別和音樂信息檢索領(lǐng)域得到了廣泛的應(yīng)用，主要用于敲擊聲音的分類(如Koniaris et al., 2010; Moon et al., 2007).

表1 描述地震噪聲的7種數(shù)據(jù)特征Table 1 List of 7 data features employed in characterizing seismic noise

我們進(jìn)一步將樣本數(shù)據(jù)集以9∶1的比例劃分為訓(xùn)練樣本集和測(cè)試樣本集，數(shù)量分別為45772386和5085821個(gè).為了保證樣本各個(gè)特征處于同一量級(jí)，我們對(duì)各個(gè)特征采用標(biāo)準(zhǔn)化處理，以確保每個(gè)特征的均值為0，方差為1.訓(xùn)練數(shù)據(jù)集的特征散點(diǎn)圖(一種數(shù)據(jù)可視化方法，即將任意一個(gè)特征作為x軸，另一個(gè)不同的特征作為y軸)和特征協(xié)方差矩陣分別如圖3和圖4a所示.特征協(xié)方差矩陣C的定義為

(1)

(2)

圖3 原始7維特征空間中的數(shù)據(jù)分布(a1)—(a21) 由任意2個(gè)樣本特征組合所得到的21種數(shù)據(jù)分布，部分樣本特征之間存在較強(qiáng)相關(guān)性，如圖(a13)中的過零率和中心頻率. (a22)—(a28) 7個(gè)樣本特征的頻率直方圖. 所有的特征已做去均值和方差標(biāo)準(zhǔn)化處理.Fig.3 Data distribution in the original 7-dementional feature space(a1)—(a21) show the 21 combinations of data distributions in the subdomains of two features. It is clear that some features present a strong linear correlation, for example the central frequency and zero up-crossing rate shown in (a13). (a22)—(a28) show the distribution of data in the single parameter space of the 7 features. All the features have been normalized by their variances after the means being removed.

圖4 樣本特征的協(xié)方差矩陣(a) 原始7維空間中的樣本特征協(xié)方差矩陣；(b) 新的5維空間中樣本特征協(xié)方差矩陣. 圖(a)的協(xié)方差矩陣存在較大的非對(duì)角元素，表明部分原始特征間存在較強(qiáng)線性相關(guān)性，而在圖(b)的協(xié)方差矩陣中不存在這種現(xiàn)象.Fig.4 Covariance matrixes of sample features(a) Covariance matrix computed from the original 7 seismic features; (b) Covariance matrix computed from the 5 independent new features. Note that some original features shown in (a) present strong linear correlations characterized by significant amplitudes of the off-diagonal elements, which are not seen in (b).

圖5 新的5維特征空間中的數(shù)據(jù)分布(a1)—(a10)表示由任意2個(gè)樣本特征組合得到的10種數(shù)據(jù)分布，樣本數(shù)據(jù)顯示特征之間不存在相關(guān)性. (a11)—(a15)表示5個(gè)新的樣本特征的頻率直方圖，所有特征已做去均值和方差標(biāo)準(zhǔn)化處理.Fig.5 Data distribution in the new 5-dementional space(a1)—(a10) show the 10 combinations of data distributions in the subdomains of two features. There is no obvious correlation between any of the two new features. (a11)—(a15) show the distribution of data in the single parameter space of the new 5 independent features. All the features have been normalized by their variances after the means being removed.

我們針對(duì)相對(duì)獨(dú)立的5個(gè)新特征進(jìn)行K-means算法聚類分析，獲取不同類型的噪聲信號(hào).具體使用的是Python開源庫scikit-learn(Pedregosa et al., 2011)中的K-means算法.在特征空間中，K-means算法首先隨機(jī)選擇已知數(shù)量的樣本點(diǎn)作為初始聚類簇中心；然后計(jì)算剩余樣本至各聚類簇中心的歐氏距離，分配樣本于距離最近的聚類簇；隨后，算法將聚類簇中心的位置更新為相同聚類簇內(nèi)的所有樣本的質(zhì)心位置.迭代執(zhí)行后兩步，直至聚類簇中心的位置不再明顯變化時(shí)，算法結(jié)束迭代.一般而言，如果迭代次數(shù)足夠多，K-means算法總能收斂，但有可能最終收斂于局部極小值，而非全局最優(yōu)解，在很大程度上取決于初始聚類簇中心的選擇.為了解決這個(gè)問題，我們隨機(jī)選擇初始聚類簇中心200次，最終以誤差平方和最小的簇中心作為真正的初始中心.同時(shí)，我們還采用了k-means++方法來進(jìn)行聚類簇初始中心的選取.k-means++算法的基本思路是盡可能選取相互距離遠(yuǎn)的初始聚類中心.具體的步驟為：1) 隨機(jī)選取一個(gè)樣本作為第一個(gè)聚類中心c1；2)計(jì)算每個(gè)樣本與當(dāng)前已有類聚中心最短距離(即與最近一個(gè)聚類中心的距離)，用D(x)表示；這個(gè)值越大，表示被選取作為聚類中心的概率較大，并用輪盤法選出下一個(gè)聚類中心；3)重復(fù)2)，直到選出K個(gè)聚類中心.這種選取初始聚類簇中心的策略雖然簡單，但在一定程度上卻可以極大增加聚類效果的穩(wěn)定性(Arthur and Vassilvitskii, 2007).其次，K-means算法還需要在應(yīng)用前已知最佳的聚類簇?cái)?shù)量(K*).我們通過計(jì)算誤差平方和(SSE)并使用“拐肘法”來獲取K*值.

SSE的計(jì)算公式為

(3)

圖6 原始7維空間旋轉(zhuǎn)至新的5維空間后的誤差平方和(SSE)函數(shù)和平均輪廓系數(shù)曲線(Average Silhouette Score). 當(dāng)K=6時(shí)，誤差平方和出現(xiàn)明顯拐點(diǎn)，如黑點(diǎn)所示，同時(shí)平均輪廓系數(shù)達(dá)到最大值Fig.6 The sum of squared error and average silhouette score computed from the 5-dementional independent features rotated from the original 7-dementional space are plotted as a function of cluster number K. When K=6, the SSE function turns to the point of inflexion, indicated by a solid black circle,meanwhile, the average silhouette score reaches its peak

同時(shí)，我們使用輪廓系數(shù)分析(silhouette coefficient analysis)來進(jìn)一步評(píng)價(jià)K*的可靠性.樣本的輪廓系考慮簇內(nèi)距離和簇間距離對(duì)聚類效果的影響，具體定義為(Rousseeuw, 1987)

(4)

其中，a(i)和b(i)分別代表樣本點(diǎn)i的內(nèi)聚度(簇內(nèi)密集程度)和分離度(簇間分離程度)，a(i)的計(jì)算見公式(5)，其中，j代表與樣本i在同一簇中的不同樣本點(diǎn)，distance(i,j)表示樣本i和j之間的距離.所以，a(i)越小，說明該簇越緊密.b(i)的計(jì)算方式與a(i)類似，只不過計(jì)算的是樣本i與其它簇之間的最小距離.

(5)

可以看出，輪廓系數(shù)的取值范圍為-1～1，當(dāng)取值越趨近于1時(shí)，表示樣本聚類越緊湊，反之，越趨近于-1時(shí)，表示聚類越松散.公式(4)給出了計(jì)算每個(gè)樣本點(diǎn)的輪廓系數(shù)值(silhouette coefficient)方法，在此基礎(chǔ)上，輪廓分?jǐn)?shù)(silhouette score)定義為所有樣本輪廓系數(shù)的平均值.我們計(jì)算了不同K時(shí)的輪廓分?jǐn)?shù)，發(fā)現(xiàn)當(dāng)K=6時(shí)，輪廓分?jǐn)?shù)最大(0.44)(圖6)，對(duì)應(yīng)各樣本點(diǎn)的輪廓系數(shù)值見附圖A1所示.可以看出，在最佳聚類簇?cái)?shù)為6時(shí)，聚類效果最好，且各類間的樣本量比較均衡.

此外，為了驗(yàn)證僅僅使用協(xié)方差矩陣(圖4a)的前5個(gè)本征向量來構(gòu)建旋轉(zhuǎn)矩陣的有效性，我們也測(cè)試了使用7個(gè)全部本征向量的情況.將原始的7維特征空間旋轉(zhuǎn)至7維新特征空間(X′1-X′7),經(jīng)過標(biāo)準(zhǔn)化處理后的新特征散點(diǎn)圖如附圖A2所示，新的7個(gè)特征在統(tǒng)計(jì)上不相關(guān).針對(duì)新特征，同樣使用“拐肘法”來分析K*值(附圖A3)，可以看出，當(dāng)K等于6時(shí)，SSE曲線也出現(xiàn)較為明顯的下降拐點(diǎn)，然而，當(dāng)僅使用前5個(gè)本征向量構(gòu)建的旋轉(zhuǎn)矩陣來處理數(shù)據(jù)時(shí)，得到的拐點(diǎn)位置更為明顯(圖6).為了進(jìn)一步驗(yàn)證K*值取決于5個(gè)獨(dú)立的特征，與非獨(dú)立的特征的數(shù)量無關(guān)，我們通過將過零率和譜峰率兩個(gè)特征簡單線性組合(公式(6))，得到第8個(gè)特征，為了方便我們稱之為綜合特征：

附圖A1 最優(yōu)聚類數(shù)為6時(shí)樣本的輪廓系數(shù)分析不同顏色的橫向條形圖對(duì)應(yīng)不同聚類簇, 每個(gè)條形圖中, 樣本的輪廓系數(shù)值從小到大依次排列. 紅色虛線為所有樣本的平均輪廓系數(shù)值.Appendix Fig. A1 Result of silhouette analysis of the optimal 6 clusters The horizontal bars correspond to different clusters, and the sample are arranged according to its silhouette coefficient, and the red dotted line corresponds to the average silhouette coefficient of all samples.

附圖A2 7維特征空間中的數(shù)據(jù)分布 (a1)—(a21) 由任意2個(gè)樣本特征組合所得到的21種數(shù)據(jù)分布, 7個(gè)新特征之間不存在相關(guān)性. (a22)—(a28) 7個(gè)新樣本特征的頻率直方圖, 所有的樣本特征已做去均值和方差標(biāo)準(zhǔn)化處理.Appendix Fig. A2 Data distribution in the 7-dementional feature space (a1)—(a21) show the 21 combinations of data distributions in the subdomains of two features. These subplots suggest that there is almost no correlation between any of the features. (a22)—(a28) show the distribution of data in the single parameter space of the 7 new features. All the features have been normalized by their variances after the means being removed.

附圖A3 新的7維特征空間中誤差平方和函數(shù)(SSE)與聚類簇?cái)?shù)之間的關(guān)系. 當(dāng)K等于6時(shí), SSE曲線出現(xiàn)明顯拐點(diǎn), 如紅點(diǎn)所示Appendix Fig. A3 Sum of squared error (SSE) as a function of number clusters computed in the new 7-D feature space. The SSE curve starts to bend at K=6, which is shown as solid red point

附圖A4 假設(shè)的8維特征空間中的數(shù)據(jù)分布 (a1)—(a28)由任意2個(gè)樣本特征組合所得到的28種數(shù)據(jù)分布, 部分特征之間存在較強(qiáng)相關(guān)性, 如圖(a28)的譜峰率和綜合特征之間的相關(guān)性. (a29)—(a36)表示8個(gè)樣本特征(原始7維+假設(shè)的綜合特征)的頻率直方圖. 所有特征已做去均值和方差標(biāo)準(zhǔn)化處理.Appendix Fig. A4 Data distribution in the hypothetic 8-dementional feature space (a1)—(a28) show the 28 combinations of data distributions in the subdomains of two features. It is clear that there are correlations among some of the features, for example, the spectral peak frequency and the comprehensive features shown in (a28). (a29)—(a36) show the distribution of data in the single parameter space of the 8-dementional feature (original 7 features plus the hypothetic comprehensive feature). All the features have been normalized by their variances after the means being removed.

附圖A5 樣本特征的協(xié)方差矩陣 (a) 8維特征空間(原始7維+假設(shè)的綜合特征)的樣本特征協(xié)方差矩陣; (b) 由8維特征空間協(xié)方差矩陣的5個(gè)主要本特向量為特征的5維空間樣本特征協(xié)方差矩陣. 圖(a)的協(xié)方差矩陣中存在較大的非對(duì)角元素, 表明部分特征之間存在較強(qiáng)線性相關(guān)性, 而在圖(b)的協(xié)方差矩陣中不存在這種現(xiàn)象.Appendix Fig. A5 Covariance matrixes of sample features (a) Covariance matrix computed from the 8-dimenional feature space (original 7 features plus the hypothetic comprehensive feature); (b) Covariance matrix computed from the 5 independent new features derived from the major eigenvectors. Note that some features shown in (a) present strong linear correlations characterized by significant amplitudes of the off-diagonal elements, which are not seen in (b).

附圖A6 由8維假想特征空間計(jì)算得到的新5維空間中的數(shù)據(jù)分布 (a1)—(a10)由任意2個(gè)樣本特征組合得到的10種數(shù)據(jù)分布, 樣本數(shù)據(jù)顯示特征之間不存在相關(guān)性. (a11)—(a15)表示5個(gè)新的樣本特征的頻率直方圖, 所有特征已做去均值和方差標(biāo)準(zhǔn)化處理. 注意該5維獨(dú)立樣本特征是由附圖A4與附圖A5的8維假想特征空間計(jì)算獲得的，而圖5中的5維新樣本特征由原始的7維特征空間得到.Appendix Fig. A6 Data distribution in the new 5-dementional space (a1)—(a10) show the 10 combinations of data distributions in the subdomains of two features. There is no obvious correlation between any of the two new features. (a11)—(a15) show the distribution of data in the single parameter space of the new 5 independent features. All the features have been normalized by their variances after the means being removed. Note that the 5 independent features are derived from the hypothetic 8-dimensional feature space shown in Appendix Fig. A4 and A5, where the 5 independent features shown Fig.5 are derived from the original 7-dimensional feature space.

綜合特征=0.2×過零率+0.8×譜峰率.

(6)

附圖A4為經(jīng)過標(biāo)準(zhǔn)化處理后的8個(gè)特征的散點(diǎn)圖，其中，綜合特征與譜峰率之間存在非常明顯的線性相關(guān)性.通過計(jì)算8個(gè)特征的協(xié)方差矩陣(附圖A5a)的本征值(λ1=4.81,λ2=1.65,λ3=1.11,λ4=0.25,λ5=0.13,λ6=0.04,λ7=0.02,λ8=-6.87×10-16)，可以看出，前5個(gè)特征占有樣本信息量的99.25%，因此，同樣僅選取前5個(gè)主要本征值對(duì)應(yīng)的本征向量來構(gòu)建樣本獨(dú)立特征空間，將原始的8維特征空間旋轉(zhuǎn)至5維新特征空間(X′1-X′5).經(jīng)過標(biāo)準(zhǔn)化處理后的5個(gè)新特征的散點(diǎn)圖和協(xié)方差矩陣分別如附圖A6和附圖A5b所示，進(jìn)一步驗(yàn)證K*=6的結(jié)論(附圖A7).因此，PCA方法能夠通過消除樣本冗余信息來降低計(jì)算復(fù)雜度，同時(shí)不影響聚類分析的結(jié)果.

3 結(jié)果

我們采用批量(mini batch)訓(xùn)練的策略，即分批次向模型輸入訓(xùn)練樣本，直至損失函數(shù)收斂.這種優(yōu)化想法非常樸素：當(dāng)總樣本量非常大時(shí)，會(huì)導(dǎo)致模型訓(xùn)練時(shí)間過長，此時(shí)，可以不斷從整體中隨機(jī)抽取一定數(shù)量的樣本來訓(xùn)練模型.具體使用的是scikit-learn中的MiniBatchKMeans算法.其中，隨機(jī)初始化質(zhì)心200次，聚類簇?cái)?shù)取6，最大迭代次數(shù)設(shè)置為1000，每次的batch大小設(shè)置為100，其余參數(shù)選用默認(rèn)值.訓(xùn)練所用的服務(wù)器CPU為Intel Xeon E5-2687W，內(nèi)存125G.當(dāng)訓(xùn)練第1069個(gè)mini batch時(shí)，損失函數(shù)收斂，此時(shí)的慣性量(mean batch inertia)1.29，總共訓(xùn)練用時(shí)78 s.最終，K-means模型將訓(xùn)練噪聲樣本聚類為6類，各類樣本數(shù)量的比例依次為：Label 0(16.23%)、1(19.20%)、2(6.85%)、3(32.69%)、4(14.85%)和5(10.18%).我們進(jìn)而利用已訓(xùn)練好的K-means模型，通過在特征空間中搜尋與待預(yù)測(cè)樣本歐式距離最近的聚類簇中心的方式，實(shí)現(xiàn)對(duì)測(cè)試樣本(5,085,821)的標(biāo)簽進(jìn)行預(yù)測(cè).

值得說明的是，我們雖然使用PCA方法對(duì)原始的7維特征降維至5維空間，然后再在降維后的新特征空間中使用K-means方法對(duì)樣本進(jìn)行聚類分析，但在5維空間中的分析仍有可能存在維度過高的問題.為了進(jìn)一步說明我們聚類分析的可靠性，我們使用t-SNE方法對(duì)樣本聚類結(jié)果進(jìn)行可視化分析.t-SNE將原始空間和嵌入空間聯(lián)合概率的Kullback-Leibler(KL)散度定義為損失函數(shù)，然后通過對(duì)損失函數(shù)進(jìn)行優(yōu)化來達(dá)到降維的目的.值得注意的是，此損失函數(shù)為非凸函數(shù)，所以降維結(jié)果有可能收斂于局部最優(yōu)解.我們采用scikit-learn中的t-SNE模塊進(jìn)行分析，其中，perplexity參數(shù)為70，其余均為默認(rèn)參數(shù).同時(shí)，考慮到t-SNE算法計(jì)算極慢，我們分別從聚類后的6類噪聲中隨機(jī)抽取10000個(gè)樣本，利用總共60000個(gè)樣本進(jìn)行降維可視化分析(圖7).可以看出，6類噪聲能夠較為清晰地聚類.我們也注意到Label 5聚類簇被分割為幾個(gè)明顯的小聚類簇，這可能是由于t-SNE算法本身優(yōu)化的是非凸損失函數(shù)，圖7的降維結(jié)果只是一個(gè)局部最優(yōu)解，進(jìn)而造成某一簇被分割為幾個(gè)不同位置的小聚類簇(van der Maaten, 2018).

圖7 基于t-SNE降維的數(shù)據(jù)可視化. 分布顯示不同的簇類可以被較好地區(qū)分Fig.7 Data visualization using the t-SNE analysis,suggesting various clusters are well separated

各類噪聲的振幅譜如圖8所示，可以看出，6類噪聲的振幅譜特征存在明顯的差異.其中，Label 0噪聲信號(hào)主頻約為2.5 Hz，且振幅較大(圖8a1)；而Label 2噪聲的主頻約為11 Hz(圖8a3)；Label 4噪聲在較低頻帶(2～11 Hz)具有較寬的頻帶范圍，且振幅較大(圖8a5)；Label 1信號(hào)振幅較小，且從12～100 Hz，振幅整體呈現(xiàn)減小的趨勢(shì)(圖8a2)；而Label 5噪聲在1～100 Hz頻帶范圍內(nèi)有三個(gè)較為明顯的主峰，其中，在50 Hz頻率處存在一個(gè)非常明顯的主峰(圖8a6)；Label 3噪聲則只在2～20 Hz頻率范圍內(nèi)存在兩個(gè)明顯的主峰(圖8a4).噪聲信號(hào)在時(shí)間域的波形如圖9所示，各類波形非常相似，難以區(qū)分.相比而言，K-means聚類模型能較好地分離和標(biāo)注不同噪聲信號(hào).

圖8 聚類后的6類噪聲信號(hào)的疊加振幅譜其中黑色實(shí)線表示中位值，兩個(gè)灰色點(diǎn)劃線分別表示25%和75%分位數(shù). 鑒于連續(xù)波形采樣間隔為0.005 s，并且采用1 Hz的高通濾波，所以僅展示1～100 Hz的頻帶范圍.Fig.8 Stacked spectra of six-class noise signalsThe median is shown with the solid black line, and the gray dashed lines correspond to the 25th and 75th percentiles, respectively. Since the waveform is sampled in 0.005 second, and the pre-processing employs a 1-Hz high-pass filter, we only plot the spectra in the frequency band of 1～100 Hz.

4 討論

4.1 標(biāo)簽物理意義

針對(duì)無監(jiān)督聚類問題，一般可以通過分析聚類樣本的特征和實(shí)際觀測(cè)環(huán)境來推斷標(biāo)簽可能的物理意義.我們通過分析不同位置的臺(tái)站記錄的各個(gè)標(biāo)簽出現(xiàn)的累計(jì)次數(shù)和樣本時(shí)間-頻率域的特征，對(duì)各類標(biāo)簽可能代表的噪聲源進(jìn)行推斷.

圖10是臺(tái)站SP54所記錄樣本的標(biāo)簽統(tǒng)計(jì)(圖10a)和臺(tái)站實(shí)際的觀測(cè)背景(圖10b).臺(tái)站SP54被安置于建筑物和高速公路附近，所記錄的波形以Label 0和Label 4噪聲為主.類似地，臺(tái)站SP51和SP53也被布設(shè)于建筑物和高速公路附近，所記錄的噪聲類型統(tǒng)計(jì)和臺(tái)站SP54也非常相似(附圖A8).同時(shí)，考慮到Label 0類噪聲信號(hào)主頻較低(2.5 Hz)，且振幅較大的特點(diǎn)(圖8)，我們推斷Label 0噪聲信號(hào)可能來自于風(fēng)擾動(dòng)低層建筑物引起的地面震動(dòng) (Johnson et al., 2019).不同于臺(tái)站SP54、SP51和SP53距離高速公路的最短距離大于200 m的布設(shè)背景，臺(tái)站SP58與附近高速公路之間的最短距離僅約36 m，同時(shí)臺(tái)站SP58所記錄的噪聲類型以Label 4噪聲(頻帶較寬，振幅較大)為主(附圖A8)，因此，可以我們推斷Label 4噪聲代表軌道交通活動(dòng)所引起的地面震動(dòng)信號(hào) (Johnson et al., 2019; Meng et al., 2019)，類似的規(guī)律也可參考臺(tái)站SP69、SP3E、SP68、SP3C和SPBC(附圖A8—附圖A9).臺(tái)站SPB2、SPA3、SPD1和SPC2則以Label 1噪聲為主(附圖A10)，同時(shí)Label 1噪聲的振幅較小(圖8)，結(jié)合臺(tái)站SPC2實(shí)際的觀測(cè)環(huán)境，我們推斷Label 1噪聲來源于風(fēng)擾動(dòng)樹木造成的地表震動(dòng)(Johnson et al., 2020).臺(tái)站SPA3和SPC2也包含Label 5噪聲(附圖A10)，且Label 5噪聲的振幅譜在50 Hz處存在非常明顯的信號(hào)(圖8)，所以，我們推斷Label 5噪聲來源于附近電力輸送線的電流電離空氣所引起的地面震動(dòng)(附圖A10).

圖9 臺(tái)站SPE3垂向分量記錄的部分波形(a) 2 h長度的連續(xù)波形記錄. 參考時(shí)刻為北京時(shí)間2017-09-12T00∶00∶00.000. 不同的顏色代表被聚類為不同的噪聲類型；(b)—(g) 展示了長度為2 s的連續(xù)波形，分別對(duì)應(yīng)6種不同的噪聲信號(hào). 6類噪聲信號(hào)在時(shí)間域的波形極為相似，難以區(qū)分.Fig.9 Noise waveforms recorded by the vertical component of the station SPE3(a) Example of a two-hour continuous record (the reference time here is 2017-09-12T00∶00∶00.000, Beijing time) and its composition of the 6 types noise signals shown in different colors. (b)-(g) correspond to the waveforms of different noise signal with a length of 2 s. These six-class noise signals are indistinguishable in time domain.

圖10 (a) 臺(tái)站SP54的連續(xù)波形中各類噪聲信號(hào)的累計(jì)次數(shù). 觀測(cè)時(shí)間內(nèi)，此臺(tái)站主要記錄到Label 0和Label 4噪聲信號(hào)； (b) 此臺(tái)站的布設(shè)環(huán)境. 圖像來源于臺(tái)站運(yùn)行時(shí)間內(nèi)的谷歌地圖歷史影像Fig.10 (a) Normalized cumulative number of six labels at station SP54. During its operation, this station is dominated by the Label 0 and Label 4 noise. (b) The local environment of this station. This image comes from the historical image of the Google Maps in the operation of the station

值得說明的是，這里有關(guān)不同類型噪聲聚類分析依然比較初步，比如識(shí)別的噪聲信號(hào)可能是幾種其它噪聲信號(hào)的疊加(圖8和圖9)，在未來的研究中，我們將嘗試使用更多的聚類簇?cái)?shù)目(K*)對(duì)噪聲信號(hào)做更詳細(xì)的剖析.其次，準(zhǔn)確可靠地識(shí)別和解釋各類噪聲信號(hào)，需要具備非常好的觀測(cè)前提，如Johnson等(2019)對(duì)與風(fēng)有關(guān)的噪聲信號(hào)的識(shí)別就得益于密集臺(tái)陣和風(fēng)速計(jì)的布設(shè)，以及對(duì)臺(tái)站布設(shè)環(huán)境的詳細(xì)記錄.最后，利用本文已訓(xùn)練的K-means聚類模型可以幫助我們?cè)谖磥淼难芯恐锌焖贅?biāo)注不同類型的噪聲信號(hào)，進(jìn)而在監(jiān)督學(xué)習(xí)中用來識(shí)別某些特定信號(hào)，如軌道交通信號(hào)等 (Johnson et al., 2019; Meng and Ben-Zion, 2018a; 劉磊和蔣一然, 2019).

4.2 噪聲信號(hào)的時(shí)空演化

文獻(xiàn)調(diào)研發(fā)現(xiàn)，關(guān)于二維密集臺(tái)陣噪聲信號(hào)的時(shí)間-空間演化規(guī)律的研究較少，尤其針對(duì)橫跨城市-農(nóng)村混合區(qū)域觀測(cè)數(shù)據(jù)的研究.圖11a和11b分別為2017年9月12日凌晨2點(diǎn)和白天14點(diǎn)時(shí)整個(gè)臺(tái)陣的噪聲分布.當(dāng)夜晚安靜時(shí)，臺(tái)陣北部區(qū)域以風(fēng)擾動(dòng)建筑物所產(chǎn)生的噪聲信號(hào)(Label 0)為主，而在白天喧鬧時(shí)刻，北部區(qū)域則以軌道交通引起的噪聲信號(hào)為主(Label 4).對(duì)于以山區(qū)為主，人口分布較為稀疏的南部區(qū)域，無論夜晚還是白天，都以風(fēng)擾動(dòng)樹木產(chǎn)生的噪聲信號(hào)為主(Label 1).類似的規(guī)律也可以從2017年9月5日、7日和12日連續(xù)24 h的噪聲類型動(dòng)態(tài)變化中看出.

為了進(jìn)一步定量分析不同臺(tái)站布設(shè)位置處的噪聲類型，我們計(jì)算了每個(gè)臺(tái)站所記錄到的各個(gè)噪聲類型所占的比重.同時(shí)，附圖A11a展示了Label 0噪聲信號(hào)在二維空間上的 分布規(guī)律，可以看出，Label 0噪聲信號(hào)主要出現(xiàn)在臺(tái)陣的北部區(qū)域，對(duì)應(yīng)于風(fēng)擾動(dòng)建筑物所產(chǎn)生的噪聲信號(hào).而由風(fēng)擾動(dòng)樹木產(chǎn)生的噪聲信號(hào)(Label 1)則主要分布于臺(tái)陣的南部區(qū)域(附圖A11b).其次，與軌道交通相關(guān)的噪聲信號(hào)(Label 4)主要出現(xiàn)于城市區(qū)域(附圖A11c).與電力輸送線有關(guān)的噪聲信號(hào)(Label 5)在各個(gè)臺(tái)站上記錄到的比重相對(duì)較小(附圖A11d).因此，噪聲信號(hào)的時(shí)間-空間演化規(guī)律與臺(tái)站布設(shè)區(qū)域的城市化建設(shè)程度密切相關(guān).由于我們使用的密集臺(tái)陣的臺(tái)站間距較大(平均5 km)，所以噪聲信號(hào)在空間上的相關(guān)性不明顯.未來在城市監(jiān)測(cè)時(shí)，如果觀測(cè)臺(tái)陣足夠密集，可以借鑒本文的分析思路，實(shí)現(xiàn)空間上可能相關(guān)的公路交通信號(hào)或者高鐵地震信號(hào)的觀測(cè)，進(jìn)而可能利用軌道交通信號(hào)的時(shí)空變化特征來對(duì)交通工具的運(yùn)行軌跡進(jìn)行監(jiān)測(cè) (Snover et al., 2020; 劉磊和蔣一然, 2019).

圖11 不同噪聲信號(hào)在臺(tái)陣上的空間分布臺(tái)站的顏色對(duì)應(yīng)于1 h范圍內(nèi)所記錄次數(shù)最多的噪聲類型. (a) 北京時(shí)間2017-09-12T02∶00∶00.000時(shí)刻的噪聲空間分布；(b) 北京時(shí)間2017-09-12T14∶00∶00.000時(shí)刻的噪聲空間分布. 北部區(qū)域夜晚主要的噪聲來源于風(fēng)吹動(dòng)建筑物產(chǎn)生的地面振動(dòng)，白天則主要來源于軌道交通；而南部鄉(xiāng)村地區(qū)與山區(qū)，無論白天還是夜晚，噪聲都主要來源于風(fēng)吹動(dòng)樹木產(chǎn)生的地面振動(dòng)信號(hào).Fig.11 Maps showing distribution of different types of noise across the dense arrayStation colors correspond to the noise labels occurring most frequently during one-hour period. (a) Local time from 2017-09-12T02∶00∶00.000; (b) Local time from 2017-09-12T14∶00∶00.000. In the northern part of the array, the nighttime recordings are dominated by noise from interactions between the wind and local buildings, and the daytime noises are predominantly induced by traffic. Meanwhile, stations in the southern part of the array are prevailed by ground shaking related to wind-driven tree motions in both daytime and nighttime.

我們利用的是傳統(tǒng)聚類方法中的K-means算法，此類傳統(tǒng)聚類方法在人工選取的特征空間中對(duì)樣本進(jìn)行聚類.需要指出的是，當(dāng)人工選取的特征不合適時(shí)，就有可能出現(xiàn)較差的聚類效果(Zhou et al., 2019).在本研究中，當(dāng)K=6時(shí)，采用K-means聚類得到的平均輪廓系數(shù)值只有0.44(圖6)，且不同類之間也存在少量樣本的混疊(圖7).產(chǎn)生這些問題的原因很有可能是與信號(hào)特征提取不夠充分有關(guān)，因此我們K-means分析得到的可能是一個(gè)比較初步的噪聲聚類結(jié)果.近些年快速發(fā)展起來的深度聚類算法，借助深度神經(jīng)網(wǎng)絡(luò)的強(qiáng)大非線性映射能力，采用深度神經(jīng)網(wǎng)絡(luò)替代人工來設(shè)計(jì)來提取更合適的樣本特征、提升聚類算法的性能，同時(shí)，聚類結(jié)果的反饋又能幫助神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)到更好的特征，進(jìn)而改善最終的樣本聚類效果.這類深度聚類算法已被應(yīng)用地震學(xué)領(lǐng)域，如Snover等(2020)使用深度嵌入聚類模型從噪聲信號(hào)的頻譜圖像中識(shí)別出了3種和軌道交通震源相關(guān)的噪聲信號(hào)，Jenkins等(2021)使用高斯混合模型和深度嵌入聚類兩種深度聚類算法，識(shí)別出了88種地震信號(hào)，并且分析發(fā)現(xiàn)地震時(shí)空演化特征與季節(jié)性和潮汐活動(dòng)之間存在相關(guān)性，Seydoux等(2020)則分別使用深散射網(wǎng)絡(luò)和高斯混合模型來進(jìn)行特征的提取和地震信號(hào)的識(shí)別工作.在深度聚類中，不同的數(shù)據(jù)類型常常被用來作為模型的輸入(Mousavi et al., 2016)，在后續(xù)研究中，我們將嘗試多種類型數(shù)據(jù)的組合來作為模型輸入，如時(shí)頻譜圖像和多道波場灰度圖像，其次，針對(duì)現(xiàn)有的數(shù)據(jù)，測(cè)試選擇合適的深度聚類模型，進(jìn)一步聚類分析本論文的噪聲類型.

5 結(jié)論

我們基于無監(jiān)督機(jī)器學(xué)習(xí)算法對(duì)布設(shè)于郯廬斷裂帶濰坊段的密集觀測(cè)臺(tái)陣的噪聲信號(hào)進(jìn)行聚類分析，得到以下結(jié)論.

(1)雖然不同噪聲信號(hào)在時(shí)間域的波形非常相似，但是基于K-means模型可以較好地聚類出6類不同噪聲信號(hào).

(2)綜合分析噪聲信號(hào)的波形、振幅譜和觀測(cè)臺(tái)站的布設(shè)環(huán)境，我們對(duì)部分噪聲信號(hào)來源進(jìn)行了初步解釋，認(rèn)為可能的來源包括軌道交通、電力輸送線和風(fēng)擾動(dòng)建筑物、樹木導(dǎo)致的地面震動(dòng).未來還需要分析噪聲信號(hào)與物理震源的時(shí)空相關(guān)性，這需要進(jìn)一步收集噪聲源的產(chǎn)生的時(shí)間、地點(diǎn)和持續(xù)時(shí)間等信息.

(3)通過分析噪聲信號(hào)的時(shí)空演化規(guī)律，我們發(fā)現(xiàn)噪聲信號(hào)的時(shí)空變化與臺(tái)陣布設(shè)區(qū)域內(nèi)的人口密集程度或城市化發(fā)展程度密切相關(guān).

致謝感謝編輯、編委與二位評(píng)審專家對(duì)完善本文提出的具有建設(shè)性的意見和建議.從國家地震數(shù)據(jù)科學(xué)中心獲取了臺(tái)陣運(yùn)行期間的地震目錄.本文使用的機(jī)器學(xué)習(xí)相關(guān)算法來自于Python開源程序庫scikit-learn(Pedregosa et al., 2011)，地圖來自于Bigemap軟件，其余圖件則基于Python程序庫matplotlib(Hunter, 2007)繪制.數(shù)據(jù)處理過程中，中國石油大學(xué)(北京)的趙楊教授、陳海潮副研究員和張樂樂同學(xué)提出了很多寶貴意見，在此一并表示感謝.

附圖A7 由假想的8維空間旋轉(zhuǎn)至新的5維空間后的誤差平方和(SSE)曲線與聚類簇?cái)?shù)之間的關(guān)系. 當(dāng)K等于6時(shí), 誤差平方和曲線出現(xiàn)明顯拐點(diǎn), 如紅點(diǎn)所示Appendix Fig.A7 The sum of squared error (SSE) computed from the 5-dementional independent features rotated from the hypothetic 8-dementional space is plotted as a function of cluster number K. When K=6, the SSE function turns to the point of inflexion, indicated by a solid red circle

附圖A8 臺(tái)站記錄連續(xù)波形中各類噪聲信號(hào)累計(jì)次數(shù)和臺(tái)站實(shí)際布設(shè)環(huán)境圖像來源于臺(tái)站運(yùn)行期間的谷歌地圖歷史影像.Appendix Fig.A8 The normalized cumulative number of six labels and local environment for stationThe image of local environment comes from the historical image of Google Maps in the operation time of stations.

附圖A9 臺(tái)站記錄連續(xù)波形中各類噪聲信號(hào)累計(jì)次數(shù)和臺(tái)站實(shí)際布設(shè)環(huán)境圖像來源于臺(tái)站運(yùn)行期間內(nèi)的谷歌地圖歷史影像.Appendix Fig.A9 The normalized cumulative number of six labels and local environment for stationThe image of local environment comes from the historical image of Google Maps in the operation time of stations.

附圖A10 臺(tái)站記錄連續(xù)波形中各類噪聲信號(hào)累計(jì)次數(shù)和臺(tái)站實(shí)際布設(shè)環(huán)境圖像來源于臺(tái)站運(yùn)行期間內(nèi)的谷歌地圖歷史影像.Appendix Fig.A10 The normalized cumulative number of six labels and local environment for stationThe image of local environment comes from the historical image of Google Maps in the operation time of stations.

附圖A11 噪聲信號(hào)在每個(gè)臺(tái)站中出現(xiàn)的累計(jì)比重 (a) Label 0噪聲; (b) Label 1噪聲; (c) Label 4噪聲; (d) Label 5噪聲.Appendix Fig.A11 Cumulative percentage of noises at each station for (a) noise Label 0; (b) noise Label 1; (c) noise Label 4; (d) noise Label 5