?

結合LSTM 和Self-Attention 的滾動軸承剩余使用壽命預測方法

2024-01-10 01:40高俊峰李周正江志農高金吉
振動工程學報 2023年6期
關鍵詞:皮爾遜頻段軸承

黃 宇,馮 坤,高俊峰,李周正,江志農,高金吉

(1.北京化工大學信息科學與技術學院,北京 100029;2.北京化工大學高端機械裝備健康監控與自愈化北京市重點實驗室,北京 100029;3.中國石油天然氣股份有限公司煉油與化工分公司,北京 100007)

引言

滾動軸承是機械設備中最常見且極其重要的關鍵零部件。在設備運行的過程中由于環境、工況等因素的影響,滾動軸承因易出現點蝕、磨損等現象,而成為工業上最易損壞、可靠性最差的零部件之一。在旋轉設備中,約有30%的故障是由滾動軸承引起的[1]。因此,對滾動軸承的運行狀態進行監測,預測其剩余使用壽命并確定最佳維修點,可以為決策者建立維修方案提供支撐,具有重要意義。

對滾動軸承的剩余使用壽命預測主要有三種主流的方法:基于模型的方法、數據驅動方法和二者混合的方法[2]?;谀P偷姆椒ㄐ枰罅康膶<抑R和先驗知識,不僅浪費時間和勞動力,而且不具有通用性。因此,研究數據驅動的剩余使用壽命預測方法,實現軸承精準的RUL 預測,是目前的研究熱點?;跀祿寗拥妮S承RUL 預測主要有三個步驟:①從加速度振動信號中提取具有單調性和趨勢性的退化特征,表征滾動軸承的衰退過程;②構建健康度指標模型,利用傳統人工智能或深度學習的方法學習退化特征與軸承健康度指標之間的復雜映射關系,得到軸承的退化曲線;③對軸承退化曲線進行擬合預測,計算軸承壽命失效點,最終得到軸承剩余使用壽命。

退化特征提取是軸承RUL 預測的前提條件,目的是從軸承原始加速度振動信號中提取出能反映其退化趨勢的特征。目前常見的方法是利用深度學習方法從時域、頻域以及時頻域中提取原始振動信號的退化特征。楊宇等[3]提出一種改進的深度信念網絡,直接以滾動軸承原始振動信號作為網絡輸入,挖掘出原始振動信號深層本質特征。李少鵬[4]將原始信號通過快速傅里葉變換轉換為頻域幅值信號,然后利用一維卷積神經網絡從頻域信號中挖掘深層退化特征。WANG 等[5]提出了通過將原始信號的時頻表示作為輸入,并利用三維深度卷積神經網絡提取退化特征的方法。CAO 等[6]將原始信號邊緣譜作為輸入,利用時間卷積網絡提取高級退化特征表示。

上述基于深度學習的特征提取方法往往需要大量標簽對深度學習網絡進行監督學習,標簽的缺乏嚴重制約了利用深度學習進行退化特征提取,而傳統的統計特征難以反映軸承的整個退化過程,單調性和趨勢性不明顯[7-9]。

此外,利用退化特征構建健康度指標是軸承RUL 預測至關重要的一步,目的是將退化特征映射為趨勢性健康度指標。如果模型輸出的健康度指標能夠很好地反映軸承的退化趨勢,便能準確地預測軸承的剩余使用壽命。軸承的振動數據是一種時間序列數據,而循環神經網絡(Recurrent Neural Network,簡稱RNN)具有強大的時間序列處理能力,因此,RNN 被廣泛應用于滾動軸承剩余使用壽命預測。GUO 等[10]利用RNN 將提取的8 個時域特征和6 個頻域特征映射為健康度指標,并用雙指數函數模型擬合預測,得到了很高的預測精度。韓林潔[11]構建了頻域幅值累計特征,并用門控循環單元網絡進行軸承剩余使用壽命預測。LI 等[12]利用核主成分分析和指數加權移動平均獲得了退化曲線,并利用分層門控循環單元網絡進行剩余使用壽命預測。WU 等[13]提出了一種動態差異技術,從原始監測數據中提取特征,利用LSTM 學習監測數據背后真正的物理退化機制,并通過實驗說明了LSTM 的性能比標準RNN 和門控循環單元(GRU)更優越。

