90°彎管存在下的供水管道泄漏定位研究

文 靜， 張敏姿， 張 恒

(重慶大學 光電工程學院 傳感器與儀器系統(tǒng)研究中心，重慶 400044)

供水管道泄漏不僅造成水資源流失，還給人們的生活環(huán)境帶來安全隱患。及時發(fā)現(xiàn)并準確定位泄漏是至關重要的?；谛孤┞曊駝訖z測的相關定位法是目前被廣泛采用的泄漏檢測定位方法[1-4]。該方法利用疑似漏點兩端管道振動信號的時延值來定位漏點。

在相關檢測定位法的檢測信號模型中，兩路信號所包含的非泄漏噪聲被假設為是互不相關的，若被測管段存在彎管，由彎管引發(fā)的聲信號[5]也會沿管道傳輸，從而在檢測信號中產(chǎn)生相關的非泄漏噪聲，此時利用相關檢測定位法將無法實現(xiàn)準確定位。文獻[6]通過研究不同源信號在管道傳播介質(zhì)中的衰減特征，來識別互相關峰值是來自泄漏還是非泄漏聲源，但該方法僅僅用于泄漏辨識，無法解決漏點準確定位問題。文獻[7]提出利用固定干擾源與泄漏信號源之間的獨立性實現(xiàn)固定干擾源存在下的管道泄漏檢測定位；而彎管與泄漏信號都是由流體激勵管道而引發(fā)的振動信號，基于流體的連續(xù)性，兩信號源不滿足相互獨立的關系。因此，有必要對彎管源產(chǎn)生機理進行詳細的分析，進而研究彎管存在下的管道泄漏準確定位方法。

本文利用數(shù)值仿真分析軟件ANSYS分別對泄漏和彎管進行模擬，在無噪聲影響下對泄漏和彎管引發(fā)的流致振動信號進行分析，在此基礎上提出一種基于EMD(Empirical Mode Decomposition)分解，以樣本熵為特征值進行篩選的彎管噪聲抑制方法，并通過仿真和實驗測試，驗證了該方法的有效性。

1 泄漏與彎管引發(fā)的管道振動仿真

1.1 計算方法

泄漏和彎管引起的管道振動都是流固耦合過程。文獻[8]比較了ANSYS軟件中不同湍流數(shù)值模擬方法，指出對90°圓形截面彎管內(nèi)的流體使用大渦數(shù)值模擬能更準確地反映出流體的流動狀態(tài)。因此，本文在仿真計算過程中，采用基于流固耦合控制方程和大渦模擬控制方程的數(shù)值方法進行研究。

流體流動要遵循物理守恒定律，在不考慮能量傳遞的情況下，需要遵守質(zhì)量守恒定律和動量守恒定律。對于一般的可壓縮牛頓流來說守恒定律通過N-S控制方程實現(xiàn)

(1)

(2)

式中：t為時間；f為體積力矢量；ρ為流體密度；v為流體速度矢量；p為流體壓強；μ為流體的動力黏度。固體部分的控制方程由牛頓第二定律導出

(3)

大渦模擬是一種典型的湍流數(shù)值模擬方法[9]， 它用瞬時的N-S方程直接模擬湍流中的大尺度渦，小尺度渦則通過大尺度渦的瞬時N-S方程體現(xiàn)出來，對瞬時N-S方程進行過濾，可得到LES控制方程

(4)

(5)

式中：τij為亞網(wǎng)格應力，定義為

(6)

(7)

式中：νt為亞格子渦黏系數(shù)。

1.2 計算模型

本文討論的管道幾何尺寸如圖1所示。

直泄漏管內(nèi)部流體介質(zhì)為水和水蒸氣的混合相，彎管內(nèi)部流體介質(zhì)為單相水，管道材質(zhì)為結(jié)構鋼，具體材料屬性如表1所示。

其中：ρ為物質(zhì)的密度；E為材料的楊氏模量；μ為材料的泊松比；d為管道的管壁厚度；D為管道的內(nèi)徑。彎管和直泄漏管的固體兩端都采用固定約束， 即ux，uy，uz都為零， 流體以1 m/s的速度從管道入口流入，管道的出口端為流體自由流出，泄漏管道的泄漏口為0 Pa的壓力出口。

圖1 管道模型幾何尺寸圖Fig.1 The geometry of pipe model

