?

數據驅動的含IIDG配電網短路電流計算多輸出模型

2022-09-14 08:52葉睿愷王慧芳張亦翔
電力自動化設備 2022年9期
關鍵詞:測量點短路準確性

葉睿愷,王慧芳,張 森,張亦翔

(浙江大學 電氣工程學院,浙江 杭州 310027)

0 引言

我國“雙碳”戰略目標促使逆變型分布式電源IIDG(Inverter Interfaced Distributed Generator)在配電網中逐漸呈現高占比的發展趨勢[1]。與傳統同步電機電源不同,IIDG 具有強非線性的故障特性,不僅要滿足低壓穿越要求[2],還受故障前運行特性、接入位置、控制策略等因素影響。大量IIDG 的接入使得傳統短路電流計算方法難以直接應用,迫切需要對含IIDG配電網短路電流計算方法進行研究[3]。

針對IIDG 故障特性以及短路電流計算方法,國內外學者已開展了很多研究。文獻[4]詳細介紹了IIDG 的低壓穿越控制策略,并給出了短路電流的表達式;文獻[5-6]根據不同目標的控制策略分析了不對稱故障下IIDG 的電流輸出特性;文獻[7]分析了IIDG 并網后的故障特性,得到了對稱故障與不對稱故障下短路電流計算模型;文獻[8]分析了IIDG 并網后對配電網的影響,給出了含少量IIDG 配電網短路電流計算通用模型;文獻[9]將IIDG 分為故障點上游和故障點下游并分別建立故障等效模型,采用基于疊加原理的迭代短路電流計算方法。上述研究均是基于機理的物理建模,不同控制策略下IIDG 模型不同,不同故障類型下不同的短路電流計算方法都需要迭代計算。當IIDG 大量接入或網絡規模增大時,上述方法因迭代計算耗時必將大增,甚至可能出現不收斂的情況。為提升計算速度,有學者對IIDG 模型做極簡處理,例如用1.2~2 倍額定電流來表征,但這必將犧牲計算準確性[10]。因此現有物理建模方法應用于IIDG 占比高的大型配電網時,會出現計算快速性與準確性之間的矛盾。

為了解決該矛盾,文獻[11]基于機器學習思想提出了機理與數據融合驅動的含IIDG 配電網短路電流計算方法,并與詳細的物理建模方法進行對比,結果表明文獻[11]所提方法的計算誤差接近物理建模方法的計算誤差或更小,且計算速度快。同時,電網規模、IIDG 接入數量對該方法的計算性能影響較?。?2]。然而,該方法使用單輸出的學習模型,當需要計算配電網中多個測量點電流值時,必須針對每個測量點訓練各自的學習模型,因此存在模型數量多的問題。

為解決模型數量問題,本文提出了含IIDG 配電網短路電流計算多輸出模型。對不同的多輸出模型進行分析與對比,并根據獲得的配電網運行狀態的實際條件,選擇2 種包含運行方式和故障信息的輸入特征,分別進行模型訓練,最終實現了同時輸出配電網所有支路短路電流的目標。

1 多輸出回歸模型選擇

傳統機器學習的目標多為建立單輸出模型,即采用訓練完成的模型對特定目標進行回歸計算。實際問題中往往需要同時解決多個任務,且多個任務之間存在不同的相關性,因此需要研究多輸出回歸計算。根據多輸出回歸的子任務關聯方式,大致可分為問題轉化和算法適應方法2類[13]。

1.1 問題轉化方法

多輸出任務目標是建立M個輸入與D個輸出之間的映射關系。問題轉化的核心思想是將相應的映射關系轉化為D個單輸出的形式,然后將D個單輸出模型拼合成1 個多輸出模型,達到同時輸出多個預測值的目的。問題轉化的算法包括單目標方法ST(Single-Target method)、多目標回歸模型融合MTRS(Multi-Target Regressor Stacking)方法和回歸鏈RC(Regressor Chain)方法等。

設有N組樣本,Xl為第l組的輸入,Yl為第l組對應的輸出,即第l組樣本的輸入、輸出為:

第二階段通過將第一階段的輸出預測值拼合作為訓練集的輸入,在一定程度上考慮了多個輸出之間的潛在關系,對問題轉化預測值偏差進行了一定的修正,可以提升預測準確性。但是,MTRS 方法本質上與ST 類似,都是直接考慮輸入和輸出之間的映射,對于輸出之間的潛在關系的考慮并不全面。此外,MTRS 方法需要進行2 次ST 計算,使得模型偏大、計算速度慢。