LSTM 是RNN 的變體,其性能比RNN 更為優越。然而,LSTM 仍然有一些不可忽視的缺點。首先,LSTM 不能進行并行處理,LSTM 必須等待前一個數據處理完成才能處理下一個數據,這不僅降低了模型的靈活度,還導致誤差逐級累積[14]。其次,LSTM 在訓練過程中存在梯度爆炸、梯度消失以及占用大量內存等問題,而這一系列問題目前沒有很好的解決方法。此外,當輸入序列超過一定長度后,LSTM 會出現記憶衰退現象,即不能有效記憶很久之前的信息[15]。而自注意力機制(Self-Attention)可以很好地解決上述問題。Self-Attention 能讓模型對關鍵信息重點關注并充分學習吸收,能學習任意長度的時間序列信息;且Self-Attention 是并行計算方式,處理速度較LSTM 大大加快[16]。

綜上,提出一種利用包絡譜特征并結合LSTM和Self-Attention 構建趨勢性健康度指標來實現軸承RUL 預測的方法。LSTM 將退化特征編碼為高階特征向量,Self-Attention 對LSTM 隱藏層的輸出計算權重系數,并且選擇性地保留中間結果,使隱藏層輸出向量之間的聯系更加緊密,強調關鍵信息和減少對無效信息的關注,實現了滾動軸承RUL 的準確預測。

1 理論背景

1.1 包絡譜特征和皮爾遜相關系數

1.1.1 包絡譜特征

包絡解調可以有效地將軸承故障信號從高頻信號中分離并提取出來[17]。軸承在出現故障時,由軸承故障引起的特征信號被調制到高頻帶段,此時時域波形和頻譜均難以明顯體現其故障特征。包絡解調方法對信號進行希爾伯特變換和快速傅里葉變換后得到原始信號的包絡譜,從包絡譜上可以獲取到清晰的故障特征。同時,包絡解調還可以發現軸承更早期的缺陷,從而提取軸承初期故障特征[18]。

軸承在退化過程中,多數情況下不止存在一種故障即單一故障模式,而是存在多種故障即復合故障模式。因此,各種故障頻率會在高頻段與多部件的固有頻率發生共振。對信號包絡解調后,可以在低頻段得到故障特征頻率。傳統的包絡解調在得到包絡譜后利用低通濾波器濾掉高頻段信號,消除高頻段信號對故障分析的干擾,但鑒于復合故障模式的復雜性,包絡譜的高頻段仍可能包含故障特征[19-20]。因此,本文嘗試將包絡譜按頻率均分成n段。這使得在分析低頻段故障特征時不會引入高頻段信號的干擾,同時還保留了高頻段中的故障信息;由于高頻段故障特征相對較微弱,按頻率分段后再分析高頻段信號,減弱了低頻段故障信號對機器學習模型特征提取時的絕對主導性,一定程度上放大了高頻段信號的故障特征。

1.1.2 皮爾遜相關系數

皮爾遜相關系數用于度量兩個向量X和Y之間的相關性,定義如下式所示:

式中Cov表示協方差;D表示方差;ρXY表示向量X和向量Y之間的皮爾遜相關系數,取值范圍為[-1,1],|ρXY|越大,說明X和Y的相關性越大。

皮爾遜相關系數消除了不同變量量綱上的差別,即兩個變量的位置和尺度變化并不會改變其皮爾遜相關系數,所以通過皮爾遜相關系數理論計算得到的相關性所衡量的是趨勢[21]。因此,將包絡譜分段后各子頻段與標準樣本的皮爾遜相關系數作為退化特征,更有利于模型準確表征軸承退化過程的健康度指標。

1.2 長短期記憶網絡

LSTM 網絡是從標準RNN 改進得來的。LSTM 通過其內部復雜的門運算和引入細胞態,有效緩解標準RNN 的長期依賴問題[22]。遺忘門ft決定上一時刻的細胞態ct-1有多少信息保留到當前時刻的細胞態ct;輸入門it決定當前時刻的輸入有多少信息存儲到當前時刻的細胞態ct中;輸出門ot控制當前細胞態ct有多少信息保留到當前時刻輸出ht中。圖1 為LSTM 單元的內部結構,xt-1,xt和xt+1分別指t-1,t和t+1 時刻的 輸入信 息;ht-1,ht和ht+1分別指t-1,t和t+1 時刻的輸出信息。