水水蒸氣結(jié)構鋼壁厚管內(nèi)徑ρ/(kg·m-2)ρ/(kg·m-2)ρ/(kg·m-2)EμdD998．20．554278502000．36219

1.3 泄漏管內(nèi)與彎管內(nèi)流體輸出分析

1.3.1 泄漏管內(nèi)流體輸出

在泄漏管道中，泄漏口的流體運動是影響管道振動的主要因素。要獲取真實的管道振動，就需要對泄漏口進行準確的流場模擬。泄漏口處，管壁內(nèi)外存在很大的壓差，又因為壁厚在毫米量級，導致泄漏口處的壓力梯度極大。流體在極大的壓差條件下，壓力對流體做功，使流體的動能增加，以噴射狀流出管外。圖2所示為泄漏口處的流體速度矢量圖。壓力與動能的轉(zhuǎn)換，致使泄漏口壁面出現(xiàn)負壓，這與文獻[10]中對泄漏口的數(shù)值仿真結(jié)果一致，泄漏口處的壓力云圖如圖3所示。當管內(nèi)的壓強下降到一定的臨界值，就會發(fā)生相的變化，出現(xiàn)空化現(xiàn)象。從圖3可以看出，泄漏口壁面附近的最低壓強已低至-1.223×105Pa，早已低于空化的臨界點，在泄漏口管壁附近出現(xiàn)空泡，泄漏口處的氣體體積分數(shù)云圖如圖4所示。

圖2 泄漏口處流體速度矢量圖Fig.2 The fluid velocity vector of leakage hole

圖3 泄漏口處的壓力云圖Fig.3 The pressure of leakage hole

圖4 泄漏口處的氣體體積分數(shù)云圖Fig.4 The gas volume fraction of leakage hole

泄漏處的聲源主要由三部分組成：泄漏空化聲、泄漏湍流聲和泄漏時流體與管壁的摩擦聲[11]。泄漏湍流聲和泄漏摩擦聲都是流體與固體間力的相互作用的體現(xiàn)，ANSYS平臺下的流固耦合對于這些相互作用力的傳遞都能很好的捕捉到。對于泄漏空化聲，不少的空化噪聲實驗發(fā)現(xiàn)[12]：空化輻射的噪聲與空泡中心到壁面的無量綱距離γ有關，γ=0.9時，其輻射的噪聲會特別小。因為當γ足夠大時，空泡潰滅時對應泡的體積發(fā)生劇烈的變化，從而產(chǎn)生很大的輻射噪聲。當γ=0.9左右時，在空泡體積變化之前形成了一股射流，射流使泡的上下表面接觸并形成一個渦環(huán)，這樣大部分的泡能轉(zhuǎn)換成渦能，泡的體積沒有發(fā)生很大的變化，所以輻射的聲音減小很多。

由圖4可知，空泡分布在靠近泄漏口管壁兩側(cè)的地方，對于泄漏口壁面附近的空泡，雖然不能定量衡量γ值的大致范圍，不確定是否會有相應的射流產(chǎn)生使泡能轉(zhuǎn)換成渦能，但從圖2的速度矢量圖可觀察出，流體在泄漏的邊界條件下，自身就會產(chǎn)生射向空泡方向的射流，同樣導致空泡上下表面接觸形成渦環(huán)，大大減小了空化噪聲。泄漏口處的渦能如圖5所示。由圖可知，管道內(nèi)的渦能量都集中到泄漏孔處了。這從一定程度上證明了泄漏口處的空泡能大量轉(zhuǎn)換成渦能的推論。

圖5 泄漏口處的渦能圖Fig.5 The eddy energy of leakage hole

綜上所述，無論壁面附近的空泡對應的γ為多少，泄漏空化聲都應是極小的，可以忽略不計。因此，泄漏引發(fā)的管道振動主要來源于湍流和摩擦。

1.3.2 彎管內(nèi)流體輸出

現(xiàn)有的對彎管內(nèi)流體的研究指出，彎管處的流體在外壁面壓力大速度小，內(nèi)壁面壓力小速度大；另外，彎管內(nèi)圓形截面上存在二次流動，這是由于離心力的作用而形成的一對反向?qū)ΨQ渦旋。彎管內(nèi)流體在各角度截面上的壓力云圖和速度流線圖分別如圖6、圖7所示，其中圖(a)～圖(e)分別為彎頭內(nèi)0°、30°、45°、60°、90°處的流體流動狀態(tài)。圖中所示結(jié)果反映出與文獻[13-18]描述相一致的彎管內(nèi)流體特性。