RC 方法是基于馬爾可夫隨機鏈的思想,根據確定的鏈順序對每個輸出建立獨立的回歸模型。當有D個輸出值時,假設指定鏈順序C為:

最終將D個預測值一起輸出。由上述過程可見,RC 方法將不同輸出考慮為不同狀態,將不同狀態以輸入的方式挖掘多個輸出之間的相關性,因此考慮了不同輸出之間的關系。針對不同的實際問題,應根據輸出任務間的相關性選擇不同的問題轉化方法,達到最佳的多輸出性能。含IIDG 配電網短路電流計算問題中,由于IIDG 與系統電源提供的短路電流相位往往不同[8],各支路穩態電流的幅值之間不再顯式地滿足基爾霍夫電流定律,但電流相量之間還是滿足基爾霍夫電流定律的,即輸出之間存在一定的相關性,因此將其他模型的輸出量納入輸入的MTRS和RC方法這2種多輸出模型方法適用于配電網短路電流多輸出計算,可充分挖掘不同支路電流間的潛在關系。

由于MTRS 和RC 方法是將單輸出模型拼合為多輸出模型,因此方法的性能與作為基學習器的單輸出模型性能密切相關。配電網短路電流計算由于輸入特征數量較多,使模型復雜程度增大;同時輸出電流數量多且有部分數值相近,因而存在一定的過擬合風險。為此本文選擇近年來在各大比賽中性能表現優秀的基學習器——極限梯度提升機XGBoost(eXtreme Gradient Boosting)[14]和輕量梯度提升機LightGBM(Light Gradient Boosting Machine)[15]。XGBoost 和LightGBM 都是梯度提升決策樹GBDT(Gradient Boosted Decision Tree)的改進與高效實現,其中XGBoost 采用的是按層生長策略,能夠同時分裂同一層的葉子,從而進行多線程優化,可有效降低過擬合風險;LightGBM 采用帶有深度限制的按葉子生長策略,每次從當前所有葉子中找到分裂增益最大的葉子進行分裂,在分裂次數相同的情況下,該策略可以降低更多的誤差。此外,LightGBM 使用了基于梯度的單邊采樣對輸入進行降維以降低模型復雜程度,因此其可以在保證計算準確度的前提下降低模型復雜程度,使得模型計算速度更快。

但是,MTRS 方法每進行一組樣本訓練都需要更新一次輸入特征,RC 方法則需要更新D-1 次輸入特征,這可能導致LightGBM 速度慢于XGBoost。原因是LightGBM需要對新輸入特征使用互斥特征捆綁EFB(Exclusive Feature Bundling),當輸入特征不變或變化次數少時,EFB 能大幅提升算法性能,當輸入特征變化次數多時,EFB 有可能減慢速度。因此,在MTRS方法下LightGBM的速度可能比XGBoost快,在RC 方法下LightGBM 的速度可能比XGBoost 慢。

綜上所述,含IIDG 配電網短路電流計算多輸出模型采用問題轉化方法時,分別采用XGBoost 和LightGBM 作為基學習器,并對比基于MTRS和RC方法的多輸出模型性能。

1.2 算法適應方法

算法適應方法是將單輸出模型根據實際問題進行修改,用1 個模型同時輸出多個預測值。算法適應方法的輸出映射不僅需要考慮輸入變量,更需要考慮不同輸出量之間的關系,以此提升多個預測值的準確性。

算法適應方法包含統計學方法、多輸出支持向量機、多輸出回歸樹、核函數方法、神經網絡等多種方法。算法在修改之前是適用于單輸出模型的,針對不同的實際問題,需要修改算法中的輸入、模型結構、輸出標簽、損失函數等參數,不同的算法適應方法修改方式并不完全相同。本文采用應用較為廣泛的神經網絡作為多輸出回歸預測的基礎算法。

