?

基于改進自注意力機制生成對抗網絡的智能電網GPS欺騙攻擊防御方法

2021-11-20 08:33李元誠楊珊珊
電力自動化設備 2021年11期
關鍵詞:防御機制注意力向量

李元誠,楊珊珊

(華北電力大學 控制與計算機工程學院,北京 102206)

0 引言

現代電力系統的規模不斷擴大,互聯程度不斷加強,結構和運行方式日趨復雜,對狀態估計[1]提出了更高的要求,傳統的數據采集與監視控制SCADA(Supervisory Control And Data Acquisition)系統所采集量測值的精度已不能滿足狀態估計的要求。相量測量裝置PMU(Phasor Measurement Unit)[2]能夠對節點電壓相量和支路電流相量等狀態量進行高精度的同步量測,有效提高智能電網狀態估計的精度[3]。但PMU成本高,經常采用基于PMU和SCADA系統混合量測的電力系統狀態估計方法[4]。

PMU 在提高智能電網狀態估計精度的同時,也帶來了一系列的安全問題。PMU 的全球定位系統(GPS)接收機以射頻方式接收來自不同衛星的GPS信號。GPS 信號包括未加密的C/A 碼和加密的P(Y)碼。PMU 接收到的GPS 信號為未加密的C/A碼[5]。C/A 碼的碼結構是公開的,很容易被攻擊者利用,通過偽造GPS信號欺騙接收機,所以未加密的GPS 信號容易受到欺騙攻擊。研究人員已經對GPS欺騙攻擊GSA(GPS Spoofing Attack)進行了現場測試,驗證了其可行性[6]。此外,欺騙者可以成功地修改PMU 量測值的時間戳。當GPS 信號受到欺騙時,相應的PMU 量測值會變得不可靠,因此必須在檢測到欺騙后立即通過刪除或校正等方式對量測值進行修正,否則,基于欺騙數據的狀態估計結果將是錯誤的,并可能誤導控制中心發起不必要的、可能破壞穩定的補救控制措施。例如,基于欺騙數據的狀態估計可能導致墨西哥當前運行系統的控制方案自動使發電機跳閘,這會導致整個電網發生級聯故障,甚至崩潰[7]。因此,如何高效地對GSA 進行防御,減輕其帶來的危害,對于保障智能電網的安全運行具有重要意義。

針對GSA 的防御機制主要分為2 類:一類是在GPS 接收器接收GPS 信號階段進行防御;另一類是在PMU 測量數據之后及進行狀態估計之前的階段進行防御。第一類防御機制通過檢測導航數據,區分真實的GPS 信號和偽GPS 信號,通過恢復偽GPS信號對GSA 進行防御。該類方法可以分為信號處理方法、天線防御方法、基于與其他時間源相關性的防御方法、加密策略四大類。信號處理方法需要從接收信號中監測導航信息的功率、質量,提取出導航信息的具體特征用于識別不相關數據,利用信號中的數學關系能夠區分真實的GPS 信號和偽GPS 信號[8]。天線防御方法大多采用到達角分辨技術[9],但通常需要高電平的附加硬件?;谂c其他衛星系統權威機構的相關性的防御方法僅對特定的GPS信號[10]有用。加密策略較昂貴,且需要修改GPS 信號結構,主要用于軍事版本的GPS 信號[11]。第二類防御機制是在PMU 測量數據之后,通過各種算法提取出量測數據的具體特征用于識別異常數據,對被攻擊的PMU 量測值進行校正。文獻[5]提出了一種基于探測技術的GSA 識別算法,文獻[12]提出了一種基于廣義似然比檢驗的校正算法,但都只適用于單一GSA 的情況。文獻[9]從物理角度提出了一種針對多GSA 的檢測機制,需在PMU 的現有接收器附近安裝另一臺商用GPS接收器。文獻[13]考慮多GSA的情況,將狀態估計、攻擊重構問題轉化為雙線性最小二乘問題,并利用交替最小化算法進行求解,但所提方法僅適用于系統可觀測的環境,并不適用基于PMU 和SCADA 系統混合量測的電力系統環境。文獻[14]提出了一種可以同時估計線參數和修正量測值的方法,但系統需具備一個預先校準的PMU。文獻[15]提出了一種多層的多接收機結構,可以增強GPS 的定時性能,以抵抗干擾、欺騙和接收機誤差,但會增加額外的硬件成本。

