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

基于Stacking策略的過程剩余執行時間預測

2020-01-14 09:36:22李帥標趙海燕陳慶奎
小型微型計算機系統 2019年12期
關鍵詞:活動方法模型

李帥標,趙海燕,陳慶奎,曹 健

1(上海理工大學 光電信息與計算機工程學院,上海市現代光學系統重點實驗室,光學儀器與系統教育部工程研究中心,上海 200093)2(上海交通大學 計算機科學與技術系,上海 200030)

1 引 言

為了提高運營效率,降低成本,提高客戶服務質量并降低人為錯誤的可能性,越來越多的公司采用過程感知信息系統(Process-Aware Information Systems,PAISs)來支持其業務流程,同時以事件日志的形式記錄每個執行活動的蹤跡.這一轉變致使企業積累大量數據,為了提高企業效益,從海量數據中提取有用信息預測業務過程實例這一趨勢得到很大發展.

目前已經提出很多用于預測業務過程實例行為的方法[1,2],比如預測業務過程風險或者業務過程剩余執行時間等.對于業務過程風險而言,這些預測可幫助過程管理者識別出具有風險的過程實例并提前給出解決方法以減少損失,比如資源重新分配或人工干預等.同樣,業務過程剩余執行時間對過程管理者或用戶體驗非常重要,比如客戶在銀行辦理業務時,提供預測結束時間對客戶滿意度有積極影響.

近年來,研究者提出了一些方法以對業務過程剩余執行時間進行預測.但是現有預測方法主要采用了黑盒模型,它們是通過歷史事件日志文件構建一個整體模型,然后將模型輸出值作為預測的剩余時間值.雖然現有技術能夠預測業務過程活動剩余時間值,但是沒有解釋后續活動的信息是如何影響預測的時間總值的.比如對一個簡單的串行業務過程實例A→B→C→D而言,假設當前業務過程執行到B,如需預測活動B的剩余時間,黑盒方法則利用歷史記錄訓練模型直接預測活動B到達活動D的剩余時間,而不考慮活動C對剩余時間的影響.一些學者在對業務過程時間預測研究的實驗結果表明考慮上下文信息對預測效果有積極影響[3],所以考慮活動C的相關特征應該能夠對預測結果準確性的提高有積極影響.另外,業務過程活動作為時間序列活動,活動之間存在的依賴關系對結果也有很大影響.

為了考慮各個活動對剩余時間的影響以及活動之間的關系對最終預測結果的影響,我們提出一種基于Stacking策略的模型融合(Stacking Strategy for Model Fusion,SSMF)的白盒方法用于業務過程剩余時間的預測.該模型由兩部分組成,第一部分引入樸素貝葉斯(Naive Bayes,NB)模型用于計算當前狀態選擇下一個活動的概率,引入支持向量回歸(Support Vector Regression,SVR)用于計算當前活動狀態到達下一活動狀態花費的時間;第二部分將當前活動預測所有可能的分支概率與分支時間作為訓練集,對長短期記憶(Long Short-Term Memory,LSTM)神經網絡模型進行訓練;最后模型預測結果即為當前活動剩余執行時間.

本文后續內容組織安排如下:第2節介紹近年來業務過程時間預測的研究方法;第3節介紹SSMF算法的基本思想、本論文中使用的相關機器學習技術以及SSMF模型的具體實現過程,第4節闡述相關實驗與分析,第5節是對全文進行總結以及未來工作方向的分析.

2 相關工作

