999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)模型優(yōu)化研究

2024-05-08 00:00:00辛超孫成田張效源孫凱孫凱進
粘接 2024年1期

摘 要:鍋爐燃燒系統(tǒng)的建模是實現(xiàn)鍋爐節(jié)能減排的基礎,為實現(xiàn)燃煤鍋爐節(jié)能減排,針對現(xiàn)有建模方法多步預測精度不足問題,采用經(jīng)驗模態(tài)分解算法把復雜的輸出信號轉(zhuǎn)化為多個具有周期性規(guī)律或趨勢相對平穩(wěn)的模態(tài)信號,降低數(shù)據(jù)復雜度,基于K近鄰聯(lián)合互信息法得出遲延時間,提出基于動態(tài)時間規(guī)整距離進行在線更新的最小二乘支持向量機算法,并進行鍋爐燃燒系統(tǒng)的建模,基于鍋爐實際運行數(shù)據(jù)的仿真結(jié)果表明,該方法可以有效提高模型的自適應能力和多步預測精度,為后續(xù)實現(xiàn)閉環(huán)燃燒優(yōu)化控制打下了良好基礎。

關(guān)鍵詞:鍋爐燃燒系統(tǒng);經(jīng)驗模態(tài)分解算法;K近鄰聯(lián)合互信息;最小二乘支持向量機

中圖分類號:

TP391.92;TQ534

文獻標志碼:

A文章編號:

1001-5922(2024)01-0117-04

Research on optimization of adaptive dynamic model of boiler combustion system based on EMD

XIN Chao1,SUN Chengtian1,ZHANG Xiaoyuan1,SUN Kai1,SUN Kaijin2

(1.CHN Energy Feixian Power Generation Co.,Ltd.,F(xiàn)eixian 276001,Shandong China;

2.Southeast University,Nanjing 210096,China)Abstract:Modelling of the boiler combustion system is the basis for energy saving and emission reduction of boilers,in order to realize the energy saving and emission reduction of coal-fired boilers and solve the problem of insufficient multi-step prediction accuracy of existing modeling methods,the empirical mode decomposition algorithm wasfirst used to transform the complex output signals into multiple modal signals with periodic rules or relatively stable trends,so as to reduce the data complexity.The delayed time was obtained based on the K-nearest neighbor joint mutual information method,a least-squares support vector machine algorithm based on the dynamic time regular distance for online update was proposed,the boiler combustion system was modeled,and the simulation results based on the actual operation data of the boiler showed that the proposed method could effectively improve the adaptability and multi-step prediction accuracy of the model,and lay a good foundation for the subsequent realization of closed-loop combustion optimization control.

Key words:boiler combustion system; empirical modaldecomposition algorithm; k-nearest neighbor joint mutual information; least squares support vector machine

燃燒優(yōu)化是實現(xiàn)燃煤鍋爐節(jié)能減排的重要手段,而建立準確的鍋爐燃燒系統(tǒng)模型是實現(xiàn)燃燒優(yōu)化的基礎。但由于鍋爐燃燒系統(tǒng)是一個復雜的動態(tài)系統(tǒng),其控制變量和影響因素眾多,且由于負荷和煤質(zhì)多變以及積灰等因素,導致鍋爐燃燒系統(tǒng)特性經(jīng)常發(fā)生變化,因此給鍋爐燃燒系統(tǒng)建模帶來很大挑戰(zhàn)[1]。目前圍繞鍋爐燃燒系統(tǒng)建模問題已經(jīng)開展了大量研究工作。最小二乘支持向量機算法(Least squares support vector machine,LSSVM)作為一種小樣本建模算法,具有較好的預測效果和較短的模型訓練時間,得到研究人員越來越多的關(guān)注。蘇濤[2]、李揚[3]和藍茂蔚[4]

使用LSSVM模型分別對煙氣含氧量和鍋爐效率進行建模[2-4]

;朱鈺森[5]、梁濤[6]和金秀章[7]