圖1 LSTM 結構Fig.1 The structure of LSTM

遺忘 門ft、輸入門it、輸出門ot、細 胞態ct和輸出ht的計算公式分別如下式所示:

式中xt為當前時刻的輸入;ht-1為上一時刻的輸出;W為權重矩陣;b為偏置;σ(x)=1/(1+ex)為Sigmoid 激活函 數;下 標“f”,“i”,“c”分別表 示矩陣W和偏置b為分別代表遺忘門、輸入門、細胞態的參數;?表示逐元素相乘法。

1.3 自注意力機制

Self-Attention 是注意力機制的改進,不僅可以快速篩選出關鍵信息,減少對其他無關信息的關注,還可以降低對外部信息的依賴,更易捕捉輸入數據的內部相關性[23]。神經網絡通過引入自注意力機制,在解決模型信息過載問題的同時,還提高了網絡的準確率和魯棒性[24]。

Self-Attention 的計算分為兩步。步驟1:計算輸入序列任意向量之間的注意力權重;步驟2:根據注意力權重計算輸入序列的加權平均值。自注意力機制如圖2 所示,a(ii=1,2,3,…,t)表示輸入序列;v(ii=1,2,3,…,t)表示由輸入序列生成的值向量;αti(i=1,2,3,…,t)表示輸入序列與各自的向量q和k做運算并經過Softmax 函數后 的結果;b(ii=1,2,3,…,t)則表示輸入序列中第i個位置信息與所有位置信息進行注意力機制運算后的結果。

圖2 自注意力機制結構Fig.2 The structure of Self-Attention mechanism

Self-Attention 具體運算如下式所示:

式中Q,K和V分別為查詢矩陣、鍵矩陣和值矩陣,由輸入X分別與相應的權重矩陣Wq,Wk和Wv相乘得到;dim表示Q,K和V的維數。

2 RUL 預測方法及流程

針對基于深度學習的退化特征提取方法的不足,提出了一種基于包絡譜特征和皮爾遜相關系數的退化特征提取方法;并結合 LSTM 和Self-Attention 各自的優點,利用退化特征構建健康度指標,準確表征軸承的退化過程,實現滾動軸承RUL 的預測,流程框圖如圖3 所示。

圖3 剩余使用壽命預測方法流程圖Fig.3 The flowchart of RUL prediction method

具體流程如下:

步驟1:計算信號的包絡譜,將包絡譜按頻率平均劃分為n段。以軸承正常運行的振動信號作為標準樣本,分別計算各樣本每個子頻段及標準樣本對應頻段的皮爾遜相關系數,將計算得到的相關系數作為軸承退化特征。

步驟2:模型的訓練數據集是軸承全壽命周期振動信號所提取的退化特征。同時,由于軸承的退化是一個漸變的過程,因此設時刻為t的訓練樣本的健康度指標HI為:

式中T為軸承的全壽命時間;t為軸承當前運行的時間。將式(11)的計算結果作為訓練集的標簽,標簽的取值范圍為0~1。

賦予訓練集標簽后,利用其訓練LSTM-SA 模型。網絡隨機初始化參數,根據預設的誤差函數,計算網絡輸出與標簽之間的相對誤差,并朝著誤差減小的梯度方向訓練參數,直至誤差值降到預設的閾值以下或訓練次數達到預設的迭代次數為止。訓練過程中,將學習率初始化為較小的值,采用Adam 優化器訓練網絡并自適應調整學習率。

步驟3:利用最小二乘法多項式擬合健康度指標曲線[25],得到軸承的RUL。LSTM-SA 模型的輸出是軸承每個時刻的健康度指標,各時刻的健康度指標形成退化曲線。將非全壽命測試軸承信號輸入模型,得到非全壽命的健康度退化曲線,擬合該曲線并計算其達到失效閾值的時間點,該時間點就是測試軸承的壽命結束點。軸承的健康度指標取值在0~1 之間,0 表示軸承失效,1 表示軸承完全正常,因此將失效閾值設為0。

3 實驗數據驗證

3.1 數據描述

為了驗證所提方法的有效性,采用IEEE 協會提供的PHM 2012 數據集進行驗證[26]。該數據集包含利用加速度傳感器采集的17 組軸承的全壽命周期振動信號,其中6 組訓練集和11 組測試集,如表1所示。信號的采樣頻率為25.6 kHz,采樣間隔為10 s,單次采樣時間為0.1 s。采集設備共采集水平和垂直兩個方向的振動信號,根據文獻[27-28]的研究,水平方向的振動信號比垂直方向的振動信號提供的有效信息更多。因此,本文采用水平方向的振動信號進行驗證。