生成對抗網絡(GAN)[16]已被應用于網絡攻擊的防御策略研究中,但大多GAN 方法在某些設定條件下仍存在樣本生成質量低或者不能收斂等問題。文獻[17]將自注意力機制引入GAN 的框架,提出了自注意力機制生成對抗網絡SAGAN(Self-Attention Generative Adversarial Network)。自注意力機制模塊在建模長期依賴方面很有效,有利于提高生成樣本的質量。另外,SAGAN 將譜歸一化應用于生成器,解決了GAN 訓練穩定性的問題,且生成器和判別器更新規則TTUR(Two-Timescale Update Rule)加速了正則化判別器的訓練。

本文首先在SAGAN 結構[17]中引入基于門控循環單元GRU(Gated Recurrent Unit)的時間注意力模塊,提出了一種改進SAGAN 模型;然后,基于改進SAGAN,實現對引入深度學習參數的新網絡-物理模型的訓練。改進SAGAN 通過學習歷史量測值,提取數據的時空特征,利用判別器對欺騙攻擊進行檢測,利用生成器實現對欺騙數據的恢復,最終得到一個基于改進SAGAN的主動防御模型。與現有研究相比,本文所提方法無需任何變電站級別的硬件增強,極易應用于實時場景或現場研究。此外,本文所提方法不僅能實現對GSA 的檢測,還能對欺騙數據進行校正。隨著現代化電網的數據量越來越大,電網拓撲結構越來越復雜,現有研究提出的大多數據校正算法可能無法保證實時性,但算例仿真結果證明本文所提防御方法能很好地保證防御算法的實時性。

1 問題描述

1.1 GSA過程

GSA 的原理為:攻擊者的信號發射器發射與衛星信號結構相同/相似而功率更強的信號,使PMU中的GPS接收機誤以為其是真實信號而進行搜索捕獲。衛星導航系統的基本原理是無線電測距定位,近地面用戶通過接收4 顆以上衛星的信號,計算偽距后通過球面交匯完成定位。所以GSA 的本質是通過發射信號誤導目標用戶得到錯誤的傳播時延、偽距,致使目標用戶定位到欺騙點上。若欺騙干擾源能靈活控制偽距的變化,則可實現任意位置欺騙。

為了欺騙GPS 接收機,需要誤導GPS 接收機獲取偽GPS 信號。由于捕獲階段是通過在載波頻率-碼相二維空間中搜索最高相關峰值(直觀而言,具有較高信噪比(SNR)的信號將具有較高的相關峰值)實現的,故存在一種兩步法欺騙策略:在第一步中,欺騙者發射某些干擾,導致GPS接收機失去跟蹤;在第二步中,當GPS 接收機進行捕獲時,欺騙者啟動偽GPS信號。由于偽GPS 信號具有更高的SNR,GPS 接收機將因搜索到較高的相關峰值而跟蹤偽GPS信號。

GSA的過程如下:攻擊者通過學習真實的GPS信號,判斷給定時間區域內攻擊目標附近的軌道衛星,然后利用公開數據庫中的公式,偽造不同衛星的C/A碼,在攻擊目標附近廣播與真實信號相同的C/A碼值信息;當攻擊目標在捕獲階段的載波頻率-碼相二維空間中進行粗略掃描以搜索功率較大的GPS信號時,攻擊者逐漸增加虛假信號的功率,致使GPS接收機鎖定偽GPS 信號,之后攻擊者可以慢慢改變偽GPS 信號的碼相,接收機會調整其信號發生器與偽信號對齊,使碼相偏離真實信號,最終將真實信號當作噪聲處理。碼相是計算傳播時間和時間偏置的關鍵,因此GSA 會通過隨機移動GPS 信號中的相角來破壞接收機與系統時間之間的時間同步,最終使接收機估計得到錯誤的衛星位置和時鐘偏置量,使PMU計算得到錯誤的相角,之后能量管理系統(EMS)中的狀態估計器利用被篡改的量測值估計得到不正確的系統狀態,從而給電網帶來嚴重威脅[18]。

1.2 受GSA影響的PMU量測值

假設在有N條母線的電力系統中安裝了p臺PMU,PMUk提供的同步數據表示如下:

式中:mk為PMUk提供的量測值向量,其元素mki(i=1,2,…,n;n為量測值數量)為PMUk提供的序號為ki的量測值;s為電網狀態向量;Sj(j=1,2,…,N)為母線j的電壓相角;Ak由系統狀態和PMUk提供的量測值決定,可以基于PMU 的位置、網絡拓撲和傳輸線參數獲得;ek為PMUk的測量噪聲向量,不失一般性,假設ek服從相同的獨立分布(方差為σ2的復高斯分布)。

將p臺PMU提供的量測值進行疊加,可得:

式中:m、A、e可以分別由對應于不同PMU 的子塊mk、Ak、ek適當地構造得到,在不失一般性的前提下,假設所有PMU的測量噪聲服從相同的分布。

假設一個GSA出現在PMUk上,且相量偏移量為θspf,根據GSA 對同步相量數據的影響特征,被欺騙的同步相量數據可表示為:

式中:mspf為被欺騙的同步相量數據;G為攻擊矩陣;I(kk=1,2,…,p)為單位矩陣,其大小由PMUk提供的量測值數量決定(即與mk同維度)。

2 GSA防御機制

為了防御GSA,本文提出了一種改進網絡-物理模型,該模型根據歷史量測數據計算當前的量測值,并基于改進SAGAN實現模型訓練。

2.1 改進網絡-物理模型

首先通過回歸算法對量測值向量進行預處理。量測值的非線性回歸模型如式(6)所示。

式中:Mt為t時刻的量測值向量;Vi、Vj分別為母線i、母線j的電壓幅值向量;Pij、Qij分別為線路ij的有功潮流、無功潮流向量;θij為母線i、母線j間的電壓相角差向量;gij、bij分別為線路ij的電導、電納向量;αt、βt、ηt、γt、δt為t時刻的測量系數;l(Pij;Qij)為有功潮流與無功潮流間的非線性關系函數;f(gij;bij;θij)為電導、電納、電壓相角差之間的非線性關系函數;et為t時刻的測量噪聲向量。

基于式(6),可以得到計算t時刻量測值的改進網絡-物理模型如式(7)所示。

式中:Mc(t)為改進網絡-物理模型計算所得t時刻的量測值向量;Vi(t-1)、Vj(t-1)、Pij(t-1)、Qij(t-1)、gij(t-1)、bij(t-1)、θij(t-1)為歷史量測值向量;α、β、η、γ、δ為可以通過深度神經網絡模型學習得到的模型參數值;ut為t時刻深度神經網絡模型預測所得量測值與實際量測值之間的偏差向量。

結合所提改進網絡-物理模型,利用t-1時刻采集的量測值對神經網絡模型進行訓練,可以計算得到t時刻的量測值。

2.2 改進SAGAN模型

SAGAN 中的生成器和判別器網絡采用卷積結構,結合自注意力機制可以有效提取輸入數據的空間特征,但并不擅長提取數據中的周期性變化規律。在現代化電網中,拓撲信息越來越復雜,數據量也越來越大,量測值不僅具有空間相關性,也具有時間相關性,而GSA 會篡改量測值,導致被篡改后的量測值具有與正常數據不同的時空特征。所以同時提取量測值的時空特征,學習正常數據和被欺騙數據的不同之處,可以有效提高修正數據的精度。

循環神經網絡(RNN)可以有效地對動態時間序列數據進行建模,學習序列的周期性變化規律。而GRU 是RNN 的一種變體,其訓練時間短,收斂速度快,可以避免梯度消失或梯度爆炸問題。

將SAGAN 生成器和判別器中的空間自注意力模塊提取所得的特征圖輸入GRU 中進行動態時間建模,同時提取數據的時間和空間特征,可以有效提高判別器的判別準確率和生成器生成樣本的質量,所以本文在SAGAN 的生成器和判別器網絡的最后一層卷積層之后分別融入基于GRU 的時間注意力模塊,提出了一種改進SAGAN 模型,通過結合多種結構有效學習輸入數據中豐富的特征和規律。

本文所提改進SAGAN 中空間自注意力模塊和時間注意力模塊的具體敘述如下。

1)空間自注意力模塊。