近年來,業務過程的剩余時間預測問題已經被很多學者所研究.Aalst等人在文獻[4]中提出了第一個側重于時間視角的預測方法,該方法從歷史事件日志信息構建帶有標注的變遷系統(An Annotated Transition System),為了預測業務過程的剩余完成時間,使用狀態表示函數將到目前為止執行的活動序列映射到變遷系統中的一個狀態,并通過基于類似狀態下的較早的過程實例的平均完成時間來進行預測.文獻[5,6]中提出了一種在文獻[4]中描述技術的擴展方法,他們根據“上下文特征”對日志文件中的蹤跡進行聚類,然后,對于每個聚類,他們使用文獻[4]中描述的方法創建預測模型.對于新的運行實例,使用預測聚類樹對其進行聚類,然后使用屬于特定聚類的模型進行剩余時間的預測.在文獻[7]中,Bolt等人則提出一種稱為查詢目錄的新方法,為了計算新的部分蹤跡尾部(Partial Trace Tail)的剩余時間,將其與查詢目錄(Query Catalogs)中存儲的部分蹤跡尾部進行比較,將匹配的部分蹤跡尾部的平均剩余時間作為剩余時間預測值.文獻[8]中,Polato等人使用了一種ε-SVR機器學習算法對業務過程進行剩余執行時間的預測.該方法使用業務過程的歷史蹤跡活動屬性信息訓練ε-SVR模型,并通過編碼技術將部分蹤跡的特征轉換為適合模型輸入的格式,輸出的值即為估計的剩余執行時間.Navarin等人[9]提出了一種基于LSTM神經網絡模型對剩余時間預測的方法,其特點是學習長期的依賴關系.它將歷史完整蹤跡中各活動信息通過編碼技術轉換為支持LSTM輸入的特征向量進行模型的訓練,最后將未完成的蹤跡信息輸入模型之中進行時間的預測.

上述方法雖然考慮了活動特征對結果產生的影響,但是并沒有考慮單個活動對最終預測結果能夠產生的影響大小.另外,文獻[9]中雖然考慮了活動之間的依賴關系,但是其學習對象為已經發生的活動,與本文有所區別.本文中的LSTM神經網絡是從后續還未發生的業務過程活動對象中學習長期依賴關系.基于上述描述,我們提出一種SSMF模型來預測當前正在執行業務過程的剩余時間.另外,實驗結果表明該模型考慮更多的上下文特征,其預測結果具有更高的預測精度.

3 基于Stacking模型融合策略的業務過程剩余時間預測算法

3.1 算法思想介紹

在真實的業務流程執行過程中,當某個活動結束后,其后續可能包含多條分支路徑,為了展示單個活動或基本元素以及活動之間存在的“上下文”對整體結果的影響,我們的算法思想可分為兩步.第一步對單個活動對最終預測效果的影響程度進行學習.為了解決這個問題,我們根據NB與SVR預測當前活動選擇后續活動的概率與到達各活動的執行時間,對當前活動后續所有分支路徑上的活動進行同樣的操作以獲取當前活動點所有可能值(概率與時間值)組合.第二步依據已獲取的信息估計當前活動的剩余時間.通過對最大概率路徑時間預測、期望值時間預測以及真實執行路徑時間預測等幾種不同方法的分析,我們發現在每一個活動處,根據已經獲取的各種可能的信息組合對活動剩余時間的預測可能是一個非線性的模型.神經網絡能夠很好的解決非線性問題[10],同時能夠學習活動之間存在的依賴關系,所以我們選用神經網絡模型用于預測剩余時間.

本文所提模型目的用于預測當前活動剩余執行時間,由于只有部分活動得到了執行,后續的執行路徑還是不確定的,因此,模型首先使用NB與SVR預測其后續可能選擇活動的分支概率與執行時間,這是一種用于統計當前活動點后續所有可能特征的方式,為訓練后續模型做準備.其次,將當前活動點所統計的所有分支概率與各個活動之間的執行時間作為LSTM模型的輸入特征,并通過訓練好的LSTM模型用于預測當前活動點的剩余執行時間.

基于上述描述,我們提出一種基于Stacking模型融合策略的業務過程剩余執行時間預測方法.為了更好的理解算法基本思想,我們通過圖1對其作出解釋.

圖1 模型執行過程解析圖Fig.1 Analysis of the model execution process