表1 軸承數據集描述Tab.1 Description of bearing datasets

3.2 滾動軸承RUL 預測

實驗首先對數據集進行退化特征提取,第一步需要確定包絡譜的分段數n。如前所述,將包絡譜分段的退化特征提取方法,可以有效提取軸承早期損傷的故障特征以及優化退化特征的單調性和趨勢性,這種性能改進與分段數n有關。因此,接下來將研究n對退化特征的單調性和趨勢性的影響。n的初始值分別設置為1,2,3,4,5,6,7 和8,利用主成分分析[29]對不同分段數所提取到的特征降維至一維,計算并比較主成分的單調性和趨勢性,選出最優的n。

趨勢性和單調性的計算分別如下式所示:

式中xi和yi分別表示時間和特征的值分別表示x和y的平均值。Trend的取值在0~1 之間,值越大,趨勢性越強。

式中N表示樣本的數量。Mono取值在0~1 之間,值越大,單調性越強。

除此之外,為了從整體上對單調性和趨勢性進行度量,定義Cori[30]作為退化特征度量指標,其值越大,說明退化特征越能反映軸承的退化趨勢。計算公式如下:

退化特征的Cori對比結果如圖4 所示。從圖4中可以看出,Cori隨著n的增加而增大,并在n=4 時達到最高點。之后,Cori隨著n的增加而減小。表2以軸承1_1 和2_2 為例,對比了n取不同值時,退化特征的Cori值。當分段數n太小時,會在退化特征所在頻段混入太多干擾信息;相反,當n太大時,會將退化特征所在頻段切碎,難以提取出單調性和趨勢性明顯的退化特征?;谝陨戏治?,在實驗中將包絡譜分段數n設置為4。

表2 退化特征單調性和趨勢性對比結果Tab.2 Comparison results of monotonicity and tendency of degradation characteristics

圖4 退化特征的Cori 對比Fig.4 Comparison of Cori of degradation features

分段數n確定后,計算原始振動信號包絡譜,隨后劃分子頻段并分別計算各頻段與標準樣本的皮爾遜相關系數,將相關系數作為退化特征。將包絡譜劃分為n個子頻段的方法如下式所示:

式中s1,s2,…,sn分別表示n個子頻段的頻率范圍;F為包絡譜的分析頻率。將各子頻段按式(1)計算皮爾遜相關系數,得到退化特征。

在本實驗中,包絡譜分析頻率為12.8 kHz,如圖5 所示。將包絡譜按頻率劃分為4 個子頻段,每個子頻段的頻率范圍分別為0~3.2 kHz,3.2~6.4 kHz,6.4~9.6 kHz 和9.6~12.8 kHz,如圖6 所示。分別計算每個樣本各子頻段、全頻段及標準樣本相應頻段的皮爾遜相關系數,計算得到的相關系數作為各樣本的退化特征。利用上述退化特征訓練LSTM-SA 模型,從而構建健康度指標模型。

圖5 全頻段包絡譜Fig.5 Envelope spectrum for full band

為了進一步說明退化特征的單調性和趨勢性,以訓練數據軸承1_1 為例,利用t-SNE 對提取到的退化特征進行可視化分析,如圖7 所示。從圖7 中可以發現,所提取的特征隨著時間有序變化,能夠反映出軸承的退化過程。

圖7 t-SNE 降維可視化退化特征Fig.7 t-SNE reduced dimensional visualization of degradation features

本文按照PHM 2012 數據集的劃分,將軸承1_1,1_2,2_1,2_2,3_1 和3_2 全壽命周期數據作為訓練集,其余軸承作為測試集。訓練集提取退化特征后輸入到LSTM-SA 模型進行訓練。所提模型由2 層LSTM 層、1 層Self-Attention 運算層和2 層全連接層組成,所有前饋神經網絡的神經元個數均為64。初始學習率為0.005,隨機初始化權重參數,以均方誤差(MSE)作為損失函數,并利用Adam 優化器進行訓練。訓練和測試的CPU 環境為Intel Core i5,內存為16 GB,深度學習框架為Pytorch 1.10.1。