將前一隱藏層輸出的數據特征x∈RC×M(C為通道數,M為前一隱藏層特征的位置數量)變換為2 個特征空間f、g來計算注意力,其中f(x)=Wfx,g(x)=Wgx,Wf∈RCˉ×C、Wg∈RCˉ×C(Cˉ為通道數,且有Cˉ=C/K,K=1,2,4,8)為學習權重矩陣,則模型在合成第j個區域時對第i個區域的關注程度βj,i可表示為:

式中:sij=fT(xi)g(xj),xi、xj分別為第i、j個區域的數據特征。自注意力層的輸出為o=[o1,o2,…,oj,…,oM]∈RC×M,其第j列元素oj如式(9)所示。

式中:Wh∈RCˉ×C、Wv∈RCˉ×C為學習權重矩陣。為了提高內存效率,本文選取K=8進行測試。

此外,將自注意力層的輸出與比例參數相乘,并重新添加輸入的特征圖,得到最終空間自注意力模塊的輸出為:

式中:γs為可學習的標量,初始值為0。2)時間注意力模塊。

將空間自注意力模塊輸出的特征圖構造為時間序列形式輸入GRU 進行動態時間建模,并引入注意力機制通過映射加權和學習參數矩陣賦予GRU 隱含狀態不同的權重,減少歷史信息的丟失并加強重要信息的影響。其中單個GRU 結構由更新門和重置門組成,具體公式如下:

式中:zt為更新門,用于更新行為時的邏輯門;rt為重置門,用于決定候選行為時是否要放棄之前的行為ht;xt為t時刻GRU 的輸入;h?t為t時刻的候選行為,接收{xt,ht-1};ht為t時刻GRU 的隱藏層輸出,接收{ht-1,h?t};Wz、Uz、Wt、Ut、W、U為權重矩陣;“°”為哈達碼積;σ(?)、tanh(?)為激活函數。

將GRU層的輸出記為H。注意力機制層的輸入為經過GRU 層激活處理的輸出向量H,根據權重分配原則計算不同特征向量對應的概率,不斷更新迭代得到的較優權重參數矩陣。注意力機制層的權重系數計算公式為:

式中:Gt為t時刻由GRU 層輸出向量ht所決定的注意力概率分布值;u、w為權重系數;B為偏置系數;at為注意力機制層在t時刻的輸出。

在改進SAGAN 模型中,通過最小化鉸鏈式的對抗性損失函數對模型交替訓練,判別器的損失函數LD、生成器的損失函數LG分別見式(20)和式(21)。

式中:E[?]為求期望值;x~pdata表示量測值向量數據x服從真實數據分布;D(x)為判別器正確判斷的概率;G(z)為生成器生成的數據;D(G(z))為判別器對生成器生成數據的判別結果;z為噪聲向量。

2.3 GSA防御模型

將改進SAGAN 模型應用到GSA 防御中,設計GSA 防御模型。防御模型分為數據采集與控制、數據檢測與防御2 個模塊,其中,數據采集與控制模塊進行數據采集、狀態估計、初步檢測壞數據等操作,并將含有欺騙數據的量測數據輸入數據檢測與防御模塊;數據檢測與防御模塊包括訓練好的生成器和判別器2 個組件,主要采用的防御機制見2.4 節。本文所提GSA防御模型如圖1所示。

圖1 GSA防御模型Fig.1 Defense model of GSA

GSA 防御模型的主要工作原理可以概括如下:攻擊者對PMU 的GPS 接收機進行欺騙,引入錯誤的時間戳,從而篡改同步相量數據;PMU 和SCADA 系統將采集的數據傳送到狀態估計器進行狀態估計和壞數據檢測;經過壞數據檢測后,將量測數據輸入所提算法中;利用所提算法對量測數據進行檢測,并對檢測所得欺騙數據進行校正。數據檢測與防御模塊的實現過程如下:

1)進行GSA檢測,檢測數據是否為欺騙數據;

2)分離欺騙數據和正常數據;

3)刪除欺騙數據并將剩余正常數據輸入生成器;

4)利用生成器從正常數據中提取特征;

5)生成器結合網絡-物理模型,根據提取的特征計算得到與損失數據特征盡可能一樣的數據;

6)將計算得到的數據加入剩余正常數據中,合成補全后的數據,將其發送給狀態估計器。

2.4 防御機制

根據所提改進SAGAN 模型設計GSA防御機制,如圖2所示。圖中,M為原始量測值數據;Mr為正常量測值數據;Mc為生成器計算得到的數據;G(Mc)表示將一個向量Mr輸入生成器G 并輸出量測值Mc的構造過程。