圖6 彎管處各角度的壓力云圖Fig.6 The pressure of 90° bend

圖7 彎管處各角度的速度流線圖Fig.7 The speed streamline of 90° bend

由此可知，采用這種雙向流固耦合仿真，能較準確的模擬出彎管內(nèi)流體的流動特性，以此作用于固體管道，激勵管道相應特征的振動，由此能真實地反映出實際彎管引發(fā)的管道振動特性。

1.4 泄漏與彎管引發(fā)的管道振動特性分析

基于以上仿真結(jié)果，提取管道振動的徑向加速度，對泄漏和彎管引發(fā)的管道振動特性進行分析。

1.4.1 泄漏管道振動特性分析

在圖1(a)的泄漏管模型上提取泄漏口上游和下游管壁上不同4點的振動信號進行分析。node-1，node-3位于泄漏口上游，離泄漏口距離分別為0.5 m和3 m；node-2，node-4位于泄漏口下游，離泄漏口距離分別為0.5 m和3 m，4個點的管壁振動信號功率譜如圖8所示。

圖8 泄漏管道仿真各點功率譜圖Fig.8 The power spectrum of each point on leaking pipeline

圖8表明，泄漏引發(fā)的管道振動頻率主要集中在低頻段500 Hz左右，高頻段的信號能量近乎為0。在node-1和node-2兩點處捕捉到了2 000 Hz左右的高頻信號，此高頻信號疑為少量的空泡噪聲，由于高頻信號衰減快的特性，在node-3、node-4處的功率譜中，近乎觀察不到2 000 Hz左右的高頻信號。

1.4.2 彎管振動特性分析

在圖1(b)的彎管模型上同樣提取上游和下游管壁上不同4點的振動信號進行分析。node-1，node-2位于彎頭上游，離彎頭入口距離分別為0.8 m和1.5 m；node-3，node-4位于彎頭下游，離彎頭出口距離分別為0.8 m和2 m。4個點的管壁振動信號功率譜如圖9所示。

圖9表明，彎頭上游的管道振動能量在各頻段上的分布是比較均勻的，而彎頭下游管道振動的能量主要集中在1 000 ～2 500 Hz頻段。這是因為流體在彎頭的作用下，被迫改變流動方向后，彎頭的下游一定距離內(nèi)，流體在主流的流動狀態(tài)下疊加上二次流形成螺旋流。這樣一種流動狀態(tài)，會在很長的距離上極緩慢的消失[19]。node-3和node-4剛好是彎頭下游不遠處的2點，螺旋流對管道的作用直接被這2點的振動信息捕捉到，而沒有經(jīng)過管道傳遞的衰減，于是node-3和node-4的功率譜在1 000～2 500 Hz頻段上的能量要明顯大于node-1和node-2的相應頻段上的能量。

圖9 彎管仿真各點的功率譜圖Fig.9 The power spectrum of each point on 90° bend

1.4.3 振動信號的樣本熵分析

樣本熵從衡量時間序列復雜性的角度來度量信號中產(chǎn)生新模式的概率大小，產(chǎn)生新模式的概率越大，序列的復雜性就越大，相應的樣本熵值也越大。

由圖8和圖9對比觀察可發(fā)現(xiàn)，彎管信號的頻率分布要遠遠廣于泄漏信號的頻率分布。在流體受力分析方面，其引力受力情況為：泄漏時，流體只受重力影響；而彎頭內(nèi)的流體是在重力和離心力共同作用下運動的。綜合以上兩方面的分析，可提出彎管信號復雜度大于泄漏信號復雜度的假設。

從流體間的相互作用和流體與固體的相互作用分析來看：泄漏信號中，包含了泄漏空化、泄漏湍流和泄漏流體與管壁的撞擊與摩擦；彎管信號中，包含了流體撞擊拐角、彎頭處的二次流和邊界層分離等。以上的理論分析不能定量的描述泄漏信號和彎管信號的復雜度，于是本文分別在泄漏管模型和彎管模型上各選取4個點，對兩模型各點處采集的數(shù)據(jù)進行樣本熵求解，定量的分析兩種信號在復雜度上的區(qū)別。