軸承在運行過程中逐漸發生退化,但其早期運行過程的退化狀態表現不明顯。以軸承1_7 和2_3為例,分別對兩組軸承全壽命時域波形圖進行分析,分別求各自的時域特征RMS,如圖8 和9 所示。從圖8 和9 中可以看出,軸承1_7 和2_3 在運行的前期和中期RMS 變化平穩,直到軸承損壞后期才出現RMS 的跳變。由此可知,軸承早期退化不明顯,難以捕捉其退化趨勢。

圖8 軸承1_7 原始振動信號均方根值Fig.8 RMS value of original vibration signal of bearing 1_7

圖9 軸承2_3 原始振動信號均方根值Fig.9 RMS value of original vibration signal of bearing 2_3

軸承1_7 和2_3 的健康度指標如圖10 和11 所示。從圖10 和11 中可知,軸承的健康度指標緩慢變化,較好地反映了軸承在運行過程中的退化狀態。由于存在局部振蕩,采用Savitzky-Golay 濾波器[31]對健康度指標進行平滑處理,消除局部振蕩影響并使健康度指標的退化趨勢更為平緩。由圖中健康度指標的變化趨勢可知,所提方法構建的健康度指標具有明顯單調性,能夠反映軸承的退化過程,且對軸承早期退化特征敏感。圖12 展示了訓練集中6 個軸承的健康度指標,針對軸承每個時刻的運行狀態,均給出0~1 之間的量化指標。

圖10 軸承1_7 的健康度指標Fig.10 Health indicator of bearing 1_7

圖11 軸承2_3 的健康度指標Fig.11 Health indicator of bearing 2_3

圖12 訓練集軸承的健康度指標Fig.12 Health indicators of bearings in training set

為驗證所提方法構造的健康度指標對提高滾動軸承RUL 預測精度的作用,利用最小二乘法多項式擬合由健康度指標得到的軸承退化曲線,預測軸承的失效點,從而得到軸承的RUL。選取三次多項式進行退化曲線擬合,其公式為:

式中y表示健康度指標;k表示第k個樣本,k=1,2,…,N;a,b,c和d表示待擬合的系數。

設t為軸承當前運行時間,當所預測的健康度指標達到失效閾值0 時,求對應失效時間t',二者之差即為預測的剩余使用壽命:

以軸承1_7 和2_3 為例,圖13 和14 分別為軸承1_7 和2_3 的剩余使用壽命預測結果。已知數據集中軸承1_7 當前壽命為15020 s,從圖13 可知,預測的失效時刻為20160 s,由式(19)計算得到RUL 的預測值為5140 s,而真實剩余使用壽命為7570 s。同理,已知軸承2_3 的當前壽命為12010 s,從圖14可知,失效時刻為18560 s,則預測RUL 為 6550 s,而真實剩余使用壽命為7530 s。

圖13 軸承1_7 的RUL 預測結果Fig.13 RUL prediction results of bearing 1_7

圖14 軸承2_3 的RUL 預測結果Fig.14 RUL prediction results of bearing 2_3

3.3 模型評估與誤差分析

本文采用誤差Ei來評估模型的預測精度,如下式所示:

式中 actRULi表示第i個測試軸承的真實剩余使用壽命;predRULi表示第i個測試軸承的預測剩余使用壽命。

除評估模型的預測精度之外,還需建立一個總體指標來評估模型的有效性,即是否超前預測或滯后預測。因此,本文采用PHM 2012 數據集設立的模型平均得分Score來評估模型的有效性。計算方法如下式所示:

式中Ai表示第i個軸承的得分。

圖15 展示了誤差Ei與得分Ai之間的關系。從圖15 中可以看出,當超前預測,即Ei>0 時,模型的得分更高;反之,當滯后預測,即Ei≤0 時,模型的得分較超前預測低。這是因為在實際生產過程中,超前預測比滯后預測更有意義,因此該得分計算方法對滯后預測進行了懲罰,具有公平性與合理性。

圖15 誤差Ei 與得分Ai 的函數關系Fig.15 Function relationship of error Ei and score Ai