圖2 GSA防御機制Fig.2 Defense mechanism of GSA

本文所提防御機制的基本思想如下:①將原始量測值數據M輸入判別器D,判別器D 的輸出為概率值D1和D2,其中D1為生成器G計算得到的數據屬于正常數據的概率,D2為該原始量測值數據屬于正常數據的概率;②分離被欺騙的量測數據mspf與正常量測值數據Mr;③將Mr輸入生成器G,生成器G計算并輸出Mc,將Mc作為輸入發送給判別器D。改進SAGAN 中的判別器D 和生成器G 被交替訓練,如果D1很大,則說明生成器G 生成樣本的質量已經符合要求;否則,說明質量未達到要求,需要繼續訓練。譜歸一化可以防止參數幅度增大,避免異常梯度。使用譜歸一化的生成器和判別器使模型在每次生成器更新時進行較少的判別器更新次數,從而大幅縮短了訓練時間,同時也提高了訓練的穩定性。

判別器D 的訓練目的是為了使EM[D(M)]最大化,以此實現從量測值數據M中提取更好的特征進行建模,增大正確判斷的概率D(M),即最大化甄別哪些向量值屬于真實數據分布的概率。生成器G 的訓練的目的是為了使EMc[D(Mc)]最大化,從而使生成器計算得到的數據Mc能接近真實的未被欺騙的量測數據,即最大化生成器輸出的向量值被判別器判斷為來自真實數據分布的概率。

本文所提改進SAGAN 的具體訓練過程見附錄A。在改進SAGAN 中,生成器和判別器以不同的學習率同時進行訓練,具體步驟如下。

1)對原始量測數據集{mi}進行預處理,將n個量測值向量處理為a×b階矩陣M,如式(22)所示。將矩陣M作為改進SAGAN的輸入。

2)將經過預處理的矩陣M和生成器生成的數據輸入判別器的空間自注意力模塊,利用其卷積層和自注意力層對輸入數據進行高維特征提取,輸出特征圖。

3)將空間自注意力模塊輸出的特征圖構造為時間序列形式,作為時間注意力模塊的輸入。

4)將步驟2)和步驟3)所得結果輸入全連接層,分離正常數據和被欺騙數據,刪除欺騙數據。

5)將正常量測數據作為先驗的輸入噪聲變量輸入生成器,經過一個全連接層,將其重塑為2維矩陣,然后將2維矩陣輸入生成器的空間自注意力模塊對輸入數據進行高維特征提取,輸出正常數據的特征圖。

6)將空間自注意力模塊輸出的特征圖構造為時間序列形式,作為時間注意力模塊的輸入,并提取數據的時間特征,得到生成數據。同時,將生成的數據輸入判別器,利用判別器進行判別。

7)經過生成器和判別器不斷地進行博弈,最終輸出最接近于正常數據特征的數據。

3 算例分析

3.1 參數設置

為了驗證所提防御方法的有效性,以IEEE 14節點、IEEE 118節點測試系統為例,采用文獻[19]中的系統拓撲結構、參數進行仿真設置。每個系統對應不同的PMU位置配置文件,具體見附錄B表B1。

在每個系統中,隨機選擇1 臺或多臺PMU,對每臺PMU隨機選擇相量偏移量θspf進行GSA,并在傳輸到狀態估計器之前修改其量測值。本文考慮對1 臺PMU 進行攻擊的情況,采取文獻[20]中的欺騙干擾策略,被欺騙者修改的PMU量測值以0.8(°)/s的速度偏離真實值,在2 min內突破PMU的IEEE C37.118標準[9]。本文所提方法通過學習量測值來恢復被欺騙的量測值數據,所以同樣適用于針對多臺PMU 進行攻擊的情況。設改進SAGAN 判別器、生成器的學習速率分別為0.0004、0.0001。

3.2 GSA下PMU量測值變化測試

為了驗證所提防御機制的性能,在進行防御機制性能測試之前,首先對受GSA 的PMU 量測值變化情況進行研究。由于GSA 只會對電壓相角值造成影響,而對電壓幅值的影響不大,本節只研究電壓相角值的變化情況。實驗分別在2 個測試系統中進行。設IEEE 14 節點系統中的PMU1(即母線2 處的PMU)遭受相量偏移量為θspf的GSA(將其記為GSA(1,θspf)),具體實驗結果如圖3 所示。針對IEEE 118節點系統的GSA結果見附錄B圖B1。