對LSSVM模型進行改進,提高了對鍋爐NOx排放量的預測精度和泛化性能[5-7]。經(jīng)驗模態(tài)分解算法(Empirical mode decomposition,EMD)是一種典型的信號分解算法,算法簡單,并且具有良好的分解效果,被廣泛應用。張雲(yún)欽等[8]

將EMD方法用于光伏功率預測[8],基于EMD的建模方法也在潮位測量[9]、水電機組振動[10]和狀態(tài)識別[11]等領(lǐng)域得到廣泛應用,并取得不錯的效果。基于此,研究提出一種基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)建模方法。

1 基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)建模方法

1.1 EMD-LSSVM預測建模

LSSVM是以二次損失函數(shù)為經(jīng)驗風險的支持向量機,將不等式約束以等式約束替代,通過線性方程組求解模型參數(shù),提高訓練效率[12]。

現(xiàn)有訓練樣本集D={(x1,y1),…,(xn,yn)}。其中xi為輸入特征且xi∈Rn,yi為輸出特征且yi∈R,i=1,2,…,n,n為訓練樣本個數(shù)。在上述基礎上構(gòu)造一決策函數(shù),如式(1)所示。其中,φ(xi)為非線性函數(shù),能夠?qū)⒃卣骺臻gRp的輸入特征xi映射至一個高維的特征空間Rh中,ω為權(quán)重向量,b為偏差。

f(x)=ωT·φ(xi)+b (1)

根據(jù)結(jié)構(gòu)風險最小化原理,選擇平方誤差為成本函數(shù),通過拉格朗日方法,推導優(yōu)化得到線性方程組,如式(2)所示:

L(ω,b,ξ,α)=12ωTω+12γ∑ni=1ξ2i-∑ni=1αi[ωTφ(xi)+b+ξi-yi] (2)

其中,γ是用于平衡模型的復雜度及準確性的懲罰系數(shù);ξi為松弛變量,表征部分樣本出現(xiàn)誤差的容許程度,ξ=[ξ1,ξ2,…,ξn]T;α=[α1,α2,…,αn]T為拉格朗日乘數(shù)向量。

根據(jù)最優(yōu)化KKT條件,可以得到線性方程組,求解方程組可以得出參數(shù)α和b,把其代入式(1),可以得到LSSVM模型:

y(x)=ωTφ(xi)+b=∑ni=1αiK(x,xi)+b (3)

基于LSSVM建模算法可以得出模型的預測值,如式(4)所示:

Y^=f(X)(4)

式中:Y^為模型輸出的預測值;f(·)為預測模型,即LSSVM模型;X為模型的輸入量。

EMD算法首先由Huang等提出,適用于處理非線性和非平穩(wěn)的時間序列數(shù)據(jù)[13]。EMD不需要預先設定基函數(shù),依據(jù)數(shù)據(jù)自身的時間尺度特征進行信號分解。其原理是通過不斷的迭代求解,將信號分解成一系列不同頻率尺度下的本征模態(tài)函數(shù)集合和一個余項,并且分解前后的序列長度不變[12]。本征模態(tài)函數(shù)表示原始信號在不同頻域下的特性,余項反映了原始信號的基本變化趨勢。

對于復雜度較大的數(shù)據(jù),LSSVM不能很好的學到數(shù)據(jù)的特征,預測性能不夠理想,而EMD算法能夠?qū)碗s的數(shù)據(jù)進行分解,降低其復雜度,因此本文將EMD與LSSVM結(jié)合起來,提高模型的學習效率,從復雜數(shù)據(jù)中準確學到數(shù)據(jù)特征,提高預測精度。

根據(jù)EMD算法分解得到的不同本征模態(tài)函數(shù),可以根據(jù)其特點分為2類進行建模:第一類是低頻信號這類變化平穩(wěn)的模態(tài)信號,模型的輸入選取為控制量、狀態(tài)量和輸出特征在歷史一段時間內(nèi)的實際值,如式(1)所示;第二類是高頻信號這類有一定周期性變化規(guī)律的模態(tài)信號,模型的輸入為歷史一段時間內(nèi)的輸出特征的實際值,如式(2)所示,這2類的輸出都是對應輸出特征下一時刻的預測值。

