• <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*
    漫畫四幅
    免费在线观看黄色视频的| 亚洲av福利一区| 亚洲精品色激情综合| 国产成人av激情在线播放| 一本大道久久a久久精品| 男人舔女人的私密视频| 日韩,欧美,国产一区二区三区| 日本黄大片高清| 少妇人妻久久综合中文| 免费观看无遮挡的男女| 国产精品一二三区在线看| 男人添女人高潮全过程视频| 男女下面插进去视频免费观看 | 一级爰片在线观看| 美女脱内裤让男人舔精品视频| 如日韩欧美国产精品一区二区三区| 欧美国产精品一级二级三级| 9191精品国产免费久久| 亚洲精品一区蜜桃| 九九在线视频观看精品| 爱豆传媒免费全集在线观看| 国产亚洲精品久久久com| 免费黄网站久久成人精品| 精品少妇黑人巨大在线播放| 午夜av观看不卡| 欧美精品av麻豆av| 99热这里只有是精品在线观看| 80岁老熟妇乱子伦牲交| freevideosex欧美| 亚洲欧洲精品一区二区精品久久久 | 在线观看人妻少妇| 中文字幕人妻丝袜制服| 亚洲国产欧美在线一区| 欧美日韩av久久| 国产日韩欧美在线精品| 边亲边吃奶的免费视频| 国产成人精品婷婷| 日韩 亚洲 欧美在线| 国产色婷婷99| 天堂中文最新版在线下载| 亚洲综合色惰| 2021少妇久久久久久久久久久| 亚洲少妇的诱惑av| 寂寞人妻少妇视频99o| 日韩大片免费观看网站| 久久国内精品自在自线图片| 人妻 亚洲 视频| 最近最新中文字幕大全免费视频 | kizo精华| 国产精品一国产av| 国产成人免费无遮挡视频| 欧美成人精品欧美一级黄| 伦精品一区二区三区| 国产精品国产三级专区第一集| 免费av中文字幕在线| 国产 精品1| 成年动漫av网址| 久久99一区二区三区| 黑人高潮一二区| 欧美国产精品va在线观看不卡| 91成人精品电影| 成人亚洲欧美一区二区av| 精品国产一区二区久久| 夜夜骑夜夜射夜夜干| 国产日韩欧美亚洲二区| 你懂的网址亚洲精品在线观看| 成人国产av品久久久| 国产精品无大码| 成人18禁高潮啪啪吃奶动态图| 97精品久久久久久久久久精品| 一区二区日韩欧美中文字幕 | 校园人妻丝袜中文字幕| 精品少妇久久久久久888优播| 亚洲成色77777| 美女大奶头黄色视频| 大香蕉久久成人网| 18禁裸乳无遮挡动漫免费视频| 国产亚洲最大av| 国产免费视频播放在线视频| av女优亚洲男人天堂| 蜜桃国产av成人99| 中文字幕免费在线视频6| 亚洲第一av免费看| 80岁老熟妇乱子伦牲交| 狠狠婷婷综合久久久久久88av| 成人国产麻豆网| 97在线视频观看| videossex国产| 韩国av在线不卡| 夜夜爽夜夜爽视频| 亚洲性久久影院| 天堂中文最新版在线下载| 国产一区二区在线观看日韩| 一级片'在线观看视频| www.色视频.com| 哪个播放器可以免费观看大片| 日韩不卡一区二区三区视频在线| 青春草亚洲视频在线观看| 又黄又粗又硬又大视频| 久久久久久人妻| 尾随美女入室| 国产女主播在线喷水免费视频网站| 一区二区av电影网| 精品一区在线观看国产| 亚洲国产精品一区三区| 大片电影免费在线观看免费| 草草在线视频免费看| 久久99精品国语久久久| 欧美精品人与动牲交sv欧美| 亚洲av福利一区| av卡一久久| 成人二区视频| 日本爱情动作片www.在线观看| 天天操日日干夜夜撸| 在线观看人妻少妇| 18禁在线无遮挡免费观看视频| 香蕉国产在线看| 国产av一区二区精品久久| 亚洲精品国产av成人精品| 国产精品成人在线| 欧美成人午夜精品| av网站免费在线观看视频| tube8黄色片| 精品少妇内射三级| 日本色播在线视频| 国产精品人妻久久久久久| 亚洲国产精品成人久久小说| 亚洲综合色惰| 久久精品人人爽人人爽视色| 亚洲成人一二三区av| 欧美少妇被猛烈插入视频| 你懂的网址亚洲精品在线观看| 亚洲美女搞黄在线观看| 精品少妇久久久久久888优播| 久久99精品国语久久久| 赤兔流量卡办理| 国产片特级美女逼逼视频| 久久久久久伊人网av| 亚洲中文av在线| 黄色 视频免费看| 久久影院123| 亚洲精品乱码久久久久久按摩| 亚洲三级黄色毛片| 春色校园在线视频观看| 日本wwww免费看| 欧美激情极品国产一区二区三区 | 成年动漫av网址| 最后的刺客免费高清国语| 天堂中文最新版在线下载| 91成人精品电影| 国语对白做爰xxxⅹ性视频网站| 成人国产麻豆网| 成人亚洲精品一区在线观看| 在线观看www视频免费| 亚洲精品国产色婷婷电影| 人体艺术视频欧美日本| 夫妻性生交免费视频一级片| 内地一区二区视频在线| 国产av国产精品国产| 亚洲天堂av无毛| 国产探花极品一区二区| 啦啦啦啦在线视频资源| 9色porny在线观看| 亚洲精品自拍成人| 精品人妻在线不人妻| 最近的中文字幕免费完整| 精品一区在线观看国产| 国产极品天堂在线| 日韩成人伦理影院| 久久精品夜色国产| 久久久精品94久久精品| 日韩av不卡免费在线播放| 亚洲精品视频女| 欧美精品国产亚洲| 我要看黄色一级片免费的| 久久午夜综合久久蜜桃| 欧美日韩成人在线一区二区| 高清视频免费观看一区二区| 欧美3d第一页| 国产麻豆69| 久久精品aⅴ一区二区三区四区 | 亚洲精品日韩在线中文字幕| 国产成人午夜福利电影在线观看| 97精品久久久久久久久久精品| 丰满迷人的少妇在线观看| 亚洲天堂av无毛| 90打野战视频偷拍视频| videos熟女内射| 国产精品嫩草影院av在线观看| 国产亚洲欧美精品永久| 国产永久视频网站| 午夜免费男女啪啪视频观看| 少妇熟女欧美另类| 啦啦啦中文免费视频观看日本| 久久精品国产自在天天线| 看十八女毛片水多多多| 97在线人人人人妻| 欧美日韩av久久| 国产午夜精品一二区理论片| 视频在线观看一区二区三区| 国产熟女午夜一区二区三区| 免费观看a级毛片全部| 最新中文字幕久久久久| 久久久久久伊人网av| 精品久久国产蜜桃| 91精品伊人久久大香线蕉| 国产男人的电影天堂91| 亚洲精品久久久久久婷婷小说| 捣出白浆h1v1| 国产免费又黄又爽又色| 欧美成人午夜精品| 好男人视频免费观看在线| 亚洲国产欧美在线一区| 在线观看免费日韩欧美大片| 大香蕉久久网| 久久精品熟女亚洲av麻豆精品| 妹子高潮喷水视频| 日韩欧美一区视频在线观看| 国产成人精品久久久久久| 日韩精品有码人妻一区| av免费在线看不卡| 国产一区亚洲一区在线观看| 亚洲欧洲日产国产| xxxhd国产人妻xxx| 99九九在线精品视频| 精品99又大又爽又粗少妇毛片| 丝袜脚勾引网站| 国产成人a∨麻豆精品| 91久久精品国产一区二区三区| 国产熟女午夜一区二区三区| 午夜福利网站1000一区二区三区| 欧美激情国产日韩精品一区| 涩涩av久久男人的天堂| 国产精品国产三级国产专区5o| 亚洲成人手机| 欧美国产精品一级二级三级| 搡女人真爽免费视频火全软件| 边亲边吃奶的免费视频| 国产成人精品福利久久| 日韩欧美一区视频在线观看| 老熟女久久久| 亚洲国产日韩一区二区| 伊人久久国产一区二区| 青春草国产在线视频| 亚洲国产最新在线播放| 国产视频首页在线观看| 精品一区二区免费观看| 在线观看免费日韩欧美大片| 免费观看av网站的网址| 国产午夜精品一二区理论片| 国产在线一区二区三区精| 青春草国产在线视频| 国产亚洲精品久久久com| 免费观看性生交大片5| 高清av免费在线| 日韩精品有码人妻一区| 日产精品乱码卡一卡2卡三| 欧美日韩综合久久久久久| 免费看不卡的av| 少妇精品久久久久久久| 久久久久久久大尺度免费视频| 欧美精品亚洲一区二区| 亚洲成人手机| 午夜免费鲁丝| 久久久国产一区二区| 视频区图区小说| 18在线观看网站| 全区人妻精品视频| 黄色配什么色好看| 成人午夜精彩视频在线观看| 国产色爽女视频免费观看| 妹子高潮喷水视频| 97人妻天天添夜夜摸| 欧美日韩国产mv在线观看视频| 亚洲情色 制服丝袜| 国国产精品蜜臀av免费| 亚洲av.av天堂| 国产熟女午夜一区二区三区| 丝袜人妻中文字幕| 永久网站在线| 黑丝袜美女国产一区| 欧美激情极品国产一区二区三区 | 国产一区二区在线观看日韩| 日产精品乱码卡一卡2卡三| 国产成人精品无人区| 熟女电影av网| 90打野战视频偷拍视频| 亚洲国产日韩一区二区| 久久亚洲国产成人精品v| 久久婷婷青草| 亚洲三级黄色毛片| tube8黄色片| 91成人精品电影| 韩国精品一区二区三区 | 亚洲第一av免费看| 99久久精品国产国产毛片| 好男人视频免费观看在线| 免费观看无遮挡的男女| xxxhd国产人妻xxx| 人人澡人人妻人| 亚洲av电影在线观看一区二区三区| 有码 亚洲区| 国产男女内射视频| 少妇的丰满在线观看| 七月丁香在线播放| 最近的中文字幕免费完整| 久久久久网色| 欧美 亚洲 国产 日韩一| 一本久久精品| 高清欧美精品videossex| 亚洲国产色片| 在现免费观看毛片| 国产精品嫩草影院av在线观看| 中文字幕精品免费在线观看视频 | 秋霞伦理黄片| 日韩欧美一区视频在线观看| 午夜福利乱码中文字幕| 欧美97在线视频| 国产成人精品一,二区| 在线观看人妻少妇| 国产精品久久久av美女十八| 国产精品麻豆人妻色哟哟久久| 一个人免费看片子| 制服丝袜香蕉在线| 久久午夜综合久久蜜桃| 亚洲精品日本国产第一区| 丝瓜视频免费看黄片| 深夜精品福利| 亚洲精品视频女| 国产男女内射视频| 国产老妇伦熟女老妇高清| 欧美精品亚洲一区二区| 国产精品国产av在线观看| 国产精品人妻久久久久久| 深夜精品福利| 久久久久久久久久久免费av| 大陆偷拍与自拍| 欧美xxⅹ黑人| 1024视频免费在线观看| 大香蕉久久成人网| 欧美bdsm另类| 亚洲欧洲国产日韩| 亚洲成国产人片在线观看| 国产国拍精品亚洲av在线观看| 免费av中文字幕在线| 99热6这里只有精品| 亚洲精品美女久久久久99蜜臀 | 男女啪啪激烈高潮av片| 欧美精品国产亚洲| 国内精品宾馆在线| 十八禁高潮呻吟视频| 伦理电影大哥的女人| 亚洲第一区二区三区不卡| 午夜福利网站1000一区二区三区| 亚洲精品乱码久久久久久按摩| 两性夫妻黄色片 | 欧美国产精品一级二级三级| 国产综合精华液| 久久韩国三级中文字幕| 亚洲综合精品二区| 国产片特级美女逼逼视频| 亚洲人与动物交配视频| 免费看av在线观看网站| 免费黄网站久久成人精品| 一个人免费看片子| 久久久欧美国产精品| 女人久久www免费人成看片| 精品熟女少妇av免费看| 亚洲一区二区三区欧美精品| 婷婷色综合大香蕉| 爱豆传媒免费全集在线观看| 激情五月婷婷亚洲| 久久97久久精品| av电影中文网址| 中文乱码字字幕精品一区二区三区| 黄色怎么调成土黄色| 人人妻人人添人人爽欧美一区卜| 色哟哟·www| 黄色怎么调成土黄色| 日韩免费高清中文字幕av| 最近的中文字幕免费完整| 亚洲精品国产av蜜桃| 亚洲三级黄色毛片| 久久精品国产亚洲av涩爱| 亚洲第一区二区三区不卡| av在线观看视频网站免费| 男人添女人高潮全过程视频| 内地一区二区视频在线| 久久婷婷青草| 岛国毛片在线播放| 交换朋友夫妻互换小说| 国产黄频视频在线观看| 七月丁香在线播放| 欧美性感艳星| 最近手机中文字幕大全| 国产精品一二三区在线看| 亚洲国产精品国产精品| 亚洲精品一二三| 老司机影院毛片| 最新的欧美精品一区二区| 哪个播放器可以免费观看大片| 日韩中文字幕视频在线看片| 亚洲精品国产色婷婷电影| 99九九在线精品视频| a级毛色黄片| 国产日韩欧美在线精品| 亚洲精品美女久久久久99蜜臀 | 久久久久久久亚洲中文字幕| 另类亚洲欧美激情| 免费高清在线观看日韩| 免费人妻精品一区二区三区视频| 精品国产一区二区久久| 精品一区在线观看国产| 亚洲高清免费不卡视频| 国产精品国产三级国产专区5o| 精品亚洲成国产av| 久久久久精品人妻al黑| 天天操日日干夜夜撸| 亚洲欧美成人综合另类久久久| 一级a做视频免费观看| 国产精品久久久久久av不卡| 麻豆乱淫一区二区| 国产高清三级在线| 国产在线视频一区二区| 成人毛片60女人毛片免费| 中文字幕人妻熟女乱码| 男女高潮啪啪啪动态图| 春色校园在线视频观看| 18禁动态无遮挡网站| 午夜免费男女啪啪视频观看| 日本免费在线观看一区| 免费在线观看黄色视频的| 亚洲av电影在线观看一区二区三区| 大香蕉97超碰在线| 久久精品国产a三级三级三级| 99国产精品免费福利视频| 永久免费av网站大全| 女人精品久久久久毛片| av女优亚洲男人天堂| 极品少妇高潮喷水抽搐| 在线精品无人区一区二区三| 日韩人妻精品一区2区三区| 七月丁香在线播放| 一二三四中文在线观看免费高清| 韩国高清视频一区二区三区| 久久婷婷青草| 日本黄大片高清| 色婷婷av一区二区三区视频| 精品少妇久久久久久888优播| 在线免费观看不下载黄p国产| 男人添女人高潮全过程视频| 一级毛片黄色毛片免费观看视频| 最新中文字幕久久久久| 久久久久久久精品精品| 丰满饥渴人妻一区二区三| 国产av一区二区精品久久| 满18在线观看网站| 一本色道久久久久久精品综合| 一二三四中文在线观看免费高清| 啦啦啦啦在线视频资源| 97超碰精品成人国产| 男女边摸边吃奶| 夫妻午夜视频| 国产精品一区www在线观看| 精品国产一区二区三区四区第35| 亚洲第一区二区三区不卡| 国产成人免费观看mmmm| av国产久精品久网站免费入址| 人人妻人人添人人爽欧美一区卜| 免费黄网站久久成人精品| 精品福利永久在线观看| 好男人视频免费观看在线| 亚洲精品一区蜜桃| 人体艺术视频欧美日本| av免费在线看不卡| 又黄又粗又硬又大视频| 亚洲第一av免费看| 国产免费一级a男人的天堂| 一边亲一边摸免费视频| 黄色视频在线播放观看不卡| 成年女人在线观看亚洲视频| 在线看a的网站| 一级a做视频免费观看| 亚洲欧美一区二区三区国产| 80岁老熟妇乱子伦牲交| 国产黄色免费在线视频| 亚洲欧洲精品一区二区精品久久久 | 男女高潮啪啪啪动态图| 亚洲精品久久成人aⅴ小说| 人人妻人人澡人人看| 99热6这里只有精品| 人体艺术视频欧美日本| 美女xxoo啪啪120秒动态图| 少妇精品久久久久久久| 亚洲色图 男人天堂 中文字幕 | av黄色大香蕉| 日本wwww免费看| 99久国产av精品国产电影| 在线观看免费视频网站a站| 丝袜人妻中文字幕| 久久久久久久精品精品| 高清毛片免费看| 天堂中文最新版在线下载| 精品久久国产蜜桃| 韩国高清视频一区二区三区| 亚洲国产日韩一区二区| 国产麻豆69| 女性被躁到高潮视频| 一区二区日韩欧美中文字幕 | 亚洲国产av影院在线观看| 欧美日韩国产mv在线观看视频| 国产日韩欧美在线精品| 亚洲av中文av极速乱| 久久综合国产亚洲精品| 久久久久久久大尺度免费视频| 大陆偷拍与自拍| 九色成人免费人妻av| 国产69精品久久久久777片| 在线天堂中文资源库| 99久久中文字幕三级久久日本| 亚洲四区av| 热99久久久久精品小说推荐| 五月伊人婷婷丁香| 性色avwww在线观看| 日韩精品免费视频一区二区三区 | 亚洲国产精品专区欧美| 国内精品宾馆在线| 宅男免费午夜| 国产精品秋霞免费鲁丝片| 在线 av 中文字幕| 国产伦理片在线播放av一区| 国产精品秋霞免费鲁丝片| 亚洲国产色片| 捣出白浆h1v1| 大码成人一级视频| 韩国高清视频一区二区三区| 亚洲成色77777| 色哟哟·www| 亚洲精品国产av蜜桃| 亚洲国产精品一区二区三区在线| 精品一区在线观看国产| 久久久久久人妻| 中文天堂在线官网| 高清不卡的av网站| 欧美97在线视频| 国产 一区精品| 亚洲精品国产av成人精品| 国产精品不卡视频一区二区| 国产麻豆69| 亚洲五月色婷婷综合| 成人亚洲精品一区在线观看| 国产精品免费大片| 男人舔女人的私密视频| 中文天堂在线官网| 国产成人免费无遮挡视频| 免费人成在线观看视频色| 波多野结衣一区麻豆| 男人爽女人下面视频在线观看| 亚洲国产日韩一区二区| 一边亲一边摸免费视频| 狠狠精品人妻久久久久久综合| 中文天堂在线官网| 国产精品久久久久久av不卡| 欧美变态另类bdsm刘玥| 亚洲国产av新网站| 中文字幕精品免费在线观看视频 | 少妇人妻久久综合中文| 麻豆乱淫一区二区| av.在线天堂| 亚洲 欧美一区二区三区| 国产淫语在线视频| 日本wwww免费看| 激情视频va一区二区三区| 肉色欧美久久久久久久蜜桃| 色视频在线一区二区三区| 成人综合一区亚洲| 99热国产这里只有精品6| 亚洲欧美日韩另类电影网站| 精品国产露脸久久av麻豆| 天美传媒精品一区二区| 久久精品aⅴ一区二区三区四区 | 精品人妻一区二区三区麻豆| 国产精品无大码| 九色亚洲精品在线播放| 国产精品久久久久久久电影| 免费观看性生交大片5| 亚洲人成网站在线观看播放| 18禁裸乳无遮挡动漫免费视频| 国产欧美日韩一区二区三区在线| 国产黄色免费在线视频| 国产精品熟女久久久久浪| 伦理电影大哥的女人| 9色porny在线观看| 午夜免费观看性视频| 一本色道久久久久久精品综合| 欧美日韩亚洲高清精品| 日本wwww免费看| 欧美97在线视频| 9色porny在线观看| 90打野战视频偷拍视频| 五月天丁香电影| 日产精品乱码卡一卡2卡三| 久久久久久人妻| 久久精品久久久久久噜噜老黄| 天堂俺去俺来也www色官网| 中文天堂在线官网| 欧美人与性动交α欧美软件 | 爱豆传媒免费全集在线观看| 亚洲图色成人| 久久久久精品久久久久真实原创| 国产无遮挡羞羞视频在线观看| 亚洲激情五月婷婷啪啪| av女优亚洲男人天堂| 永久网站在线| 又黄又粗又硬又大视频|