• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    MorPhological Undecimated Wavelet DecomPosition Fusion Algorithm and Its APPlication on Fault Feature Extraction of Hydraulic PumP

    2015-11-21 07:56:25SunJian孫健LiHongru李洪儒WangWeiguo王衛(wèi)國YePeng葉鵬
    關(guān)鍵詞:孫健

    Sun Jian(孫健),Li Hongru(李洪儒),Wang Weiguo(王衛(wèi)國),Ye Peng(葉鵬)

    Mechanical Engineering College,Shijiazhuang 050003,P.R.China

    MorPhological Undecimated Wavelet DecomPosition Fusion Algorithm and Its APPlication on Fault Feature Extraction of Hydraulic PumP

    Sun Jian(孫健),Li Hongru(李洪儒)*,Wang Weiguo(王衛(wèi)國),Ye Peng(葉鵬)

    Mechanical Engineering College,Shijiazhuang 050003,P.R.China

    Since vibration signals of hydraulic pump are mostly nonlinear and traditional fusion algorithm cannot satisfyingly process them,a morphological undecimated wavelet decomposition fusion(MUWDE)algorithm is proposed.Eirstly,under the framework of morphological undecimated wavelet decomposition(MUWD),multichannel signals are decomposed.Approximate signals of all decomposition layers are selected by feature energy factor and fused according to the presented fusion rules.Eurthermore,specific method for optimal selection of MUWDE parameters is presented to avoid subjective influences.Einally,the proposed algorithm is verified by simulation signals and pump vibration signals.

    morphological undecimated wavelet decomposition(MUWD);information fusion;feature extraction;hydraulic pump

    0 Introduction

    Mechanical diagnosis and prognostics require effective information detected by multi-sensors to extract precise fault feature.Vibration signals of hydraulic pump caused by early fault are mostly nonlinear[1]due to structural characters.And the feature information ingredient is relatively weak. As a result,appropriate fusion algorithm is required for extracting precise fault feature effectively.

    Currently,various fusion algorithms have been proposed for mechanical feature extraction[2-3].However,weighted fusion algorithm is subjected to fusion weights[4-5].Kalman filtering algorithm lacks strict filtering functions for nonlinear system[6-7].The wavelet analysis fusion algorithm may neglect the feature information during the sampling operation[8-9].Consequently,fusion performance based on conventional algorithm is not satisfying and the fault feature is hard to effectively extract,which fails to meet the requirements of precise diagnosis and prognostics.

    As an efficient nonlinear signal processing method,morphological undecimated wavelet decomposition(MUWD)omits the sampling operation during decomposition and reconstruction,which effectively avoids distortion[10-11].MUWD has been widely used for mechanical signal processing[12-13].However,in conventional MUWD,the approximate signal in the highest decomposition layer is considered as result.Others are neglected,thus feature information is missed.Each sensor typically sends one bit of information about fault feature rather than send the observed requirements[14].Therefore,an appropriate fusion rule is required to modify conventional MUWD.

    Therefore,Based on the framework of conventional MUWD,the novel morphological undecimated wavelet decomposition fusion(MUWDE)algorithm is proposed.Multi-channel signals are respectively decomposed firstly.Approximate signals in decomposition layers are selected according to feature energy factor(EEE).Based on the presented fusion rules,the feature information is generally used and the signals are reconstructed.Eurthermore,we use the fusion energy entropy(EEE)to measure fusion performance of MUWDE with various combinations of parameters,and demonstrate the optimization method in detail.Simulation and experimental results show the validity and superiority of the proposed algorithm.

    1 MUWDF Algorithm

    1.1 Introduction of MUWDF

    Assume that collections of V and W are the i th signal space and i th detail space,respectively,and T(·)denotes the morphological operator. Eramework of conventional MUWD can be describe by Eqs.(1—3)[15]

    where Eq.(1)is to use signal analysis operator φ↑ito decompose signal from Vito Vi+1;Eq.(2)is to use detail analysis operatorω↑ito decompose signal from Vito Vi+1;Eq.(3)denotes that the composition operatorand the analysis operators meet wavelet reconstruction and pyramid principle,which avoid information missing during signal analysis and composition[16-17].

    Through the framework analysis,we can identify T(·)as the key calculation of MUWD[18].Normal ones are mean combination operator,morphological difference operator,morphological gradient operator and hybrid operator[19-20].Considering pump structural characters,the morphological difference operator is selected.It contains both black and white Top-Hat,which can effectively extract signal pulses[21].Therefore,basic calculation of MUWDE can be described as

    where f(n)is the initial signal and g(n)the structure element;"and"·denote operations of morphological open and morphological close;(i+1)g means the i th swelling operation on flat element.Based on operator,fault feature information is extracted from the former layer and restored in current approximate signal,which is to be the initial signal for decomposition next.

    1.2 Fusion of multi-channel signals based on MUWDF

    Consider that X1—X3respectively refer to the signal in Channel1—Channel3 probed by Sensor1—Sensor3.Suppose that the number of decomposition layers is N.

    Signals are decomposed by Eqs.(4—6). Then approximate signals in all N decomposition layers of X1are obtained as,…,N};for X2,the approximate signals arej=1,…,N};and for X3,they are x3=.To measure various contributions for feature information ofandEEE is defined.EEE can be described as the ratio between partial energy of the prior h doubling frequency of feature frequency and the total energy in frequency domain.Its expression is described as

    where Ehdenotes the energy of the h th feature frequency.In this article,we set h=3.Obviously,EEE can better reflect the feature information proportion.The higher the EEEis,the larger the proportion is.Therefore,for the i th layer,we can respectively calculateofofandofby Eq.(7).Consequently,fusion rule for signals in thesame decomposition layer is demonstrated as follows.

    Fusion Rule 1 Eor the j th decomposition signalsthe one with the maximum EEE will be selected.

    Based on fusion rule 1,parts of approximate signals of X1remain and parts of X2and X3are absorbed as well.Therefore,feature redundant information and complementary information are utilized initially.We can select approximate signals in all decomposition layers of X1—X3according to fusion rule 1.Define the collection of selected approximate signals C={y1,y2,…,yN},and the collection W={w1,w2,…,wN}of EEE values.Means for calculating fusion weights kiis described as

    Consequently,selected approximate signals can be fused by

    Therefore,the second fusion rule is:

    Fusion Rule 2 Eor the selected approximate signals in collection C,calculate the corresponding fusion weights and reconstruct the signal according to Eq.(9).

    Based on fusion rule 1,fusion rule 2 realizes abundant use of feature information and reconstructs the signal by fusion of approximate signals for achieving the global decision.

    2 OPtimization of Parameters

    According to the theory of MUWDE,two parameters are involved.One is N,the number of decomposition layers;the other is L,the initial length of flat structure element.

    If N is too large,information in detail signal will be extremely few,so that the ingredient in approximate signal almost stays invariant.

    If N is too small,information in detail signal will be too much,so that the extracted feature will be unclear.

    If L is too large,extracted pulses will decrease,leading to the loss of feature information.

    If L is too small,noises may be retained in the reconstruction signal.

    Therefore,specific method for optimal selection of parameters of N and L will be detailed.

    Eirstly,ranges of N and L are supposed to be confirmed.Suppose that lengths of signal X1—X3are Len1—Len3.After searching the whole sequence of X1,the local maximum collection is LY1={LY1ii=1,…,NLY1}and the minimum collection is LV1={LV1ii=1,…,NLV1},where NLY1and NLV1are the number of local maximum dots and minimum dots,respectively.Similarly,for X2and X3,we can get collections LY2,LV2,LY3and LV3.Correspondingly,collections of distances of local maximum are DY1={d Y1i= LY1i+1-LY1ii=1,…,NLY1-1},DY2and DY3. Collections of distances of local minimum are DV1={d V1i=LV1i+1-LV1ii=1,…,NLV1-1},DV2and DV3.Then the range of N can be calculated by Eqs.(10—12).

    where fsdenotes the sampling frequency and f0the feature frequency.Calculation for the number of combination of N and L is

    In order to search for the optimal combination,we present EEE to evaluate fusion performance.Eor each combination(suppose N=n,L= l),consider that Eidenotes the energy of approximate signal yiselected by fusion rule 1.Then EEE is defined as

    According to Eq.(16),we can know that EEE can describe signal complexity and fusion performance.The lower the EEE is,the simpler the signal ingredient will be and the more obvious the feature information will be.Consequently,by comparing values of EEE of each combination,we can select the one with the minimum EEE to be the optimum.

    Above all,the fusion process of multi-channel vibration signals of hydraulic pump is described in Eig.1.

    Eig.1 Eusion process of multi-channel hydraulic pump vibration signals

    3 Simulation Analysis

    In order to verify the fusion performance of MUWDE,simulation signals are analized.Consider that fs=1 024 Hz,and sampling time t=1 s.Simulation signals are y1(t),y2(t)and y3(t)

    It is clear that the signals mainly consist of three parts,fault signal x1(t),harmonic wave signal x2(t)and white noise n(t).x1(t)means the simulation of pulse signal with periodic index decrease because of rolling bearing fault.The striking frequency is f0=16 Hz,the striking function in every period is e-200tsin(2π×256t)and the resonance frequency is 256 Hz.x2(t)represents the harmonic wave signal cos(2π×40t)+ cos(2π×50t),which contains the frequencies of 40 Hz and 50 Hz.And the difference frequency is 10 Hz.y1,y2and y3are various combinations of x1(t),x2(t)and n(t),which imitate the multichannel vibration signals probed by various sensors.Eigures of time domain and frequency domain of the simulation signals are shown in Eig.2 and Eig.3,respectively.

    In Eig.3,we can see that modulation in the resonance frequency 256 Hz exists in all signals. And fault feature frequency 16 Hz emerges by noises.Therefore,the proposed MUWDE algorithm is used to fuse y1,y2and y3for effective feature information extraction.

    (1)Optimization for MUWDE parameters N and L

    Based on the analysis of local maximums and minimums of the simulation signals,we can calculate Nmin=2 and Nmax=8 by Eqs.(10—12). Eor each N=2,3,…,8,there is a corresponding region for L.As fs=1 024 Hz and f0=16 Hz,various ranges of L can be achieved by Eqs.(13— 14),and listed in Table 1.

    Eig.2 Time domain of multi-channel simulation signals

    Table 1 Various combinations of Parameters N and L

    According to Eq.(15),there are 97 combinations of N and L.Each combination has its own EEE calculated by Eq.(16),which is shown in Eig.4.And Eig.5 is the amplification of Eig.5 in area 0—0.05.

    Eig.4 shows 97 values of EEE with various parameter combinations.Each curve refers to the changing of EEE along with L for every settled N.It is clear that each curve exists a minimum. On this base,we can find the minimum EEE for each N,which is{0.182 9,0.038 4,0.0132,0.008 9,0.004 2,0.007 7,0.008 1}.Eig.5 clearly shows the values of EEE ranging from 0 to 0.05.Conclusions can be drawn that the best fusion performance is achieved by MUWDE with combination of N=6 and L=7.The minimumEEE is 0.004 2.

    Eig.3 Erequency domain of multi-channel simulation signals

    Eig.4 EEE change with various combinations of N and L

    Eig.5 Amplification of EEE ranging from 0 to 0.05

    (2)Analysis of MUWDE fusion performance

    Eirstly,simulation signals y1,y2and y3are decomposed according to Eqs.(4—6).g0=[0 0 0 0 0 0 0]and N=6.According to Eq.(7),EEE of each decomposition layer of y1is{0.252 0,0.042 1,0.246 2,0.073 2,0.031 8,0.007 7}. Similarly,EEE of y2is{0.217 1,0.137 4,0.170 3,0.095 7,0.005 3,0.068 4}and EEE of y3is{0.303 2,0.109 8,0.061 8,0.084 1,0.009 4,0.028 1}.Then they are selected based on fusion rule 1.We take the first decomposition layer for example,and envelops of the corresponding approximate signals are shown in Eig.6.

    Eig.6 Approximate signals in the first decomposition layer of multi-channel signals

    Eig.6 shows the envelops of approximate signals in the first decomposition layer.In Eig.6(a,b),the feature frequency 16 Hz still emerges by noises.However,in Eig.6(c),16 Hz and its double frequency are relatively obvious.The phe-nomenon corresponds with EEE of these approximate signals y1—y3,which are 0.252 0,0.217 1 and 0.303 2.According to fusion rule 1,as 0.303 2>0.252 0>0.217 1,the approximate signal of y3in the first decomposition layer is remained.Similarly,remained approximate signals in other layers are 2nd of y2,3rd of y1,4th of y2,5th of y1and 6th of y2.The corresponding values of EEE are{0.303 2,0.137 4,0.246 2,0.095 7,0.031 8,0.068 4}.

    Eusion weights are achieved by Eq.(8),which are k1=0.343,k2=0.156,k3=0.280,k4= 0.108,k5=0.036,k6=0.077.Signal is reconstructed by fusing 6 approximate signals using Eq.(9):

    where yij(i=1,2,3 and j=1,2,…,6)denotes the approximate signal in the jthdecomposition layer of yi;yfinalthe reconstructed signal.Eusion result is shown in Eig.7.

    Eig.7 Eusion performance of multi-channel simulation signals based on MUWDE(N=6,L=7)

    Eig.7 shows that by the general use of feature information in multi-channel signals,the feature frequency 16 Hz and its doubling frequency are obvious.Eurthermore,the direct frequencies 40 Hz and 50 Hz are restrained and the modulation in 256 Hz is settled.To better indicate the advantages of MUWDE,we take conventional methods,such as weighted averaging fusion algorithm,Kalman filtering algorithm and wavelet analysis algorithm,to fuse the same multi-channel signals for comparison.Results are shown in Eig.8.

    Eig.8 Comparison of fusion performance of multichannel signals based on conventional fusion algorithms

    Eig.8(a)shows the fusion performance of multi-channel signals based on weighted averaging algorithm.It is clear that there are still harmonic waves and modulation ingredient,and the feature information is not effectively extracted. Eig.8(b)describes the fusion performance based on Kalman filtering algorithm.As the inherent disadvantages of the method,the fusion result is not satisfying.Although the harmonic waves and modulation are restrained,the feature information is still not clear and the difference frequency 10 Hz exists.Eig.8(c)indicates the fusion performance based upon wavelet analysis algorithm. The feature frequency and its doubling frequency are basically extracted.However,some feature information is neglected because of threshold se-lection and re-sampling process.Eurthermore,the difference frequency 10 Hz still exists in the fused signal.Consequently.It is concluded that compared with conventional algorithms,the proposed MUWDE algorithm has advantages in both disturbance restraining and feature information extracting.

    4 ExPerimental Results

    In order to verify the effectiveness and utility of the proposed algorithm,MUWDE is applied in the fusion of loose slipper(Eig.9)multi-channel vibration signals of hydraulic pump.The tested hydraulic pump is SY-10MCY14-1EL with 7 pistons.The driving dynamo is Y132M-4,with the settled speed of 1 480 r/min and the period of 0.041 s.Installation of the vibration sensors is shown in Eig.10.Sampling frequency is 12 000 Hz.The sampled vibration signals are saved in computer by DH-5920 dynamic signal testing and analyzing system.

    Eig.9 Eault of loose slipper

    Eig.10 Installation of vibration sensors on hydraulic pump

    Signals observed by Sensor1—Sensor3 are respectively saved as X1—X3.Eig.11 shows the time domain.And frequency domain is shown in Eig.11.

    Eig.11 Erequency domain of multi-channel signals of hydraulic pump

    In Eig.11,it is clear that all the three signals have modulation phenomena and the amplitude in 517 Hz is relatively strong.Since the real speed of pump is 1 480 r/min,the rotating bearing frequency is 1 480/60=24.6 Hz.As there are 7 pistons,the inherent striking frequency is 24.6× 7=172.2 Hz.Therefore,this frequency is supposed to exist in each vibration signal.Eor single loose slipper fault,its feature frequency should be equal or approximate to the rotating bearing frequency 24.6 Hz.However,both the inherent frequency and the fault feature frequency submerge by noises.Therefore,the multi-channel signals will be fused by MUWDE.

    According to Eqs.(10—12),we can get that Nmin=2 and Nmax=10.Eor each N=2,3,…,10,there will be a corresponding region for L.As fs=12 000 Hz and f0=24.6 Hz,the variousranges of L can be achieved by Eqs.(13—14),which are listed in Table 2.

    Table 2 Various combinations of N and L

    According to Eq.(15),there are 949 combinations of N and L.And each combination has its own EEE which is calculated by Eq.(16).Consequently,there are 9 curves consist of 949 groups of EEE.Collection of the minimum EEE of each curve is{0.052 0,0.041 6,0.041 4,0.035 3,0.037 2,0.040 9,0.041 3,0.042 3,0.042 9}. Eig.12 is the amplification of EEE ranging from 0.03 to 0.05.

    Eig.12 Amplification of the curves in the region of EEE ranging from 0.03 to 0.05

    In Eig.12,it is clear that the global minimum EEE(EEE=0.035 3)is achieved when N=5 and L=4.Therefore,with the selected parameters,the multi-channel vibration signals are fused by the proposed MUWDE.g0=[0 0 0 0]and N=5.Signals X1,X2and X3are decomposed by multi-scaled difference operator according to Eqs.(4—6).According to Eq.(7),EEE of each decomposition layer of X1is{0.028 7,0.022 6,0.030 9,0.035 6,0.026 8}.Similarly,EEE of X2is{0.028 1,0.047 7,0.054 4,0.031 8,0.031 5}and EEE of X3is{0.027 3,0.038 9,0.061 2,0.032 1,0.030 4}.Based on fusion rule 1,the selected approximate signals are{1st of X1,2nd of X2,3rd of X3,4th of X1,5th of X2}.Values of EEE are{0.028 7,0.047 7,0.061 2,0.035 6,0.031 5}.Then fusion weights are achieved by Eq.(8),which are k1=0.140,k2=0.233,k3=0.299,k4=0.174,k5=0.154. Consequently,the signal is reconstructed based on the fusion of 5 approximate signals by Eq.(9)

    where yij(i=1,2,3 and j=1,2,…,5)denotes the approximate signal in the jthdecomposition layer of yi;and yfinalthe reconstructed signal.Eusion result is shown in Eig.13.

    Eig.13 Experimental fusion results of multi-channel vibration signals based on MUWDE when N=5 and L=4

    In Eig.13,feature frequency 24.28 Hz(24.28 Hz≈24.46 Hz)and its doubling frequency are effectively extracted,as well as the inherent vibration frequency 170.9 Hz(170.9 Hz≈172.2 Hz).Noises and modulation phenomenon are better processed.To further indicate the advantages of MUWDE,we take the multi-scaled wavelet fusion(MWE)algorithm in Ref.[9]for comparison.Eusion result is shown in Eig.14.

    Eig.14 Experimental results of multi-channel vibration signals processed by MWE algorithm

    By comparison between Eig.14 and Eig.13,we can see that MWE is only able to extract in-herent frequency information.However,parts of feature information are neglected due to threshold selection and re-sampling.Therefore,feature information is not effectively utilized.Eurthermore,we can respectively calculate values of EEE and EEE after MWE fusion.Results are EEE= 0.051 5 and EEE=0.049 8.It is obvious that EEE of MUWDE is lower than that of MWE(0.035 3<0.051 5).And EEE of MUWDE is higher than that of MWE(0.069 4>0.049 8). Conclusions can be drawn that compared with MWE,feature information are effectively utilized by MUWDE so that the extracted feature information is more obvious.

    5 Conclusions

    An effective fusion algorithm for vibration signals of hydraulic pump based on MUWDE is proposed in this article.On one hand,approximate signals in various decomposition layers are selected by feature energy factor and fused for reconstruction.On the other hand,method for parameters optimal selection of MUWDE is demonstrated in detail.Simulation and experimental results show that compared with conventional algorithms,MUWDE can achieve better fusion performance and effectively extract feature information of hydraulic pump.It will help condition based maintenance.

    Acknowledgement

    This work was supported by the National Natural Science Eoundation of China(No.51275524).

    [1] Wang H W,Wu H Q.Reliability evaluation model based on data fusion for aircraft engines[J].Transactions of Nanjing University of Aeronautics&Astronautics,2012,29(4):318-324.

    [2] Ciuonzo D,Papa G,Romano G,et al.One-bit decentralized detection with a rao test for multi-sensor fusion[J].IEEE Signal Processing Letters,2013,20(9):861-864.

    [3] Safizadeh M S,Latifi S K.Using multi-sensor data fusion for vibration fault diagnosis of rolling element bearings by accelerometer and load cell[J].Information fusion,2014,18:1-8.

    [4] Angelov P,Yager R.Density-based averaging—A new operator for data fusion[J].Information Science,2013,222:164-174.

    [5] Wei C M,Blum R S.Theoretical analysis of correlation-based quality measures for weighted averaging image fusion[J].Information Eusion,2010,11:301-309.

    [6] Rodger J A.Toward reducing failure risk in an integrated vehicle health maintenance system:A fuzzy multi-sensor data fusion Kalman filter approach for IVH MS[J].Expert Systems with Applications,2012,39(10):9821-9835.

    [7] Yang W B,Li S Y.A switch-mode information fusion filter based on ISRUKE for autonomous navigation of spacecraft[J].Information Eusion,2014,18:33-42.

    [8] Ren Y E,Ke X Z.Multi-wavelet-basis multi-scale multi-sensor data fusion[J].Transducer and Microsystem Technologies,2010,29(9):77-79.

    [9] Lv X Y,Gu X H.Eiltering algorithm for multi-microphones based on wavelet multi-scale information fusion[J].Chinese Journal of Science Instrument,2012,33(4):789-794.(in Chinese)

    [10]Zhang J E,Smith J S,Wu Q H.Morphological undecimated wavelet decomposition for fault location on power transmission lines[J].IEEE Transactions on Circuits and Systems,2006,53(6):1395-1402.

    [11]Lin Y,Yang Y D,Liu J.Eeature extraction methods of vibration signal in automobile main reducer based on morphological un-decimated wavelet[J].Journal of Agricultural Engineering,2010,41(2):209-214.(in Chinese)

    [12]Huang B E,Shen L,Zhou X J,et al.Eault feature extraction of rolling element bearing based on morphological undecimated wavelet decomposition[J]. Journal of Agricultural Engineering,2010,41(2):204-207.(in Chinese)

    [13]Sun Z H,Lv W Q.Signal analysis of rolling mill vibration based on morphological un-decimated wavelets and S-transform[J].Journal of University of Science and Technology Beijing,2013,35(3):366-370.(in Chinese)

    [14]Higger M,Akcakaya M,Erdogmus D.A robust fusion algorithm for sensor failure[J].IEEE Signal Processing Letters,2013,20(8):755-758.

    [15]Zhang L J,Yang J H,Xu J W,et al.Morphological undecimated wavelet and its application to feature extraction of impulsive signal[J].Journal of Vibration and Shock,2007,26(10):56-59.(in Chinese)

    [16]Morenilla A J,Carmona R M,Romero J L S.Mathematical morphology for design and manufacturing[J].Mathematical and Computer Modeling,2011,54(7):1753-1759.

    [17]Cousty J,Najman L,Dias E,et al.Morphological filtering on graphs[J].Computer Vision and Image Understanding,2013,117(4):370-385.

    [18]Li B,Zhang P L,Wang Z J.Gear fault detection using multi-scale morphological filters[J].Measurement,2011,44(10):2078-2089.

    [19]He W,Jiang Z N,Qin Q.A joint adaptive wavelet filter and morphological signal processing method for weak mechanical impulse extraction[J].Journal of Mechanical Science and Technology,2011,24(8):1709-1716.

    [20]Hao R J,Chu E L.Morphological un-decimated wavelet decomposition for fault diagnostics of rolling element bearings[J].Journal of Sound and Vibration,2009,4(1):1164-1177.

    [21]Wang B,Li H R,Xu B H.Motor bearing forecast feature extracting and degradation status identification based on multi-scale morphological decomposition spectral entropy[J].Journal of Vibration and Shock,2013,32(22):124-128.(in Chinese)

    (Executive editor:Zhang Bei)

    TH165 Document code:A Article ID:1005-1120(2015)03-0268-11

    *CorresPonding author:Li Hongru,Professor,E-mail:lihr168@sohu.com.

    How to cite this article:Sun Jian,Li Hongru,Wang Weiguo,et al.Morphological undecimated wavelet decomposition fusion algorithm and its application on fault feature extraction of hydraulic pump[J].Trans.Nanjing U.Aero.Astro.,2015,32(3):268-278.

    http://dx.doi.org/10.16356/j.1005-1120.2015.03.268

    (Received 4 January 2015;revised 5 April 2015;accepted 25 April 2015)

    猜你喜歡
    孫健
    孫 健 作品欣賞
    藝術(shù)品(2019年7期)2019-08-23 12:13:52
    漫畫天地
    清風(fēng)(2019年1期)2019-04-20 11:42:44
    《化妝標(biāo)準(zhǔn)》《貼心小棉襖》
    三月三(2017年8期)2017-09-02 21:38:19
    量臀定碼
    三月三(2017年6期)2017-07-01 08:13:42
    量臀定碼
    三月三(2017年6期)2017-07-01 07:31:42
    心痛不如行動(dòng)
    心理與健康(2017年3期)2017-05-30 10:48:04
    微信時(shí)代
    三月三(2016年12期)2016-12-27 18:09:28
    也有會(huì)吹的
    三月三(2016年12期)2016-12-27 18:09:10
    Modelling hydrodynamic processes in tidal stream energy extraction*
    漫畫四幅
    精品99又大又爽又粗少妇毛片 | 精品日产1卡2卡| 国产精品电影一区二区三区| 国产一区二区三区视频了| 搡老岳熟女国产| 国产一区二区在线av高清观看| 美女cb高潮喷水在线观看 | 十八禁网站免费在线| 欧美xxxx黑人xx丫x性爽| 两人在一起打扑克的视频| 国产av麻豆久久久久久久| 国产极品精品免费视频能看的| 久久精品夜夜夜夜夜久久蜜豆| 国产精品自产拍在线观看55亚洲| 夜夜爽天天搞| 日本在线视频免费播放| 熟女电影av网| 色吧在线观看| 又大又爽又粗| 国产成人av激情在线播放| 又紧又爽又黄一区二区| 美女高潮喷水抽搐中文字幕| 亚洲色图 男人天堂 中文字幕| 国产人伦9x9x在线观看| 成年女人看的毛片在线观看| 99在线人妻在线中文字幕| 亚洲五月婷婷丁香| 国产高清三级在线| 亚洲成人久久性| 五月玫瑰六月丁香| 欧美日韩福利视频一区二区| 午夜日韩欧美国产| 无遮挡黄片免费观看| 舔av片在线| 美女 人体艺术 gogo| 99热精品在线国产| 亚洲av成人不卡在线观看播放网| 极品教师在线免费播放| 婷婷丁香在线五月| av天堂中文字幕网| 男女床上黄色一级片免费看| 欧美日韩福利视频一区二区| 亚洲成av人片在线播放无| 国产视频一区二区在线看| 欧美日韩综合久久久久久 | 亚洲人成伊人成综合网2020| 欧美成人性av电影在线观看| 色精品久久人妻99蜜桃| 成人性生交大片免费视频hd| 欧美色欧美亚洲另类二区| 中文字幕久久专区| 国产精品亚洲一级av第二区| 国产一区二区在线观看日韩 | 国产精品一区二区三区四区久久| 午夜福利在线观看吧| 欧美另类亚洲清纯唯美| 熟妇人妻久久中文字幕3abv| 亚洲av成人不卡在线观看播放网| 国产三级黄色录像| 精品一区二区三区视频在线 | 成人永久免费在线观看视频| 日韩欧美国产一区二区入口| 国产一区二区激情短视频| 真人一进一出gif抽搐免费| 中文字幕最新亚洲高清| 又黄又粗又硬又大视频| 一级黄色大片毛片| 国产真人三级小视频在线观看| 免费大片18禁| 久久久久性生活片| 亚洲中文日韩欧美视频| 一进一出抽搐gif免费好疼| 亚洲精品美女久久久久99蜜臀| 丝袜人妻中文字幕| 黄色丝袜av网址大全| 激情在线观看视频在线高清| www日本黄色视频网| 黄频高清免费视频| 黄色片一级片一级黄色片| 欧美极品一区二区三区四区| 亚洲av成人精品一区久久| 精品国产亚洲在线| 99精品欧美一区二区三区四区| 51午夜福利影视在线观看| 后天国语完整版免费观看| 99久久成人亚洲精品观看| 免费观看精品视频网站| 亚洲自拍偷在线| 在线观看美女被高潮喷水网站 | 日日干狠狠操夜夜爽| 精品午夜福利视频在线观看一区| 国产毛片a区久久久久| 亚洲美女视频黄频| 午夜精品在线福利| 波多野结衣巨乳人妻| 婷婷精品国产亚洲av在线| 国产精品自产拍在线观看55亚洲| 午夜成年电影在线免费观看| 99久久久亚洲精品蜜臀av| 久久久水蜜桃国产精品网| 成人av在线播放网站| 一区二区三区国产精品乱码| www日本黄色视频网| 中亚洲国语对白在线视频| 日韩欧美国产在线观看| 久久精品人妻少妇| 亚洲黑人精品在线| 日日摸夜夜添夜夜添小说| av中文乱码字幕在线| 午夜亚洲福利在线播放| 最新美女视频免费是黄的| 动漫黄色视频在线观看| 久久久国产精品麻豆| e午夜精品久久久久久久| 免费看光身美女| 99久久精品热视频| 99久久国产精品久久久| 亚洲午夜精品一区,二区,三区| 久久精品国产亚洲av香蕉五月| 国产亚洲av嫩草精品影院| 看免费av毛片| 国产野战对白在线观看| 久久久久久久精品吃奶| 制服人妻中文乱码| 悠悠久久av| 噜噜噜噜噜久久久久久91| 久久草成人影院| 国产乱人伦免费视频| 黄片小视频在线播放| 99在线视频只有这里精品首页| 亚洲中文日韩欧美视频| 久久久成人免费电影| 亚洲成av人片免费观看| 欧美高清成人免费视频www| 在线视频色国产色| 国产亚洲精品久久久com| 可以在线观看毛片的网站| 国产精品影院久久| 这个男人来自地球电影免费观看| 中文字幕高清在线视频| 国内揄拍国产精品人妻在线| 热99re8久久精品国产| 又黄又粗又硬又大视频| 欧美黑人欧美精品刺激| 国产亚洲欧美98| 手机成人av网站| 岛国视频午夜一区免费看| 一进一出抽搐动态| 亚洲av五月六月丁香网| 亚洲精品粉嫩美女一区| 国产精品亚洲av一区麻豆| 国产精品一区二区精品视频观看| 亚洲精品美女久久久久99蜜臀| aaaaa片日本免费| 757午夜福利合集在线观看| 91麻豆精品激情在线观看国产| 国产真实乱freesex| 欧美乱色亚洲激情| 长腿黑丝高跟| 欧美另类亚洲清纯唯美| 成人国产一区最新在线观看| 久久久久国内视频| 国产成+人综合+亚洲专区| 色综合亚洲欧美另类图片| 熟女人妻精品中文字幕| 国产精品98久久久久久宅男小说| 亚洲精品粉嫩美女一区| 欧美国产日韩亚洲一区| 女人高潮潮喷娇喘18禁视频| 午夜免费激情av| 中文字幕高清在线视频| 久久婷婷人人爽人人干人人爱| 欧美色视频一区免费| 久久久久免费精品人妻一区二区| 美女黄网站色视频| 国产v大片淫在线免费观看| 九九在线视频观看精品| 18禁美女被吸乳视频| 国产高清三级在线| 久久国产乱子伦精品免费另类| 午夜福利在线观看吧| 国产激情偷乱视频一区二区| 国产蜜桃级精品一区二区三区| 成人av在线播放网站| 国产精品乱码一区二三区的特点| 99热只有精品国产| 老司机午夜福利在线观看视频| 亚洲国产高清在线一区二区三| 一个人免费在线观看的高清视频| 成人18禁在线播放| 色在线成人网| 麻豆国产av国片精品| 在线观看舔阴道视频| 国产精品久久久人人做人人爽| 熟女人妻精品中文字幕| 99视频精品全部免费 在线 | 国产黄片美女视频| 叶爱在线成人免费视频播放| 99热精品在线国产| 亚洲国产看品久久| 桃色一区二区三区在线观看| 免费在线观看影片大全网站| 中亚洲国语对白在线视频| 99国产精品99久久久久| 精品乱码久久久久久99久播| 亚洲人成伊人成综合网2020| 麻豆久久精品国产亚洲av| 18禁美女被吸乳视频| 天堂网av新在线| 在线视频色国产色| 在线免费观看不下载黄p国产 | 欧美日韩一级在线毛片| 精品不卡国产一区二区三区| 亚洲欧美一区二区三区黑人| 中文字幕高清在线视频| www.999成人在线观看| 国产成人欧美在线观看| 大型黄色视频在线免费观看| 黄片大片在线免费观看| 久久精品亚洲精品国产色婷小说| 波多野结衣高清作品| 亚洲第一欧美日韩一区二区三区| 欧美+亚洲+日韩+国产| 中亚洲国语对白在线视频| 色综合婷婷激情| 国产蜜桃级精品一区二区三区| 国内精品久久久久久久电影| 国产高清视频在线播放一区| 成人无遮挡网站| 国产亚洲精品av在线| 久久精品国产清高在天天线| 最新在线观看一区二区三区| 伊人久久大香线蕉亚洲五| 色精品久久人妻99蜜桃| 欧美乱色亚洲激情| 久久这里只有精品中国| 国产真实乱freesex| www日本在线高清视频| 欧美av亚洲av综合av国产av| 欧美乱妇无乱码| 天天一区二区日本电影三级| 欧美性猛交╳xxx乱大交人| 久久久久久久久免费视频了| 婷婷六月久久综合丁香| 免费看日本二区| 身体一侧抽搐| 日韩高清综合在线| 亚洲av五月六月丁香网| 日本a在线网址| 国产av麻豆久久久久久久| 久久99热这里只有精品18| xxxwww97欧美| 国产乱人视频| 精品欧美国产一区二区三| 岛国在线观看网站| 国内毛片毛片毛片毛片毛片| 亚洲国产欧美网| 岛国视频午夜一区免费看| 成人三级做爰电影| 国产乱人视频| 嫁个100分男人电影在线观看| 亚洲va日本ⅴa欧美va伊人久久| 国产一区二区在线av高清观看| 99久久成人亚洲精品观看| 久久国产乱子伦精品免费另类| 国产真实乱freesex| 小说图片视频综合网站| cao死你这个sao货| www日本在线高清视频| 欧美性猛交黑人性爽| 观看免费一级毛片| 国产亚洲精品久久久com| 亚洲午夜精品一区,二区,三区| 嫩草影院精品99| 天堂网av新在线| 久久婷婷人人爽人人干人人爱| 国产69精品久久久久777片 | 国内久久婷婷六月综合欲色啪| 2021天堂中文幕一二区在线观| 国产欧美日韩一区二区精品| 午夜两性在线视频| 又大又爽又粗| 麻豆成人av在线观看| 日韩欧美免费精品| 亚洲成av人片在线播放无| 一二三四社区在线视频社区8| 色尼玛亚洲综合影院| 午夜福利18| 国产真实乱freesex| 熟女电影av网| a在线观看视频网站| 亚洲va日本ⅴa欧美va伊人久久| 久久久久久久午夜电影| 99久久精品一区二区三区| 国产精品,欧美在线| 欧美极品一区二区三区四区| 国产一区在线观看成人免费| 狂野欧美激情性xxxx| 麻豆久久精品国产亚洲av| 日韩高清综合在线| 国产视频一区二区在线看| 久久久久久久精品吃奶| 欧美激情久久久久久爽电影| 日日干狠狠操夜夜爽| 久久久久免费精品人妻一区二区| 亚洲性夜色夜夜综合| 丁香六月欧美| 国产午夜精品久久久久久| 观看免费一级毛片| 97超视频在线观看视频| 欧美中文日本在线观看视频| 男人舔奶头视频| 亚洲欧美日韩高清专用| 精品国产超薄肉色丝袜足j| 欧美成人一区二区免费高清观看 | 午夜精品在线福利| 国内精品美女久久久久久| 性色avwww在线观看| 激情在线观看视频在线高清| 免费人成视频x8x8入口观看| 成人国产综合亚洲| 日韩av在线大香蕉| 精品福利观看| 老鸭窝网址在线观看| 超碰成人久久| 日本免费a在线| netflix在线观看网站| 久久香蕉国产精品| 一区二区三区国产精品乱码| 90打野战视频偷拍视频| 国产欧美日韩精品一区二区| 亚洲精品国产精品久久久不卡| 91老司机精品| 亚洲18禁久久av| 国产 一区 欧美 日韩| 热99在线观看视频| www.熟女人妻精品国产| 久久天堂一区二区三区四区| 国产高清videossex| avwww免费| 老司机午夜福利在线观看视频| 2021天堂中文幕一二区在线观| 三级毛片av免费| 亚洲 国产 在线| 18禁黄网站禁片午夜丰满| 麻豆久久精品国产亚洲av| 给我免费播放毛片高清在线观看| 身体一侧抽搐| 操出白浆在线播放| 欧美又色又爽又黄视频| 欧美黄色淫秽网站| 久久人人精品亚洲av| 99热6这里只有精品| 久久精品人妻少妇| 国产视频内射| 村上凉子中文字幕在线| 久久精品国产99精品国产亚洲性色| 2021天堂中文幕一二区在线观| 在线看三级毛片| 午夜福利在线在线| 成人鲁丝片一二三区免费| 午夜影院日韩av| 国产97色在线日韩免费| 97人妻精品一区二区三区麻豆| 美女高潮的动态| 日本与韩国留学比较| 欧美色欧美亚洲另类二区| 国产av在哪里看| 久久这里只有精品19| 999久久久国产精品视频| 国产精品美女特级片免费视频播放器 | 国产综合懂色| 久久久久久久久久黄片| 精品国产乱子伦一区二区三区| 国产伦精品一区二区三区视频9 | 国产成人精品久久二区二区免费| a级毛片a级免费在线| 日韩欧美在线乱码| 欧美日韩乱码在线| 色老头精品视频在线观看| 91久久精品国产一区二区成人 | 亚洲成人久久性| 亚洲人成电影免费在线| 久久久久亚洲av毛片大全| 可以在线观看的亚洲视频| 一进一出抽搐gif免费好疼| 国产成人aa在线观看| 日本免费一区二区三区高清不卡| 久久久久性生活片| 午夜影院日韩av| 久久精品国产99精品国产亚洲性色| 日韩欧美精品v在线| 亚洲第一欧美日韩一区二区三区| 久久亚洲精品不卡| 国产又黄又爽又无遮挡在线| 99在线人妻在线中文字幕| 国产1区2区3区精品| 欧美日韩中文字幕国产精品一区二区三区| 国产亚洲精品av在线| 免费大片18禁| 欧美日韩精品网址| 午夜福利视频1000在线观看| 999久久久精品免费观看国产| 1024手机看黄色片| 嫩草影视91久久| 两人在一起打扑克的视频| 在线视频色国产色| 欧美国产日韩亚洲一区| 波多野结衣高清作品| 丁香六月欧美| 在线观看美女被高潮喷水网站 | 色尼玛亚洲综合影院| 国产91精品成人一区二区三区| 九色国产91popny在线| 国产精品一区二区三区四区久久| 中文字幕高清在线视频| 性色av乱码一区二区三区2| 欧美激情久久久久久爽电影| 免费看美女性在线毛片视频| 成人欧美大片| 嫩草影院精品99| 欧美乱码精品一区二区三区| 久久久久久久午夜电影| 怎么达到女性高潮| 噜噜噜噜噜久久久久久91| 熟妇人妻久久中文字幕3abv| 色精品久久人妻99蜜桃| 我的老师免费观看完整版| 18禁黄网站禁片午夜丰满| 国产男靠女视频免费网站| 欧美一级毛片孕妇| 18禁国产床啪视频网站| 99久久精品一区二区三区| 黄色视频,在线免费观看| 国产精品乱码一区二三区的特点| 这个男人来自地球电影免费观看| 国产高清视频在线播放一区| 亚洲成人精品中文字幕电影| 中出人妻视频一区二区| 国产精品久久久人人做人人爽| 精品久久久久久久久久免费视频| 国产精品爽爽va在线观看网站| 国产精品亚洲美女久久久| 久久久久久久精品吃奶| 久久精品91蜜桃| 欧美日韩一级在线毛片| 亚洲国产精品sss在线观看| 蜜桃久久精品国产亚洲av| 丰满人妻一区二区三区视频av | 岛国在线观看网站| 男人的好看免费观看在线视频| 欧美一区二区精品小视频在线| 久久久久性生活片| 亚洲欧美日韩高清在线视频| 麻豆国产97在线/欧美| 韩国av一区二区三区四区| av视频在线观看入口| 婷婷精品国产亚洲av| 免费观看的影片在线观看| 日本免费一区二区三区高清不卡| 亚洲精品456在线播放app | 亚洲 欧美一区二区三区| 成在线人永久免费视频| 一个人看的www免费观看视频| 亚洲成av人片免费观看| 免费看a级黄色片| 成年免费大片在线观看| 日韩欧美一区二区三区在线观看| 免费在线观看日本一区| 午夜亚洲福利在线播放| 中文亚洲av片在线观看爽| 午夜精品久久久久久毛片777| 男女午夜视频在线观看| 高清在线国产一区| 琪琪午夜伦伦电影理论片6080| 欧美绝顶高潮抽搐喷水| 男插女下体视频免费在线播放| 好男人电影高清在线观看| 亚洲欧美精品综合久久99| 久久中文看片网| 天堂av国产一区二区熟女人妻| 国产黄色小视频在线观看| 一级黄色大片毛片| 99久久国产精品久久久| 久久亚洲精品不卡| 美女被艹到高潮喷水动态| 少妇裸体淫交视频免费看高清| 亚洲欧美精品综合一区二区三区| 成人av在线播放网站| 丁香欧美五月| 精品一区二区三区av网在线观看| 精品国产乱码久久久久久男人| 国语自产精品视频在线第100页| 男插女下体视频免费在线播放| 男人舔奶头视频| 欧美另类亚洲清纯唯美| 日韩三级视频一区二区三区| 欧美绝顶高潮抽搐喷水| 青草久久国产| 女人高潮潮喷娇喘18禁视频| 亚洲欧洲精品一区二区精品久久久| 亚洲中文av在线| 国产毛片a区久久久久| 精品欧美国产一区二区三| 亚洲成人精品中文字幕电影| 久久国产精品人妻蜜桃| 国产精品av久久久久免费| 免费av不卡在线播放| 在线观看免费午夜福利视频| 国产不卡一卡二| 高潮久久久久久久久久久不卡| 俄罗斯特黄特色一大片| 中文亚洲av片在线观看爽| x7x7x7水蜜桃| 在线永久观看黄色视频| 亚洲午夜理论影院| 这个男人来自地球电影免费观看| 国产精品影院久久| 亚洲天堂国产精品一区在线| 久久久精品欧美日韩精品| 久久国产精品影院| 亚洲欧美激情综合另类| 久久精品国产清高在天天线| 91九色精品人成在线观看| 99riav亚洲国产免费| 亚洲精品在线观看二区| 午夜久久久久精精品| 久久久精品欧美日韩精品| 麻豆成人av在线观看| 免费在线观看成人毛片| 国产私拍福利视频在线观看| 国产亚洲精品av在线| 99riav亚洲国产免费| 亚洲真实伦在线观看| 一本精品99久久精品77| 国产精品女同一区二区软件 | 国产69精品久久久久777片 | 久久精品亚洲精品国产色婷小说| 禁无遮挡网站| 国产精品久久久av美女十八| 一级毛片女人18水好多| 成人午夜高清在线视频| 国产高清视频在线播放一区| 国产成人啪精品午夜网站| 久久久久久久午夜电影| 变态另类成人亚洲欧美熟女| 久久久久性生活片| 在线免费观看的www视频| 一本综合久久免费| 精品欧美国产一区二区三| 巨乳人妻的诱惑在线观看| 中文字幕久久专区| 免费搜索国产男女视频| 嫩草影视91久久| 亚洲欧洲精品一区二区精品久久久| 天天添夜夜摸| 亚洲电影在线观看av| 老司机午夜福利在线观看视频| 国产欧美日韩精品亚洲av| 国内精品一区二区在线观看| 一级作爱视频免费观看| 男女做爰动态图高潮gif福利片| 亚洲av中文字字幕乱码综合| 欧美xxxx黑人xx丫x性爽| 免费看a级黄色片| 国产精品一区二区免费欧美| 两个人视频免费观看高清| 悠悠久久av| 午夜成年电影在线免费观看| 舔av片在线| 亚洲精品在线美女| 欧美成人免费av一区二区三区| 亚洲精品国产精品久久久不卡| av福利片在线观看| 夜夜夜夜夜久久久久| 99精品欧美一区二区三区四区| www日本在线高清视频| 女警被强在线播放| 亚洲激情在线av| 熟妇人妻久久中文字幕3abv| 99精品久久久久人妻精品| 国产精品av视频在线免费观看| 中文字幕av在线有码专区| 久久久久国内视频| 91av网一区二区| 看片在线看免费视频| 性色avwww在线观看| 搡老岳熟女国产| 波多野结衣巨乳人妻| 午夜免费激情av| 99re在线观看精品视频| 中文字幕高清在线视频| 免费观看的影片在线观看| 男人舔女人的私密视频| 亚洲第一欧美日韩一区二区三区| 精品国产美女av久久久久小说| 国产探花在线观看一区二区| 午夜福利免费观看在线| 国产毛片a区久久久久| 久久久色成人| 草草在线视频免费看| 国产伦一二天堂av在线观看| 别揉我奶头~嗯~啊~动态视频| 99久久精品热视频| 熟妇人妻久久中文字幕3abv| 可以在线观看的亚洲视频| 久久中文看片网| 婷婷亚洲欧美| 色视频www国产| 国产午夜精品久久久久久| 国产亚洲精品久久久久久毛片| 女人被狂操c到高潮| 久久精品综合一区二区三区| 亚洲中文字幕日韩| x7x7x7水蜜桃|