神經網絡一般是指全連接的人工神經網絡ANN(Artificial Neural Network),它將輸入層的神經元全部連接至隱藏層,利用所有輸入特征計算得到輸出值,因此是利用輸入特征最充分、計算結果較為準確的模型。然而,全連接人工神經網絡對訓練樣本數量要求高、模型計算量大。卷積神經網絡CNN(Convolutional Neural Network)可以解決以上問題,其卷積層可以通過更少的計算步驟實現類似效果,對于復雜模型更有優勢。CNN 的卷積層通過滑窗實現輸入特征的局部感知,即相鄰輸入數據可以為輸出提供更高價值。電網具有基爾霍夫電流、電壓定律,因此將輸入的電氣量特征相鄰排列,可使CNN 的卷積核獲得有益信息,進而提升短路電流計算準確性。

本文選用ANN 和CNN 作為算法適應方法的代表算法,驗證算法適應多輸出模型應用于含IIDG 配電網短路電流計算的有效性。

2 多輸出模型的短路電流計算方法

2.1 輸入特征選擇

目前,IIDG 高占比的配電網示范工程已配置有微型同步相量測量單元μPMU(Micro-synchronous Phasor Measurement Unit),可直接獲得電壓幅值、電壓相位、有功、無功等數據,這些數據量反映了配電網的運行方式。然而,大部分配電網不具備上述條件,但一般在各支路上安裝有電流互感器,能獲得支路電流的有效值,并在系統電源和IIDG 接入處安裝有電壓互感器獲得電壓的有效值,這些測量值也能反映配電網的運行方式??紤]配電網實際測量條件,分別采用上述2 類電氣量作為配電網運行方式特征信息,再結合故障位置及過渡電阻,構成樣本的輸入特征。故障類型則不作為輸入特征,因為不同故障類型下配電網電氣量特征差異較大,樣本混合訓練會降低模型性能,為此,與傳統短路電流計算習慣一致,針對不同故障類型調用不同的計算模型。樣本的標簽為對應上述輸入特征的各條線路短路電流穩態值。

綜上所述,對于包含K個節點、E條支路的配電網,假設有1 個系統電源接入節點和H個IIDG 接入節點,則有2類短路電流計算模型輸入特征,即X1對應配置μPMU 的配電網,X2對應未配置μPMU 的配電網,分別如式(8)和式(9)所示。

式中:下標k為節點編號,節點1 為系統電源接入節點;下標e為支路編號;|V|、θ、P、Q、I、V分別包含三相電壓幅值、電壓相位、有功、無功、支路電流、電壓有效值數據,因此能應用于不對稱配電網;Fline為故障線路編號;Floc為故障位置距離線路首端的位置百分比,Floc∈[0,100%),當Floc=0 時,認為故障發生在首節點上;Ron為過渡電阻,發生接地故障時為相與地之間的電阻,僅相間故障時為相與相之間的電阻;SDGh為第h個IIDG的容量。

采用基于CNN 的算法適應方法時,為了更好地發揮CNN 的局部敏感特性,將輸入數據特征由一維修改至二維。具體是將反映運行方式的特征數據分為不同數據層,式(8)中分為了電壓幅值層、電壓相位層、有功層和無功層,這些層以節點或支路編號對齊,且每層有3 行分別對應三相。故障信息與IIDG信息單獨為1層,該層為稀疏矩陣,僅有故障和IIDG對應的節點或支路編號有相應信息,其余均為0。式(9)類似,只是將反映運行方式的數據分為了電流層、電壓層。

2.2 評價指標

計算準確性一般由真實值yj與計算值y?j的差異來衡量,其中真實值通常由仿真獲得。通常采用平均絕對誤差MAE(Mean Absolute Error)、平均絕對百分比誤差MAPE(Mean Absolute Percentage Error)來衡量模型準確度。MAE 和MAPE 指標都受真實值影響:當真實值大時,MAE 受影響大;當真實值小時,MAPE受影響大。因此用MAE和MAPE這2個指標綜合反映算法性能,避免真實值過大或過小帶來的影響。

多輸出模型中第j個輸出的MAE值δj定義為:

式中:λAPE為APE值;λj為第j個輸出的MAPE。

上述評價指標的數值越小,模型準確性越高。多輸出模型會產生D個δj和λj,為方便展示,文中選 用D個δj、λj指標中的最大值δMAX、λMAX和平均值δMEAN、λMEAN。

2.3 模型參數尋優

對于機器學習而言,超參數是模型訓練過程中需要重點調整的內容,獲得最優的超參數可以提升模型性能,使評價指標δMEAN和λMEAN達到最小。不同類型的多輸出模型采用的超參數尋優方式不完全相同。