Y^s(t)=f(x(t),x(t-1),…,x(t-wx),Ys(t-1),Ys(t-2),…,Ys(t-wy))(1)

Y^q(t)=f(Yq(t-1),Yq(t-2),…,Yq(t-wy))(2)

式中:Y^s(t)和Y^q(t)分別為低頻和高頻信號的預測值;f(·)為預測模型,即LSSVM模型;x(·)為不同時刻的輸入量;wx為輸入量的階次;Ys(t)和Yq(t)分別為不同時刻低頻和高頻信號的真實值;wy為輸出量的階次。

將分解后的各信號的預測值進行相加,得到對原信號的預測值:

Y^=IMF^_1+IMF^_2+…+resid^ual(3)

式中:Y^為模型輸出的預測值;IM^F為模型輸出分解后各本征模態(tài)信號的預測值;resid^ual為余項的預測值。

EMD可以自適應的對復雜的數(shù)據(jù)進行分解,直至余項無法繼續(xù)分解才停止,如果人為的選取模態(tài)分解數(shù)量,當模態(tài)分解數(shù)量滿足要求時,EMD就會停止分解。因此,采用Pearson相關(guān)系數(shù)來確定余項和原始數(shù)據(jù)的相關(guān)程度,為了選取合適的模態(tài)數(shù)量,且余項能夠保留原始數(shù)據(jù)的趨勢,選取余項與原始數(shù)據(jù)的

當Pearson相關(guān)系數(shù)大于閾值P時的最大模態(tài)分解數(shù)量。

1.2 基于K近鄰聯(lián)合互信息的動態(tài)模型遲延分析方法

針對鍋爐燃燒優(yōu)化的目的,分別建立再熱汽溫模型、鍋爐效率模型和NOx排放量模型,為了能夠使燃燒優(yōu)化起到效果,所建立的燃燒系統(tǒng)模型應正確反映實際生產(chǎn)過程中的變量關(guān)系。通過對燃燒機理和燃燒過程的分析,選定模型的輸入特征為3個部分:給煤量、二次風門開度、燃盡風門開度和省煤器后氧量為控制變量、機組負荷為狀態(tài)變量和過去時刻輸出量[14-15]。

鍋爐燃燒系統(tǒng)是一個復雜的動態(tài)系統(tǒng),各輸入特征和輸出特征之間存在明顯的時滯,當輸入量發(fā)生變化,經(jīng)過一段時間之后,輸出量才會發(fā)生變化,因此需要引入輸入變量的階次,模型輸入中加入過去時刻的輸入量輸出量,建立鍋爐燃燒系統(tǒng)的動態(tài)模型。目前在計算鍋爐控制量的遲延時間時,大部分方法是單獨考慮每個控制量對應輸出特征的遲延[16-17],并沒有考慮到控制量之間的相互影響,但是鍋爐作為一個強耦合性的系統(tǒng),且運行過程中多個控制量同時動作,控制量之間會互相影響。

K近鄰聯(lián)合互信息方法是一種用于衡量多個時間序列相關(guān)程度的方法。假設X1,X2,…,Xm均為N長度的時序數(shù)據(jù),則這m組數(shù)據(jù)的K近鄰聯(lián)合互信息(JMI)為 [18]:

KNN_

JMI(X1,X2,…,Xm)=ψ(k)-1N∑Ni=1ψ[nX1(i)+1]+…+ψ[nXm(i)+1]+(m-1)ψ(N)(4)

式中:ψ(·)為雙Γ函數(shù);nx(i)是與樣本xi之間距離嚴格小于0.5ε(ni,K)的樣本數(shù)。互信息值越大,這m組數(shù)據(jù)之間的相關(guān)性越大。