如圖1所示,假設業務過程活動執行到活動A,且其后續可能的活動為B,C,D.第一步:根據已經訓練好的NB與SVR模型預測其后續可能選擇活動的分支概率與執行時間,其中P1,P2,P3值表示活動A選擇后續活動B,C,D的概率,T1,T2,T3值用于衡量從活動A到達后續活動B,C,D的估計時間.第二步:根據第一步計算的結果,將{P1,P2,P3,T1,T2,T3}作為輸入特征輸入活動A所對應的已訓練結束的神經網絡模型,其輸出結果即為活動A的預測剩余執行時間.

3.2 機器學習技術介紹

本節主要介紹在本篇論文中運用到的機器學習技術與相關知識.

3.2.1 樸素貝葉斯分類器

樸素貝葉斯分類器是一個基于貝葉斯定理應用的概率分類器,它可以預測給定樣本屬于特定類的概率.貝葉斯定理數學公式如公式(1):

(1)

其中P(A)表示A的先驗概率,P(B)表示B的先驗概率.樸素貝葉斯則是借助上述思想計算概率,具體實現過程如下:

1)假設T表示訓練樣本集合,每個樣本都有對應的類別,類別集合Y={y1,…,yk}.每個樣本由n維向量表示,X={x1,…,xn}用于描述n個屬性的n個測量值.

2)給定樣本X,樸素貝葉斯分類器借助公式(2)預測X屬于具有最高后驗概率的類.

(2)

其中,P(X)對于訓練集中所有類別而言,其值是一樣的,因此計算過程只需考慮P(X|Yi)P(Yi)的結果.另外,樸素貝葉斯分類器假設屬性值對給定類的影響與其他屬性值無關.故P(X|Yi)可轉化為公式(3):

(3)

3)對于樣本X,根據公式(4)找出所有類別中最大后驗概率值,并將該類別作為樣本X的預測類別.

P(yi|X)=max{P(y1|X),…,P(yk|X)}

(4)

3.2.2 支持向量回歸

為了說明支持向量回歸機的基本思想,通過線性回歸函數對其簡要介紹.表達形式見公式(5).

(5)

(6)

一般引入拉格朗日函數將優化問題轉為其對偶形式:

(7)

(8)

對于不能用線性函數進行回歸的訓練數據,利用核函數將訓練集映射到高維空間之后再進行線性回歸.

3.2.3 LSTM神經網絡

LSTM神經網絡是一種改進的循環神經網絡(Recurrent Neural Networks,RNN)人工智能算法,LSTM對時間序列相關數據的處理有很好的效果[12].LSTM模型可以記憶更長的上下文[13].目前圖2中的LSTM細胞應用最為廣泛[14].

圖中123標號分別表示遺忘門、輸入門以及輸出門.

圖2 LSTM隱藏層細胞結構Fig.2 LSTM cell structure in hidden layer

遺忘門:其功能用來決定從細胞狀態中拋棄何種信息.該門讀入當前時間點的輸入xt以及上一層細胞的輸出ht-1,通過公式(9)計算ft值并將其賦值給當前細胞狀態Ct-1.其中1表示保留其信息,0表示舍棄其信息.

ft=σ(Wf*[ht-1,xt]+bf)

(9)

輸入門:其功能用來決定多少新的信息添加到細胞狀態中.該過程包含兩部分,第一部分:通過sigmoid函數決定哪些輸入需要添加到細胞狀態中,另外通過tanh函數生成候選向量值用來更新細胞狀態.第二部分:更新舊細胞狀態Ct-1為新狀態Ct.計算公式見公式(10)-公式(12).

it=σ(Wi*[ht-1,xt]+bi)

(10)

Ctemp=tanh(WC*[ht-1,xt]+bC)

(11)

Ct=ft*Ct-1+it*Ctemp

(12)

輸出門:其功能用來決定最終輸出哪些信息.該門運行過程分為兩步,第一步:通過公式(13)決定細胞狀態Ct中的哪些信息被輸出;第二步:對細胞狀態使用tanh處理得到結果并與第一步的結果相乘,以確定最終輸出的信息,計算公式見公式(14).