圖3 IEEE 14節點系統遭受GSA時PMU量測值的變化曲線Fig.3 Change curves of PMU measurement values when IEEE 14-bus system suffers GSA

由圖3 和附錄B 圖B1 可以看出,在2 個測試系統中,攻擊者在50 s 后開始攻擊,50 s 之前的數據為正常數據,50 s 后的數據為受欺騙的數據,IEEE 14節點系統中的PMU2受到欺騙,IEEE 118 節點系統中的PMU9受到欺騙,欺騙在2 min 內使相角偏離真值10°。因此,GSA的相量偏移量會對系統母線的電壓相角值產生一定的影響。

3.3 防御機制性能測試

3.2 節的測試結果表明,GSA 會對系統的電壓相角值造成一定的影響。在防御機制中,生成器需要從正常數據中提取數據特征,計算得到與正常數據特征盡可能一樣的數據。生成器在生成數據之前需要進行迭代訓練,迭代的次數不同,則訓練的效果不同。測試中設置不同的迭代次數,比較生成數據與未受攻擊的正常數據之間的差別。在IEEE 14 節點、IEEE 118節點系統中進行防御測試,驗證所提防御機制的性能,結果分別見圖4和附錄B圖B2。

圖4 迭代次數不同時的電壓相角值(IEEE 14節點系統)Fig.4 Voltage phase angle values with different iteration times(IEEE 14-bus system)

由圖4 可看出,在IEEE 14 節點系統中,當訓練迭代次數為10 次時,生成器生成的數據略高于未受攻擊的正常量測值;當迭代次數增加到100 次時,生成器生成的數據幾乎接近正常量測值。由圖B2 可看出,在IEEE 118 節點系統中,當訓練迭代次數超過100 次后,生成器生成的數據與正常量測值接近。上述結果證明了本文所提防御機制在不同測試系統中的可行性。

關于測試算法的計算時間,IEEE 14節點系統和IEEE 118節點系統的計算時間分別為0.047、1.092 s,可見本文所提防御機制具有較短的延遲,從而保證了實時性防御。

為了進一步驗證本文所提防御機制的有效性,設置生成器的迭代次數為100次,在IEEE 14節點系統和IEEE 118 節點系統中進行不同算法的對比測試,結果分別如圖5和附錄B圖B3所示。由圖可知,相比于交替最小化算法、GAN 算法和SAGAN 算法,改進SAGAN 算法生成的數據更接近于正常量測值,從而證明了所提改進SAGAN算法的有效性。

圖5 不同算法所得電壓相角值對比(IEEE 14節點系統)Fig.5 Comparison of voltage phase angle values obtained by different algorithms(IEEE 14-bus system)

不同算法的計算時間對比如表1 所示。由表可知,相較于其他防御算法,本文所提算法的計算時間最短,能更好地保證實時防御。

表1 不同算法的計算時間對比Table 1 Comparison of computation time among different algorithms

4 結論

為了檢測和防御GSA,本文提出了一種新的網絡-物理模型,該模型在原模型的基礎上引入了一個深度學習參數,具有更高的計算準確率?;诖司W絡-物理模型,提出了基于改進SAGAN 的GSA 檢測和防御機制。在該機制中,利用改進SAGAN 的空間自注意力模塊對歷史量測數據進行特征提取,基于GRU 的時間注意力模塊對所提特征進行動態時序建模,不斷地通過生成器、判別器對量測數據進行學習,通過數據生成、對比和替換實現對GSA 的檢測和防御。將本文所提防御算法與交替最小化算法、GAN 算法進行對比,測試結果表明本文所提防御方法生成得到的數據更接近正常數據,防御效果更好。

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

猜你喜歡
防御機制注意力向量
向量的分解
讓注意力“飛”回來
防御機制在醫學生抑郁的性別差異中的中介作用*
聚焦“向量與三角”創新題
“揚眼”APP:讓注意力“變現”
A Beautiful Way Of Looking At Things
向量垂直在解析幾何中的應用
向量五種“變身” 玩轉圓錐曲線
軀體形式障礙患者治療前后防御機制的對照觀察
大學生防御機制及人際信任的相關研究
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合