對于問題轉化方法,本文采用隨機采樣網格搜索與5 折交叉驗證結合的方式進行超參數尋優,以最終輸出的λMEAN作為評價標準,驗證超參數尋優結果。XGBoost 最佳超參數組合為:基分類器為400個,樹最大深度為10,子節點權重為2,子采樣系數為0.8,特征采樣系數為0.8。LightGBM 最佳超參數組合為:基分類器為430 個,葉子數為20,子采樣系數為0.8,特征采樣系數為0.9。

對于基于神經網絡的多輸出模型,超參數更為復雜,為避免網格搜索陷入局部最優,調參選用效率更高的貝葉斯調參法。本文使用3層全連接的ANN,每一層都需要使用激活函數,激活函數選擇ReLU函數,學習率為10-5,單次訓練樣本100 個。CNN 最優參數組合為:卷積核大小為3×3,填充半徑為1,激活函數為ReLU,池化層大小為2,學習率為10-5,單次訓練樣本150 個。CNN 輸出層需要將前層卷積輸出展開,使用1 層全連接層連接,采用ReLU 激活獲得最終結果?;谏窠浘W絡的多輸出模型均使用λMEAN作為損失函數。

2.4 多輸出模型計算流程

含IIDG 配電網的短路電流多輸出計算模型包括訓練和應用兩部分。其中,訓練需要大量樣本,通過仿真獲得樣本的方法同文獻[11],有所不同的是,提取的輸入特征有差異,且該文獻的樣本標簽為1個測量點,而本文的樣本標簽為所有支路測量點。具體流程如下:

1)通過MATLAB/Simulink 建立所需的配電網模型;

2)在常見運行方式區間內,隨機設置配電網運行方式和故障信息;

3)進行故障仿真,并記錄輸入特征和各條支路的短路電流穩態值作為標簽值。

重復上述步驟2)、3)即可獲取1 組訓練樣本。訓練樣本的數量可以根據訓練效果進行補充或調整。為了減少不同輸入特征數值大小對模型訓練的影響,對輸入數據進行歸一化預處理。模型訓練過程本質是調參過程,尋找模型最優參數使得輸出結果誤差最優,達到預期則訓練結束。模型應用時,輸入配電網的運行方式特征信息、故障信息,根據故障類型調用訓練好的對應模型,則可輸出所有支路的短路電流值。

3 算例驗證

3.1 算例情況

含IIDG 配電網結構圖如圖1 所示,按照圖1 以IEEE 34 節點配電網系統為基礎接入IIDG,IIDG 容量設定為:SDG1=500 kW,SDG2=200 kW,SDG3=500 kW,SDG4=400 kW。節點800 為系統電源接入點,考慮最大和最小2種基礎運行方式,等值阻抗分別為j0.5 Ω和j1 Ω。配電網電壓等級為24.9 kV,系統線路共有31 條,負荷參照文獻[16]設定,在最大、最小運行方式下不發生改變。

圖1 含IIDG配電網結構圖Fig.1 Structure of distribution network with IIDGs

隨機設置配電網運行方式,模型參數設置同文獻[11]:等值電源阻抗在系統最大、最小等值阻抗區間內隨機產生;IIDG 容量及負荷大小則在其原值的[80%,120%]范圍內隨機產生;設置一定概率的隨機切除部分線路或IIDG。隨機設置故障線路編號Fline、故障發生位置Floc和過渡電阻Ron。上述運行方式和故障信息的設置,均會通過式(8)或式(9)的輸入特征反映出來。需說明的是,過渡電阻的設置范圍應與故障類型相吻合,例如配電網單相接地短路故障的過渡電阻范圍可設置為[0,1 000]Ω[17],兩相相間短路故障的過渡電阻可設置更小的取值范圍,三相短路故障可不設置。

本算例以三相短路電流為例,通過仿真獲得10 000組樣本。由于仿真建模時IIDG采用較為精確的模型,因此仿真速度很慢,獲取10 000組樣本耗時約10 d。按照8∶2 的比例將樣本分為訓練集和測試集。根據式(8),該配電網共有316 個輸入特征,根據式(9),則有91個輸入特征,輸出結果為80個。

3.2 多輸出模型有效性驗證