ot=σ(Wo*[ht-1,xt]+bo)

(13)

ht=ot*tanh(Ct)

(14)

3.3 SSMF模型框架

模型融合通過對不同模型結果的融合得到一個新的結果,該結果將取代各個原始模型的預測結果.Stacking策略是由Wolpert[15]等人提出的一種對分類器疊加組合思想的方法.該方法首先提取有效特征進行第一次的訓練和預測,并將該操作視為第一層,其次將第一層訓練的結果采用合適的融合方法進行二次訓練與預測,即第二層的訓練與預測,以獲得更為精確的預測結果[16].本文利用Stacking策略的思想,出于對單個活動對結果的影響以及活動之間存在的依賴關系的考慮,提出一種業務過程剩余時間預測算法,其模型結構如圖3所示.

圖3 Stacking模型融合結構圖Fig.3 Stacking model fusion structure

圖3表明本文所提模型共包含兩層結構:第0層結構與第1層結構.第0層結構由樸素貝葉斯分類器與支持向量回歸構成,第0層結構里包含的模型通常被稱為基分類器或基回歸器,基分類器與基回歸器的輸出結果將作為第1層的輸入數據.第1層結構由LSTM神經網絡構成,通過對第0層的輸出結果進行訓練以預測剩余時間.因該層輸入的數據是一種通過融合方法產生的元數據,故第1層也被稱為元數據分類器或回歸器.在獲取第0層的輸出結果后,需要采用一定的方法對輸出結果進行融合,融合方法可以分為兩種,固定融合方法(Fixed Rules,FR)和可訓練融合方法(Trained Rules,TR)[17].本文模型中采用的是一種流行的元學習(Meta-learning)可訓練融合方法,它是一種將上一層輸出結果作為中間特征[18](元數據)的方法.即對第0層的各個輸出結果進行組合作為第1層模型的輸入特征.

上述框架共經過兩次訓練,兩次訓練即相互獨立又相互關聯.獨立指的是兩次訓練都是單獨訓練互不影響,關聯指的是上一層模型的輸出作為下一層模型的輸入.

3.4 SSMF模型實現過程

本小節對模型的實現過程進行詳細介紹,首先介紹SSMF中各個子模型的運行過程,然后對其整體運行流程進行描述.

3.4.1 NB模型的使用

樸素貝葉斯分類器應用于SSMF模型第0層,首先將已經預處理的數據劃分為訓練集與測試集;然后以業務過程活動開始執行到當前活動獲取的信息作為訓練輸入特征,將當前活動真正選擇下一活動作為訓練輸出目標,訓練當前活動NB模型,并用于預測當前活動選擇下一可能活動的概率.業務過程活動繼續執行,為了提高預測的精度,需要將已執行的活動特征添加到數據集中,采用上述方式統計各個活動點可能選擇下一活動的概率,直至業務過程活動執行結束.

3.4.2 SVR模型的使用

支持向量回歸器應用于SSMF模型第0層,首先將已經預處理的數據劃分為訓練集與測試集;然后以業務過程活動開始執行到當前活動獲取的信息作為訓練輸入特征,將當前活動執行到下一活動的真實時間作為訓練輸出目標,訓練當前活動SVR模型,并用于預測當前活動執行到下一可能活動的時間.業務過程活動繼續執行,為了提高預測的精度,需要將已執行的活動特征添加到數據集中,采用上述方式統計各個活動點可能到達下一分支活動的執行時間,直至業務過程活動執行結束.

3.4.3 LSTM模型的運用

長短期記憶神經網絡應用于SSMF模型第1層,該層主要應用于預測當前活動的剩余執行時間.根據3.4.1小節與3.4.2小節中所述,在當前業務過程活動狀態處可統計其后續所有可能選擇的分支概率與執行時間,并將其組合結果作為訓練LSTM模型的輸入特征,將業務過程當前活動真正結束時間作為訓練LSTM模型的輸出目標.通過上述方式,訓練不同活動點的時間預測模型.

3.4.4 整體過程描述