樣本熵的計算方法可參考文獻[20]，其中關于m，r，N的選值，根據(jù)已有的研究，本文取值選為m=2，r=0.2Std，N=2 500。模型各點數(shù)據(jù)樣本熵如表2所示。

從表2可看出，直泄漏管道的樣本熵都在0.6左右，而彎管各點的樣本熵大致都為1.8左右，幾乎是管道泄漏信號的三倍，表明彎管信號的復雜度要遠遠大于泄漏信號的復雜度，可以此作為區(qū)別泄漏信號和彎管信號的特征值。

表2 不同模型各點樣本熵

2 用EMD方法去除彎管噪聲

經(jīng)驗模態(tài)分解(EMD)是依據(jù)數(shù)據(jù)自身的時間尺度特征來進行信號分解的，因此對于不同的信號具有很強的自適應性。它的關鍵在于能使復雜信號分解為有限個本征模態(tài)(Intrinsic Mode Function, IMF)，所有分解出來的各IMF分量包含了原信號的不同時間尺度的局部特征信號。EMD的分解過程其實是一個“篩分”的過程，在“篩分”的過程中，消除了模態(tài)波形的疊加，而且使波形輪廓更加明顯對稱。本文將對仿真采集的泄漏信號與彎管信號分別進行EMD分解，其原始信號與兩信號各自的分解結(jié)果如圖10所示。

圖10 泄漏與彎管的原始信號圖和EMD分解圖Fig.10 The original signal and it’s EMD decomposition signal of leaking pipeline and 90° bend

圖中10(a)為泄漏信號與彎管信號的原始信號，圖10(b)和圖10(c)分別為泄漏信號和彎管信號的EMD分解前6個IMF分量，由圖10可觀察到，各個IMF分量體現(xiàn)了不同的時間尺度特征，使得彎管和泄漏的特征分別在不同的時間尺度下顯露出來。另外，對比泄漏信號和彎管信號的IMF分量，可說明兩信號在時間特征尺度上有顯著的差異，因此通過EMD分解，能有效的將泄漏信號與彎管信號從兩者的混合信號中分離開。

根據(jù)以上分析可得出如下結(jié)論:對于混雜了彎管信息的泄漏信號，先對信號做EMD分解，根據(jù)泄漏信號和彎管信號各自不同的時間尺度特征將泄漏信號和彎管信號分別分解到不同的IMF中，又因為兩信號的樣本熵存在明顯的不同，于是可用樣本熵為特征值從IMFs中篩選出一定閾值內(nèi)的IMF分量進行重構，從而去除泄漏信號中的彎管噪聲。

3 仿真實驗與實際測試結(jié)果

3.1 仿真實驗

在ANSYS平臺下，對泄漏和彎頭同時存在的管道振動信號進行模擬。管道模型如圖11所示。管內(nèi)流體介質(zhì)為水和水蒸氣的混合相，管材料為結(jié)構鋼，具體材料屬性如表1所示。

圖11 泄漏與彎頭同時存在的管道剖面圖Fig.11 The geometry of leakage along with 90° bend pipe model

采集泄漏口上游1 m處node-1和彎頭下游2 m處node-2的管道徑向加速度振動信號，分別記為x1(t)和x2(t)， 對x1(t)和x2(t)做互相關處理，其互相關譜圖如圖12所示。

圖12 x1(t)和x2(t)的互相關譜圖Fig.12 The cross-correlation shifting of x1(t) and x2(t)

從圖12中能觀察到兩個明顯的峰值，這是由泄漏和彎頭兩個振動源共同作用的結(jié)果。據(jù)經(jīng)驗值，泄漏信號的頻率不低于200 Hz，對采集得來的x1(t)和x2(t)分別進行EMD分解，將低于200 Hz的IMF分量直接舍去。于是，x1(t)和x2(t)分解后的各IMF分量，只提取了前6個。分別對x1(t)和x2(t)的前6個IMF分量求樣本熵，相應的樣本熵值如表3所示。

表3 x1(t)和x2(t)的6個IMF分量的樣本熵

圖和的互相關譜圖Fig.13 The cross-correlation shifting of and

比較圖12和圖13，可明顯觀察到，兩路信號的互相關在沒有去除彎管噪聲的時候存在著兩個峰值，但經(jīng)過EMD和樣本熵結(jié)合去噪的處理后，就僅有單一的明顯峰值了，從而說明EMD和樣本熵結(jié)合的方法能有效地去除信號中的彎管噪聲。