文獻[11]采用XGBoost和LightGBM 模型進行預測,并已與物理建模方法進行了對比,結果表明機器學習計算短路電流的準確性。但文獻[11]為單輸出模型,若要輸出圖1虛線框中測量點①、②、③的短路電流,則需要建立3個模型,而本文所提方法僅需1個模型即可輸出上述3處測量點的短路電流。表1為不同方法預測節點834 發生三相短路故障時C 相短路電流的結果與模型的性能指標。其中,XGBoost(單輸出)模型的數據來源于文獻[11],為測量點①的電流,MTRS-XGBoost(多輸出)和MTRS-LightGBM(多輸出)模型的結果為測量點①、②、③處的預測結果。因本文與文獻[11]的系統等值阻抗、負荷、IIDG容量等均在相同范圍內隨機生成,故而運行方式差異不會太大但無法保證完全一致,因此選擇了實際值較為接近的場景來進行性能比較。

表1 單輸出模型與多輸出模型結果Table 1 Result of single-output model and multi-output model

由表1 可知,XGBoost(單輸出)模型只能預測單個測量點的短路電流,λAPE在可接受范圍內;而多輸出模型預測的短路電流,因考慮各輸出間的相關性,λAPE較單輸出模型更低。以測量點①為例,MTRSXGBoost(多輸出)和MTRS-LightGBM(多輸出)模型比XGBoost(單輸出)模型的λAPE小,分別降低了80.7%和73.4%。此外,對于不同測量點的短路電流,使用不同基學習器模型時短路電流誤差情況略有不同。綜上所示,采用多輸出模型計算含IIDG 配電網短路電流不僅有效,且準確性更高。

此外,某運行方式下測量點①、②、③的三相短路電流相量如附錄A圖A1所示。由圖可見,IIDG 的接入使得各支路穩態電流的幅值不再顯式地滿足基爾霍夫電流定律,但相量依然滿足基爾霍夫電流定律。

3.3 多輸出模型性能對比

下面對基于問題轉化方法和算法適應方法的2類多輸出模型的性能進行對比,同時對比式(8)和式(9)所示的輸入特征對計算性能的影響,結果如表2所示。表中:訓練時間為完成模型訓練的時間;測試時間為完成2000組測試的總耗時;加粗數值表示每種多輸出模型下的性能最優項。對表2 進行分析可得到下列結論。

表2 不同多輸出模型、不同基學習器、不同輸入特征的性能對比Table 2 Performance comparison of different multi-output models,base learners and input characteristics

1)6 種模型和2 種輸入特征均取得了較好的性能;δMAX、δMEAN、λMAX和λMEAN均在可接受范圍內;最長測試時間為11.750 62 s,即最長平均測試時間約為5.875 ms,而單輸出模型的平均預測時間為0.049 s[11],因此多輸出模型預測速度快,約為單輸出模型的8.4倍。

2)對比同一種多輸出模型使用不同輸入特征的計算性能可知,6 種模型選用式(8)所示的輸入特征時,其MAPE 指標優于選用式(9)所示的輸入特征時的MAPE 指標;采用問題轉化方法時,多輸出模型的MAE 比較接近,而采用基于神經網絡的算法適應方法時,多輸出模型的MAE 有一定差異,且選用式(8)、(9)所示的輸入特征均有可能獲得相對較優的結果。因此,2 類輸入特征對多輸出模型的準確度影響不大,在尚未普及μPMU 配置時,可以使用式(9)所示的輸入特征。

3)在問題轉化方法中,2 種基學習器的性能在準確度方面差異并不大,但是在速度卻有快慢之分。采用MTRS 方法時,在相同的輸入特征下,LightGBM的速度快于XGBoost,這證明在輸入特征變動不頻繁的情況下,LightGBM 速度更快,但也因為EFB 特性導致與XGBoost 的速度并沒有明顯的區分。采用RC 方法時,EFB 特性反而成為了LightGBM 提升速度的障礙。