1)獲取有效特征.根據已有的事件日志文件,通過數據預處理、特征篩選、特征提取等步驟獲取對活動時間有影響的活動特征.

2)訓練模型第0層中樸素貝葉斯分類器.根據3.4.1小節運行過程所述對不同活動點處樸素貝葉斯分類器進行訓練,并將訓練好的模型應用于測試集中,預測當前活動選擇下一個可能執行活動的概率.

3)訓練模型第0層中支持向量回歸器.根據3.4.2小節運行過程所述對不同活動點處支持向量回歸器進行訓練,并將訓練好的模型應用于測試集中,預測當前活動到達下一個可能執行活動的時間.

4)訓練模型第1層中LSTM神經網絡模型.當步驟(2-3)執行結束后,采用3.4.3小節所介紹訓練方式,對業務過程不同活動點處訓練各自的神經網絡模型.

5)預測剩余時間.對測試集的業務過程實例,利用對應活動點已經訓練好的NB模型與SVR模型用于預測其可能選擇下一活動的概率與到達下一活動的執行時間,并將統計的信息作為對應活動點神經網絡模型的輸入特征,其輸出結果即為預測當前活動剩余執行時間值.

4 實驗與結果分析

4.1 實驗數據

本文所使用的數據集來自于2017年業務過程智能(Business Process Intelligence,BPI)挑戰賽,該數據集包含了2016年1月-2017年2月期間荷蘭金融機構貸款申請流程的真實業務流程,其特征包含活動類型,活動執行者,時間戳以及其他信息,比如申請的金額和貸款目的等.該數據集中包含多次創建報價的業務過程實例,為了簡化實驗過程,故只考慮創建一次報價且信息驗證不超過兩次的業務過程實例.整理后的實驗數據包含22900條貸款申請實例與169457個活動.

4.2 實驗數據分析

貸款申請過程由一系列活動構成,數據集中貸款申請可能的執行路徑如圖4所示.

圖4表示當前數據集中所有可能執行的路徑及結果,其中相關活動名稱介紹見表1,后續相關實驗以該過程模型為基礎.

在實驗部分,我們假設當前業務過程活動執行到活動A_Concept,在此基礎上預測過程實例在活動A_Concept的剩余執行時間以及后續活動節點的剩余執行時間.

圖5描述的是數據集中業務過程活動中兩個相鄰活動時間間隔的分布,其中Duration_1表示活動A_Concept到活動A_Accepted的時間間隔分布,后續活動間隔以此類推.由圖觀察可知,不同活動間隔的時間分布有所差異,這種不同活動時間間隔差異可看做活動之間存在某種關系.神經網絡在不同活動之間學習長短期依賴信息有很好效果,這也是本文使用神經網絡作為模型訓練的原因.

圖5 活動間隔時間分布Fig.5 Activity interval distribution

4.3 評價指標

為了檢驗模型的有效性,本文使用均方根誤差(Root Mean Squared Error,RMSE)作為評測指標.RMSE計算實際觀測值與預測值之間的標準差,其值越小表示模型預測的結果要好.計算公式如下:

(15)

4.4 方法比較

為了驗證所提方法的預測精度,將其與部分業務過程剩余時間預測方法進行比較.在本文中使用歷史平均值(Mean Value,MV)方法、最大概率路徑(Maximum Probability Path,MPP)預測方法、真實執行路徑(Real Execution Path,REP)預測方法、期望值(Expected Value,EV)預測方法以及支持向量回歸方法與所提出方法進行比較.

表1 活動名稱說明
Table 1 Activity name description

活動名稱說明A_Submitted客戶提交新的貸款申請A_Concept系統自動完成第一次評估A_Accepted對申請的客戶提供貸款金額A_Complete貸款金額發送給客戶,銀行等待客戶退回已簽名文件(工資單或身份證等)A_Validating銀行對收到的文件進行檢查A_Incomplete文件信息缺失,需要客戶重新發送文件A_Pending銀行收到所有文件并且評估結果為正,為客戶提供貸款A_Cancelled客戶從未發送過文件或電話咨詢客戶被告知不需要貸款,貸款申請被取消A_Denied客戶申請條件不符合銀行檢驗標準,貸款申請被拒絕