3.2 實際測試結(jié)果

實際測試選取直徑為100 mm的鑄鐵管道，用消防栓流水模擬泄漏點，聲速取1 300 m/s。實驗環(huán)境選取的是周圍噪聲很小的校園環(huán)境，這樣才能使得采集信號中主要包含泄漏噪聲和彎管噪聲。采用壓電加速度傳感器進行實際信號采集，傳感器的參數(shù)如下：頻率響應區(qū)間為0.1～2.5 kHz，靈敏度為5.982 V/g，采用磁性基座便于其吸附在管道外壁上，用橡膠套做防水處理，傳感器輸出信號利用電磁屏蔽電纜與其他部件相連。壓電傳感器采集信號時的安裝實物圖和實際采集管道模型分別如圖14(a)、圖14(b)所示。

在圖14(b)所示的模型上分別在兩個不同的位置進行兩次測量，A點采集的信號記為x1(t)，B點采集的信號記為x2(t)。

圖14 實驗信號采集示意圖Fig.14 The model of experimental signal acquisition

圖15 第一次測試數(shù)據(jù)互相關譜圖Fig.15 The cross-correlation shifting of first test data

圖16 第二次測試數(shù)據(jù)互相關譜圖Fig.16 The cross-correlation shifting of second test data

從實驗的兩組對比圖可看出，EMD分解與樣本熵的結(jié)合使用，能有效的抑制彎管噪聲引起的多相關峰，突出泄漏信號的相關峰。經(jīng)過去除彎管噪聲的處理后，第一次測量數(shù)據(jù)的定位誤差為0.4 m，第二次測量數(shù)據(jù)的定位誤差為0.3 m，誤差均在1 m以內(nèi)，從而說明該方法在一定程度上能有效的抑制彎管噪聲對泄漏定位的影響。

4 結(jié) 論

本文提出了一種針對金屬材質(zhì)管道的基于EMD與樣本熵的泄漏信號去彎管噪聲方法，對仿真和實際振動信號進行分析表明：

(1) 泄漏引發(fā)的管道振動信號的樣本熵值在0.5～0.7內(nèi)，而彎管誘發(fā)的管道振動信號的樣本熵值在1.7～2內(nèi)，彎管信號的樣本熵要遠遠大于泄漏信號的樣本熵，說明了樣本熵作為區(qū)分泄漏信號與彎管信號的特征值的可靠性。

(2) 對泄漏信號和彎管信號做EMD分解，對比兩信號的各IMF分量，發(fā)現(xiàn)兩者存在明顯不同的時間尺度特征，這樣的顯著差異表明，EMD分解能從混合信號中自適應的分解出分別含有泄漏信號和彎管信號的IMFs。

(3) 結(jié)合EMD分解與樣本熵特征篩選，能有效的從含有彎管噪聲的信號中提取出泄漏信號，從而解決了彎管噪聲引起的多個相關峰的定位不確定問題。

[ 1 ] FUCHS H V, RIEHLE R. Ten years of experience with leak detection by acoustic signal analysis[J]. Applied Acoustics, 1991, 33(1): 1-19.

[ 2 ] GAO Y, BRENNAN M J, JOSEPH P F. A comparison of time delay estimators for the detection of leak noise signals in plastic water distribution pipes[J]. Journal of Sound & Vibration, 2006, 292(3): 552-570.

[ 3 ] ALMEIDA F, BRENNAN M, JOSEPH P, et al. On the acoustic filtering of the pipe and sensor in a buried plastic water pipe and its effect on leak detection: an experimental investigation[J]. Sensors, 2014, 14(3): 5595-5610.

[ 4 ] LI S, WEN Y, LI P, et al. Leak location in gas pipelines using cross-time-frequency spectrum of leakage-induced acoustic vibrations[J]. Journal of Sound & Vibration, 2014, 333(17): 3889-3903.

[ 5 ] 李麗, 沈壯志, 沈建中, 等. 90°彎管的流場中渦相互作用對聲場的影響[J]. 陜西師范大學學報(自然科學版), 2014, 42(5): 27-31.