4)對比同一種輸入特征下不同多輸出模型的性能可以發現:2 類問題轉化方法中,XGBoost 基學習器的誤差略小于LightGBM 基學習器,計算速度方面更多地取決于所采用的問題轉化方法,應用時RC方法明顯快于MTRS 方法;算法適應方法中,ANN 比CNN 有更低的誤差,而計算速度不如CNN,上述結論與理論分析一致。算法適應方法的準確度不如問題轉化方法,但應用速度快于問題轉化方法。綜上,針對含IIDG 配電網短路電流計算問題,采用問題轉化方法的準確性都非常高,計算速度最快的是選用XGBoost 基學習器的RC 方法。而神經網絡多輸出方法適用于對準確性要求一般,而對應用速度要求極高的場合。此外,選取測量點①、②、③處的A 相電流的真實值與不同算法輸出進行對比,結果見附錄A圖A2。由圖可見,故障線路上基于ANN和CNN的短路電流預測值偏差較為明顯,其他差異不太明顯。

3.4 模型抗干擾能力

在實際應用時,數據采集過程可能出現部分特征數據丟失的情況。選用表2 中綜合性能較優的方法,即分別使用MTRS、RC 拼合的XGBoost 算法和ANN 多輸出算法,輸入特征選用式(8),數據丟失以隨機置0 的方式表示。3 種模型在不同數據丟失比例rloss下的性能表現如表3 和附錄A 圖A3 所示,表3中加粗數值為較小誤差數據。

表3 表明:在丟失1%的數據量時,3 種模型的誤差均較小,當丟失數據量大于等于3%時,ANN 多輸出模型的誤差增大,準確性明顯降低;使用問題轉化的多輸出模型在丟失5%、10%的數據量情況下仍能保持較優的性能,其中MTRS 模型誤差略優于RC模型。

表3 不同數據丟失比例下多輸出模型的性能對比Table 3 Performance comparison among multi-output models with different rates of data loss

圖A3 為數據丟失比例高于10%的情況下,不同多輸出模型的λMEAN。由圖可見:當數據丟失比例在30%以下時,MTRS-XGBoost 模型和RC-XGBoost模型的λMEAN均較??;當數據丟失比例超過30%后,MTRS-XGBoost 模型和RC-XGBoost 模型的λMEAN快速上升;當數據丟失比例為60%時,MTRS-XGBoost模型和RC-XGBoost模型λMEAN與ANN 多輸出模型相近。由此可知,當數據少量丟失的情況下,MTRS 和RC 方法均具有較好的抗干擾能力,當數據丟失比例較高時,2 種方法的誤差隨數據丟失比例的升高而快速增大。

4 結論

含IIDG 配電網發生短路故障時,IIDG 的輸出具有較強非線性特征,使得短路電流的迭代計算耗時增長,而采用簡化模型又影響計算的準確性。在IIDG 高占比接入配電網的發展趨勢下,本文提出了數據驅動的含IIDG 配電網短路電流計算多輸出模型,解決了計算準確性與計算速度之間的矛盾,得到如下結論。

1)針對含IIDG 配電網短路電流計算,選用多輸出模型,不僅解決了單輸出模型的模型數量過多的問題,且具有更高的準確性和更快的計算速度,例如采用MTRS-XGBoost 模型時λAPE降低了80.7%,最長測試時間下也比單輸出模型快8.4倍,因此本文所提出多輸出模型能滿足在線應用要求。

2)無論在常規的配電網數據采集條件,還是配置有μPMU 的條件下,都可以滿足所提多輸出模型的輸入特征條件;即使在輸入數據丟失30%的情況下,λMEAN仍能達到小于5%的要求。

3)相比較而言,基于MTRS 方法的多輸出模型具有準確性高和抗干擾能力強的優點,雖然計算速度不及其他多輸出模型,但也遠快于單輸出模型和物理建模計算速度。RC 方法的準確性和抗干擾性能力與MTRS方法相近,但計算速度較MTRS方法快1 個數量級,是較為穩定的多輸出方法。使用神經網絡的算法適應方法具有最快的應用速度,但是其準確度相對略低,適用于對計算快速性要求極高的場景。

附錄見本刊網絡版(http://www.epae.cn)。

猜你喜歡
測量點短路準確性
飛機部件數字化調姿定位測量點的優選與構造算法
淺談如何提高建筑安裝工程預算的準確性
理解語境與名句的關系,提高默寫的準確性
淺析沖壓件測量點的規劃
熱電偶應用與相關問題研究
基于CAD模型的三坐標測量機測量點分布規劃
為橋梁領域的示值準確性護航
短路學校
影響紫外在線監測系統準確性因子分析
短路學校
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合