對上述部分比較方法進行說明:

1)MPP方法對業務過程活動中每個活動使用NB與SVR模型預測當前活動選擇下一個可能活動的概率與達到下一個活動的執行時間.在每個活動點上利用概率乘積的方式統計每個活動點所有可能執行路徑的概率,并選取概率值最大的活動組合作為當前活動執行路徑,并將該路徑中各個活動預測執行時間求和作為其預測值.

2)REP方法利用SVR模型預測活動執行時間,將當前業務過程實例真正的執行路徑作為預測執行路徑,并計算各個活動執行時間總和作為預測總值.

3)EV方法首先統計當前業務過程實例在各個活動點選擇不同結果的概率與剩余執行時間,到達不同結果的概率與剩余執行時間計算采用加權求和方式,然后根據前者計算的結果,再次利用加權求和方法求取不同活動點的剩余執行時間.雖然兩次計算方法均采用加權求和,但是所執行的對象不同.

4)SVR方法將業務過程執行開始到當前活動已獲取的特征作為輸入,將當前活動執行至業務過程結束的真實時間作為輸出進行模型訓練.對訓練結束的模型,輸入已獲得活動特征,輸出值即為當前活動的預測剩余時間.

4.5 對比試驗分析

本文搭建的SSFM模型運用到樸素貝葉斯、支持向量回歸以及LSTM神經網絡.其中LSTM參數部分對結果影響較大,因此對LSTM神經網絡部分進行介紹.本文搭建的LSTM神經網絡模型包括輸入層、隱藏層以及輸出層,其中誤差函數使用均方誤差(mean squared error,mse),優化算法使用適應性動量估計(Adaptive moment estimation,Adam)算法,它是由Kingma等人提出[19].Adam是一種基于梯度的隨機優化算法,其優點是能夠自適應性更新學習率,且在實際應用中效果很好[20].

為了防止模型學習不充分或過擬合,我們在不同迭代次數處輸出損失函數值以選擇合適的迭代次數,各個活動迭代次數損失函數輸出值如圖6所示.

圖6 不同迭代次數的損失函數值Fig.6 Loss function values for different iterations

圖6表示不同活動點迭代次數與損失函數值的關系,通過圖6(a)發現在活動A_Concept處當迭代次數達到17500左右,損失函數值趨于穩定,故將17500設置為當前活動的迭代次數.同理,由圖6(b)-圖6(f)觀察可知,后續活動迭代次數分別設置為15000,15000,8000,6000,3000.

本文通過上述參數學習不同活動點網絡模型,并對數據集按照8:2劃分為訓練集與測試集對算法進行試驗,本文所提方法與對比實驗結果如表2所示.

表2 不同預測方法的RMSE值
Table 2 RMSE values for different prediction methods

活動名稱方法MPPREPEVMVSVRSSMFA_Concept13.2998.79114.77111.12910.8068.701A_Accepted12.8538.30710.69710.64710.2888.269A_Complete12.2898.0929.37210.54410.0028.160A_Validating5.6904.9375.4975.7586.4365.140A_Incomplete5.0634.9075.0595.0866.2274.756A_Validating2.4562.1872.4572.4172.5972.360

對比實驗REP將業務過程的真正執行路徑作為其預測執行路徑,而其他對比實驗則不知道真正的執行路徑,故其實驗結果可作為其他方法的參考標準, 即其他方法的實驗結果與REP結果越接近表明該方法預測效果越好.由表2可知,本文提出的SSFM算法在不同活動點的預測精度都比其他對比算法較高.對所有方法而言,其預測誤差隨著活動執行逐漸變小,進一步說明考慮更多活動上下文對預測結果有積極影響.另外實驗結果表明SSMF預測結果與REP預測結果最為接近,且在個別活動點其預測效果更好,即表明SSFM預測結果與業務過程實例的真正執行路徑時間較為接近,更好的解釋了SSFM模型的學習活動之間依賴信息的能力以及考慮各個活動對整體結果具有積極作用.