表3 列出了測試集中11 個軸承的平均絕對預測誤差以及模型平均得分,采用Transformer Encoder模 型[7]、SOM模型[32]以 及CNN-LSTM模型[33]進 行對比實驗。其中文獻[7]和[32]先提取退化特征建立健康度指標,再進行RUL 預測,文獻[33]屬于“端到端”的RUL 預測方法。從對比結果可以看出,本文提出的方法相較于文獻[7],[32]和[33],平均絕對誤差分別降低了43.18%,62.57%和59.44%,平均得分分別提高了10.87%,45.71%和34.21%。此外,為了進一步說明利用包絡譜分段并計算皮爾遜相關系數提取退化特征的有效性,需要增加對比實驗,即不對包絡譜分段,對整個包絡譜計算皮爾遜相關系數提取退化特征,并進行RUL 預測,對比結果如表4 所示。由表4 可知,包絡譜分段比不分段的平均絕對預測誤差降低了35.56%,Score提高了49.02%。

表3 RUL 預測結果與比較Tab.3 RUL prediction results and comparisons

表4 對比實驗結果Tab.4 Comparison of experimental results

4 工程實際數據驗證

本節將所提方法應用于實際現場設備數據,進一步驗證方法的有效性。該數據來源于某石化企業離心泵,該泵為交流異步電機驅動、滾動軸承雙支撐結構,額定轉速2980 r/min。如圖16 所示,泵兩端的軸承座(3#軸承和4#軸承)上均安裝了IEPE 振動加速度傳感器,型號為PCB 608A11,線性頻響范圍為0.5~10 kHz,量程為±50g。圖16 中,Ha 表示水平方向的加速度;Va 表示垂直方向的加速度。每個軸承均采集了垂直和水平方向的振動信號。每個加速度傳感器的采樣頻率均被設定為25.6 kHz,采樣間隔為10 s,單次采樣時間為0.64 s。在監測的過程,泵兩端的軸承發生了退化,最終軸承徹底損壞。圖17 為3#軸承振動信號波形圖,從圖17 中可以看出,信號的幅值在軸承損傷的最后階段隨著時間而增加。

圖16 泵的測點布局圖Fig.16 Survey point layout diagram of a pump

圖17 3#軸承時域振動信號Fig.17 Time domain vibration signal of 3# bearing

在現場設備數據的驗證過程中,采用了與第3 節相同的操作流程和相同的n值,在相同的實驗環境下進行訓練和驗證。為了突出LSTM-SA 模型的優勢,文中用基于標準遞歸神經網絡(Standard-RNN)和基于卷積神經網絡(CNN)的RUL 預測模型進行比較。圖18為本文提出的基于LSTM-SA 模型的RUL 預測結果;圖19 和20 分別為基于Standard-RNN 和CNN的RUL 的預測結果。結合表5 三種模型的比較結果可以看出,LSTM-SA 模型的平均預測誤差分別比Standard-RNN 和CNN 模型降低了39.58% 和74.86%。該結果也進一步說明LSTM-SA 模型可以有效地預測滾動軸承的剩余使用壽命。

表5 RUL 預測結果與比較Tab.5 RUL prediction results and comparisons

圖18 LSTM-SA 模型的RUL 預測結果Fig.18 RUL prediction results of LSTM-SA model

圖19 Standard-RNN 模型的RUL 預測結果Fig.19 RUL prediction results of Standard-RNN model

圖20 CNN 模型的RUL 預測結果Fig.20 RUL prediction results of CNN model

5 結論

RUL 的預測精度很大程度上依賴于所構建的HI。所提出的方法結合包絡譜分段特征和LSTM-SA 模型提高了滾動軸承RUL 的預測準確率。在HI 的構建過程中,提出了基于包絡譜特征的退化特征提取方法。將退化特征輸入到LSTM-SA中構建HI。為了驗證所提出方法的有效性,使用公開的實驗數據和真實現場數據進行驗證。在實驗數據集的驗證中,所提出方法比文獻[7],[32]和[33]中的方法表現更好。此外,真實現場數據的驗證表明,所提出方法比基于Standard-RNN 和基于CNN的方法更能有效地預測泵軸承的RUL。

猜你喜歡
皮爾遜頻段軸承
軸承知識
軸承知識
軸承知識
軸承知識
5G高新視頻的雙頻段協同傳輸
gPhone重力儀的面波頻段響應實測研究
現代統計學之父:卡爾·皮爾遜
現代統計學之父:卡爾·皮爾遜
Excel在水文學教學中的應用
卡方分布的探源
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合