對于鍋爐燃燒系統(tǒng),采用式(4)的聯(lián)合互信息能夠同時考慮到每個輸入變量對輸出變量的有用信息貢獻[19],以便準確分析遲延變量和輸出變量之間的相關(guān)性。

KNN_JMI(X1,X2,…,Xm;Y)=KNN_JMI(X1,X2,…,Xm,Y)-KNN_JMI(X1,X2,…,Xm)(5)

式中:X1,X2,…,Xm是各輸入變量;Y是輸出量。KNN_JMI(X1,X2,…,Xm;Y)是考慮到輸入變量之間互相影響的情況,輸出量Y與多個輸入變量X1,X2,…,Xm之間的聯(lián)合互信息;KNN_JMI(X1,X2,…,Xm,Y)是沒有考慮到輸入變量之間互相影響的情況,輸出量Y與多個輸入變量X1,X2,…,Xm之間的聯(lián)合互信息,KNN_JMI(X1,X2,…,Xm)是各輸入變量之間的聯(lián)合互信息。KNN_JMI(X1,X2,…,Xm,Y)越大,輸出量對應的所有的輸入變量包含更多的有用信息,KNN_JMI(X1,X2,…,Xm)越小,輸入變量之間的相關(guān)性也就越小,由式(5)可知,此時各輸入變量與輸出變量之間的互信息就越大。

由于不同的輸入變量Xj(t)與輸出Y之間的時滯不同,所以嵌入不同的時滯τj,如式(10)所示:

maxτi=τ*iKNN_DJMI(X^1,X^2,…,X^m;Y^)

s.t.τmin≤τi ≤τmax

X^i=Xi(t-τi-w+1),…,Xi(t-τi-1),Xi(t-τi)

Y^=Y(t-w+1),…,Y(t-1),Y(t)

i∈[1,m](6)

式(6)為一個受約束的非線性優(yōu)化問題,對于m個輸入變量,有m個參數(shù)需要優(yōu)化,采用優(yōu)化算法尋優(yōu),獲得最優(yōu)的τ=[τ1,τ2,…,τm],τj∈[τmin,τmax],其中τmin和τmax通常由經(jīng)驗確定。KNN_DJMI互信息值的大小與Xi和Y的歷史數(shù)據(jù)長度w有關(guān),應選取合適的歷史數(shù)據(jù)長度w以覆蓋最相關(guān)的數(shù)據(jù)。

1.3 基于動態(tài)時間規(guī)整距離的在線更新策略

動態(tài)時間規(guī)整[20]是一種通過扭曲時間軸來更好地對時間序列的整體趨勢進行匹配的方法。假設有2個相同長度的時間序列W=w1,w2,…wn和L=l1,l2,…ln,為了計算2個序列的DTW距離,采用動態(tài)規(guī)劃的思路構(gòu)建代價矩陣Dn×n。該矩陣的組成元素為兩個序列元素的歐式距離,如式(11)所示:

Di,j=d(wi,lj)=wi-lj(7)

將求解序列的DTW轉(zhuǎn)化為動態(tài)規(guī)劃中求解代價矩陣中累計距離最小的路徑,通過式(8)所示的遞歸矩陣求解最短路徑:

ρ(1,j)=∑jk=1dw1,lk"""" j=1,2,…,nρ(i,1)=∑ik=1dwk,l1"""" i=1,2,…,nρ(i,j)=dwi,lj+min{ρ(i-1,j-1),ρ(i-1,j),ρ(i,j-1)}"""" i,j=2,…,n(8)

式中:ρ(i,j)為代價矩陣中(i,j)位置的累計最短距離, ρ(n,n)即為2序列的DTW距離。因此,采用DTW距離作為多步預測的評判指標更為合適。

在模型更新時,從式(3)可以得知,α值最小所對應的樣本在計算模型輸出時權(quán)重最小,其綜合作用最小,與其它樣本相比,最需要被新樣本替換,因此采用新樣本替換當前樣本集中的對應α值最小的樣本的在線更新策略。具體的,當模型需要更新的時候,將訓練集的數(shù)據(jù)與測試集中直至當前時刻的數(shù)據(jù)構(gòu)成一組新的時間序列數(shù)據(jù),對該數(shù)據(jù)進行分解,并得到當前時刻實際值的分解值,用當前時刻實際值的分解值替換對應的LSSVM的樣本集中α值最小的樣本,并更新模型。

1.4 算法流程

(1)首先對原始數(shù)據(jù)進行輸入變量的篩選以及數(shù)據(jù)預處理,并通過K近鄰聯(lián)合互信息法計算各輸入變量的遲延時間;

(2)然后通過Pearson相關(guān)系數(shù)確定分解模態(tài)數(shù)量,并用EMD對原始數(shù)據(jù)進行分解;

(3)針對不同頻域下的本征模態(tài)函數(shù)建立不同的LSSVM模型;

(4)將得到的各個模型的預測值重構(gòu),得到原始數(shù)據(jù)的預測值;

(5)計算預測值與實際值之間的誤差,如果預測誤差超過預設的閾值,對模型進行更新。對于基于EMD的自適應動態(tài)模型,模型的每個輸出是分解值而不是原始值,因此在更新的時候,將訓練集的數(shù)據(jù)與測試集中直至當前時刻的數(shù)據(jù)構(gòu)成一新的時間序列數(shù)據(jù),對該數(shù)據(jù)進行分解,并得到當前時刻實際值的分解值,把該分解值更新入各個模型中。

2 仿真實驗與結(jié)果分析

為驗證基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)預測建模方法在鍋爐燃燒系統(tǒng)上的效果,將所提的方法與LSSVM方法進行對比,分別從單步和多步預測2個方面進行比較。仿真對象為某650 MW的超臨界燃煤機組鍋爐燃燒系統(tǒng),并以NOx排放量模型為例,對其進行仿真分析。模型的輸入特征包括6層給煤量、6層二次風門開度、2層燃盡風門開度、1個省煤器后氧量共13個控制變量和機組運行負荷1個狀態(tài)變量以及過去時刻輸出量。選取Pearson相關(guān)系數(shù)閾值P為0.9,模態(tài)分解數(shù)選擇結(jié)果如表1所示,NOx排放量分解為4個本征模態(tài)和1個余項。通過K近鄰聯(lián)合互信息計算各控制變量的遲延時間,選擇NOx排放量分解信號的前兩個本征模態(tài)采用時間序列建模方法,預測模型結(jié)構(gòu)如圖2所示。

2.1 單步預測

圖1提供了鍋爐燃燒系統(tǒng)NOx排放量的單步預測結(jié)果。

由圖1可知,不進行分解的LSSVM模型的預測結(jié)果總是跟隨上一時刻的真實值,這是因為鍋爐燃燒系統(tǒng)的數(shù)據(jù)復雜度過高,模型未正確學習到輸入特征和輸出特征之間的對應關(guān)系,導致模型在訓練過程中易出現(xiàn)欠擬合現(xiàn)象,產(chǎn)生一個“滯后”現(xiàn)象,而引入EMD方法對輸出特征進行分解,NOx排放量預測值曲線和實際值基本保持一致,基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)模型未出現(xiàn)“滯后”現(xiàn)象。

2.2 多步預測

考慮在實際應用中提供長期預測結(jié)果的需求,還進行了多步預測的仿真實驗,預測步數(shù)設

定為10步,并采用DTW作為評判模型多步預測精

度的性能指標,將預測模型在線更新的觸發(fā)條件更改為預測得到的時間序列與實際值比較時DTW相似度超過預定的閾值,基于DTW指標的更新條件更側(cè)重于預測未來一段時間內(nèi)輸出趨勢變化的準確

性。未來10步的鍋爐燃燒系統(tǒng)預測結(jié)果如圖2所示。

從表3可以看出,基于EMD的自適應動態(tài)模型的DTW指標值相比于LSSVM模型有較大的降低,在多步預測中基于EMD的自適應動態(tài)模型仍能表現(xiàn)出較高的預測精度。

由圖2可知,LSSVM模型無法準確估計NOx排放量變化趨勢,有時甚至出現(xiàn)預測趨勢和實際趨勢相反的情況,而基于EMD的自適應動態(tài)模型則保持相對較好的預測結(jié)果。

3 結(jié)語

針對鍋爐燃燒系統(tǒng)這一復雜對象,提出了一種基于EMD的鍋爐燃燒系統(tǒng)自適應動態(tài)建模方法。使用EMD方法將原始的復雜時間序列信號分解成一系列子序列,降低數(shù)據(jù)的復雜度以匹配模型的復雜度;然后針對不同頻率的分解信號,分別不同的方法進行建模。為了兼顧預測精度和計算量,采用余項與原信號的Pearson相關(guān)系數(shù)來選取分解模態(tài)數(shù)量,針對鍋爐燃燒系統(tǒng)的時滯特性以及強耦合特性,采用K近鄰聯(lián)合互信息法綜合考慮控制量之間的相互影響,得出各控制量的遲延時間。為了準確判斷模型的多步預測能力,采用了DTW作為多步預測的性能評價指標。仿真結(jié)果表明,相比與已有的LSSVM方法,基于EMD的燃燒系統(tǒng)建模方法無論是單步預測還是多步預測,都具有更好的預測性能,為后續(xù)實現(xiàn)燃燒優(yōu)化控制打下了模型基礎,具有較高的工程實用價值。

【參考文獻】

[1] 鄧龍強,吳健成,何珦,等.生物質(zhì)散料噴燃鍋爐的燃燒特性研究[J].工業(yè)鍋爐,2022(6):23-26.

[2] 蘇濤,潘紅光,黃向東,等.基于改進PSO-SVM的燃煤電廠煙氣含氧量軟測量[J].西安科技大學學報,2020,40(2):342-348.

[3] 李楊,藍茂蔚,趙國欽,等.基于PCA-PSO-LSSVM的電站鍋爐效率預測模型研究[J].熱力發(fā)電,2021,50(12):43-50.

[4] 藍茂蔚,李楊,趙國欽,等.基于MAPSO優(yōu)化LSSVM的鍋爐燃燒建模研究[J].中南大學學報(自然科學版),2022,53(4):1506-1515.

[5] 朱鈺森,金曉明,張泉靈.基于聚類與加權(quán)連接的鍋爐NOx排放量多模型建模[J].控制工程,2019,26(4):688-693.

[6] 梁濤,靳云杰,姜文,等.基于改進MVO和WLSSVM的燃煤鍋爐NOx排放優(yōu)化[J].中國測試,2021,47(10):148-154.

[7] 李雪林,孫玉坤.基于改進果蠅算法優(yōu)化的最小二乘支持向量機動力鋰電池荷電狀態(tài)預測[J].南京工業(yè)大學學報(自然科學版),2023,45 (6):676-681.

[8] 張雲(yún)欽,程起澤,蔣文杰,等.基于EMD-PCA-LSTM的光伏功率預測模型[J].太陽能學報,2021,42(9):62-69.

[9] 胡小雷,劉小喜.經(jīng)驗模態(tài)分解法在PPK潮位測量中的應用[J].中國港灣建設,2020,40(9):17-21.

[10] 賈春雷,黃景濤.經(jīng)驗模態(tài)分解在水電機組振動信號分析中的應用[J].電工技術(shù),2020,No.

526(16):77-79.

[11] 李貴紅,杜昕,趙麗麗,等.基于EEMD-SVM的銑刀磨損狀態(tài)監(jiān)測系統(tǒng)設計[J].自動化與儀器儀表,2019,No.236(6):30-32.

[12] LV Y,YANG T,LIU J.An adaptive least squares support vector machine model with a novel update for nox emission prediction[J].Chemometrics amp; Intelligent Laboratory Systems,2015,145:103-113.

[13] HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition method and the Hilbert spectrum for non-stationary time series analysis[J].Proceedings of mathematical physical amp; engineering sciences,1998,454:903-995.

[14] 姚大春,張宏偉.火電機組節(jié)能改造技術(shù)研究[J].粘接,2021,47(9):155-159.

[15] 鄭亞妮.基于K氣體成分的積灰初始沉積層數(shù)學模型構(gòu)建研究[J].粘接,2020,41(5):47-51.

[16] 黃埔.四角切圓鍋爐的降維自適應燃燒優(yōu)化控制方法研究及應用[D].南京:東南大學,2021.

[17] 趙征,李悅寧,袁洪.燃煤機組NOx生成量動態(tài)軟測量模型[J].動力工程學報,2020,40(7):523-529.

[18] KRASKOV A,STOGBAUER H,GRASSBERGER P.Estimating mutual information [J].Physical Review E,2004,69(6):66-138.

[19] 康俊杰.電站鍋爐燃燒和SCR脫硝系統(tǒng)一體化建模與優(yōu)化控制研究[D].北京:華北電力大學(北京),2021.

[20] MYERSC,RABINERLR,ROSENBERGAE.Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[J].IEEE Transactions on Acoustics,Speech,and Signal Processing,1980,28(6):623-635.

主站蜘蛛池模板: 天天躁夜夜躁狠狠躁图片| 久久综合国产乱子免费| 欧美伦理一区| www精品久久| 亚洲免费播放| 亚洲AV无码一二区三区在线播放| 97国产精品视频自在拍| 国产香蕉一区二区在线网站| 国产精品美女在线| 免费在线成人网| 色婷婷电影网| 全午夜免费一级毛片| 亚洲欧洲综合| 日韩精品毛片| 国产无吗一区二区三区在线欢| 日韩大片免费观看视频播放| 欧美激情网址| 九色国产在线| 亚洲另类色| 国产成人精品一区二区不卡| 在线精品亚洲一区二区古装| 国产精品亚洲а∨天堂免下载| 美女裸体18禁网站| 国产成人欧美| 国产一级做美女做受视频| 国产精品高清国产三级囯产AV| 欧美亚洲另类在线观看| 欧美第一页在线| 亚洲精品无码AV电影在线播放| 国产成年无码AⅤ片在线 | 久久综合伊人 六十路| 极品尤物av美乳在线观看| 国产精品白浆无码流出在线看| 国产亚洲欧美另类一区二区| 久久国产精品波多野结衣| 69av在线| 久久精品日日躁夜夜躁欧美| 四虎国产永久在线观看| 色婷婷综合激情视频免费看| 欧美综合一区二区三区| 精品久久久久久成人AV| 天堂av综合网| 91视频首页| 欧美色香蕉| 国产综合亚洲欧洲区精品无码| 国产午夜福利亚洲第一| 亚洲综合婷婷激情| aⅴ免费在线观看| 无码精品国产VA在线观看DVD | 中文字幕无码电影| 日本欧美成人免费| 欧美一级高清免费a| 欧美色丁香| 在线va视频| 波多野结衣一区二区三区AV| 免费国产小视频在线观看| 在线日本国产成人免费的| 国产在线无码av完整版在线观看| 操操操综合网| 色网站在线免费观看| 免费人成在线观看视频色| 区国产精品搜索视频| 日韩不卡免费视频| 欧美午夜网| 视频二区中文无码| 中文字幕亚洲另类天堂| 亚洲中文字幕日产无码2021| 亚洲av无码久久无遮挡| 国产精品夜夜嗨视频免费视频| 国产中文在线亚洲精品官网| 亚洲成人在线网| 一本大道东京热无码av | 2021国产精品自拍| 久久国产精品波多野结衣| 影音先锋丝袜制服| 亚洲第一福利视频导航| 精品91自产拍在线| 亚洲国产成人久久精品软件| 亚洲第一福利视频导航| 欧美精品在线视频观看| 国产精品爽爽va在线无码观看| 一本无码在线观看|