LI Li, SHEN Zhuangzhi, SHEN Jianzhong, et a1. The interactions between vortexes influence on the sound field in the flow field of 90° curved tube [J]. Journal of Shaanxi Normal University (Natural Science), 2014, 42(5): 27-31.

[ 6 ] 楊進, 文玉梅, 李平, 等. 非泄漏固定聲源干擾下的管網(wǎng)泄漏定位技術研究[J]. 儀器儀表學報, 2012, 33(2): 248-254.

YANG Jin, WEN Yumei, LI Ping, et al. Research on leak location in the presence of non-leak noises in water distribution pipelines [J]. Chinese Journal of Scientific Instrument, 2012, 33(2): 248-254.

[ 7 ] 文玉梅, 由原, 楊進, 等. 固定干擾源存在下的管道泄漏檢測和準確定位[J]. 聲學學報, 2013, 38(1): 80-84.

WEN Yumei, YOU Yuan, YANG Jin, et al. Leak detection and accurate location in the presence of fixed source of interference [J]. Acta Acustica, 2013, 38(1): 80-84.

[ 8 ] 龔斌, 劉喜興, 楊帥, 等. 90°圓形截面彎管內(nèi)流動的大渦模擬[J]. 過程工程學報, 2013, 13(5): 760-765.

GONG Bin, LIU Xixing, YANG Shuai, et al. Simulation on large eddy turbulent flow in a circular-sectioned 90°bend [J]. The Chinese Journal of Process Engineering, 2013, 13(5): 760-765.

[ 9 ] 王福軍. 計算流體動力學分析[M]. 北京：清華大學出版社, 2004.

[10] 楊進. 供水管道泄漏檢測定位中的信號分析及處理研究[D]. 重慶：重慶大學, 2007.

[11] 馮雪松, 文玉梅, 甄錦鵬, 等. 管道泄漏聲振動信號的特征分析[J]. 聲學技術, 2015, 34(5): 413-418.

FENG Xuesong, WEN Yumei, ZHEN Jinpeng, et al. Feature analysis of pipeline leakage acoustic signals for leak identification [J]. Technical Acoustics, 2015, 34(5): 413-418.

[12] 戚定滿. 空泡潰滅及空化噪聲研究[J]. 力學季刊, 1999(1): 1-9.

QI Dingman. Empty bubble breaking and cavitation noise research [J]. Chinese Quarterly of Mechanics, 1999(1): 1-9.

[13] HURST J M. Note on the motion of fluid in a curved pipe[J]. Mathematika, 1959, 6(1): 77-85.

[14] SUDO K, SUMIDA M, HIBARA H. Experimental investigation on turbulent flow in a circular-sectioned 90-degree bend[J]. Experiments in Fluids, 1998, 30(3): 246-252.

[15] KIM J, YADAV M, KIM S. Characteristics of secondary flow induced by 90-degree elbow in turbulent pipe flow[J]. Engineering Applications of Computational Fluid Mechanics, 2014, 8(2): 229-239.

[16] R?HRIG R, JAKIRLI S, TROPEA C. Comparative computational study of turbulent flow in a 90° pipe elbow[J]. International Journal of Heat & Fluid Flow, 2015, 55: 120-131.

[17] 江山, 張京偉, 吳崇健, 等. 基于FLUENT的90°圓形彎管內(nèi)部流場分析[J]. 中國艦船研究, 2008, 3(1): 37-41.

JIANG Shan, ZHANG Jingwei, WU Chongjian, et al. Numerical simulation of inner flow in 90° bending duct of circular-section based on FLUENT[J]. Chinese Journal of Ship Research, 2008, 3(1): 37-41.

[18] 湛含輝, 朱輝, 陳津端, 等. 90°彎管內(nèi)二次流(迪恩渦)的數(shù)值模擬[J]. 鍋爐技術, 2010, 41(4): 1-5.

ZHAN Hanhui, ZHU Hui, CHEN Jinduan, et al. Numerical simulation of secondary flow (dean vortices) in 90° curved tube [J]. Boiler Technology, 2010, 41(4): 1-5.

[19] 董志勇. 彎頭(彎管)阻力系數(shù)比較與流動特性分析[C]//全國水動力學研討會. 北京：北京海洋出版社，2001.

[20] PINCUS S M. Approximate entropy as a measure of system complexity[J]. Proceedings of the National Academy of Sciences of the United States of America, 1991, 88(6): 2297-3301.