5 總結與展望

本文結合樸素貝葉斯、支持向量回歸以及LSTM神經網絡模型,提出基于Stacking策略模型融合的業務過程剩余執行時間預測方法.該方法具有以下優點:1)該方法考慮更多的上下文特征對預測結果的影響;2)該方法考慮基本組件或單個活動對整體結果所做的貢獻與影響;3)該方法考慮各個活動之間存在依賴關系.另外,將本文所用方法與相關業務活動剩余時間預測方法進行比較,在以RMSE作為評價指標的前提下,發現該模型具有更高的預測精度,具有更好的表現效果.本文使用的模型只應用于預測當前活動的剩余執行時間,業務過程執行路徑由一系列活動組成,其執行路徑可看做一系列序列活動,所以未來的計劃將研究擴展該模型,使其能夠完成業務過程執行路徑的預測.

猜你喜歡
活動方法模型
一半模型
“六小”活動
少先隊活動(2022年5期)2022-06-06 03:45:04
“活動隨手拍”
行動不便者,也要多活動
中老年保健(2021年2期)2021-08-22 07:31:10
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
三八節,省婦聯推出十大系列活動
海峽姐妹(2018年3期)2018-05-09 08:20:40
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 日韩午夜福利在线观看| 又粗又大又爽又紧免费视频| 国产va免费精品观看| 91国语视频| 91青青视频| 日韩无码视频专区| 亚洲一级色| 一级毛片免费的| 亚洲大尺度在线| 岛国精品一区免费视频在线观看 | 国产精品v欧美| 欧美国产日韩在线播放| 欧美福利在线观看| 亚洲第一天堂无码专区| 午夜毛片免费看| 亚洲性一区| 欧美丝袜高跟鞋一区二区| 中文字幕亚洲乱码熟女1区2区| 国产女人18水真多毛片18精品| 欧美成人一级| 午夜啪啪福利| 欧美一区精品| 国产欧美日韩另类| 日本在线免费网站| 日韩免费无码人妻系列| 69视频国产| 在线免费a视频| 亚洲国产精品日韩av专区| 呦女精品网站| 久久综合九九亚洲一区| 97视频免费在线观看| 国产成在线观看免费视频| 久久无码免费束人妻| 亚洲AⅤ综合在线欧美一区| 91无码人妻精品一区| 国产美女精品人人做人人爽| 国产主播在线一区| 国产成人精品优优av| 婷婷综合缴情亚洲五月伊| 国产精品白浆无码流出在线看| 欧美不卡在线视频| 亚洲青涩在线| 超清无码一区二区三区| 天堂av综合网| 欧洲亚洲一区| 九九这里只有精品视频| 亚洲日本在线免费观看| 国产乱子伦视频在线播放| 激情无码字幕综合| 欧美日本视频在线观看| 欧美日本二区| 精品视频一区在线观看| 国产情侣一区二区三区| 成人字幕网视频在线观看| 久久婷婷国产综合尤物精品| 色综合久久88| 久久成人18免费| www成人国产在线观看网站| 视频国产精品丝袜第一页| 无码人妻免费| 久久男人资源站| 久久99蜜桃精品久久久久小说| 国产一在线| 国产白浆一区二区三区视频在线| 久久人体视频| 欧美国产三级| 在线视频亚洲色图| 啊嗯不日本网站| 精品视频第一页| 久青草国产高清在线视频| 国产主播在线一区| 无码国产偷倩在线播放老年人| 欧美激情综合| 午夜国产在线观看| 无码中文AⅤ在线观看| 午夜国产大片免费观看| 自拍中文字幕| 亚洲一区波多野结衣二区三区| 97成人在线观看| 亚洲va视频| 亚欧乱色视频网站大全| 国产欧美亚洲精品第3页在线|