?

起伏目標相干檢測性能分析

2021-04-19 12:38黃子芹
信號處理 2021年4期
關鍵詞:虛警信噪比脈沖

王 杰 劉 劍 黃子芹 肖 楠

(空軍工程大學信息與導航學院, 陜西西安 710077)

1 引言

雷達的目標檢測主要解決的是目標有無的判決問題,這始終是衡量雷達性能高低的最重要因素之一[1]。隨著隱身技術的日趨成熟,飛機、導彈等目標的雷達橫截面(Radar Cross Section, RCS)減小了一到兩個數量級。目標RCS變小使得目標回波能量減小,雷達的目標檢測性能隨之大為降低[2]?,F代戰爭中,隱身飛機外、巡航導彈和無人機等高速機動的低可探測目標(Low observable targets),除RCS極低外,其RCS也是起伏的,這就對雷達的檢測性能提出了嚴峻的挑戰[3]。

要提高雷達對于隱身飛機等目標的檢測性能,本質上就是要在極低的信噪比(Signal to Noise Ratio, SNR)的條件下,盡可能提高雷達目標回波能量[4]。針對隱身飛機等目標RCS極低的情況,通常采用脈沖積累的方法。脈沖積累又分為相參積累和非相參積累,非相參積累未利用回波信號的相位信息,在面對信噪比極端低下的目標時,往往由于積累增益不夠高而失效[5]。而相參積累能將信號相干相加以獲得相參積累增益,從而提高檢測性能[6]。但相參積累需要使用所有積累脈沖的相位信息,所以必須明晰積累脈沖間的相位變化,保證相位一致性[7]。文獻[8]對于無線通信系統中的載波同步做了大量研究,基于雷達與無線通信的系統結構高度相似性,可在雷達系統中實現通信技術的優勢利用[9],以保持回波信號的相干性。對于雷達,文獻[10]提出一種新的補償方法,計算復雜度較低,能補償回波信號相位的變化,補償后脈沖間信號的相位一致,可以相干相加。雖然現有雷達一般都采用相參積累技術[11-12],但這種相參積累是對信號的包絡進行積累,積累后的信號再進行非相干檢測,仍舊損失了信號的部分相位信息。

對于目標的RCS起伏特性,通常使用起伏目標模型來估計目標起伏的影響并進行數學上的分析。目前的研究,如文獻[4]和[13]等,大都是針對非起伏目標進行相參積累。但實際上,目標往往都具有起伏特性,針對隱身飛機等起伏目標的相參積累研究,其實際應用價值也更高[14]。大多數起伏目標的統計特性由卡方分布表述,現有模型中,目標統計特性服從卡方分布且在理論上比較成熟的有4個常用的Swerling[15]模型。該模型對于一些導彈、飛機等目標的擬合性很好,常被用于分析雷達對目標的檢測性能[16]。非相干檢測可以運用于Swerling的四種起伏模型,但是相干檢測只能用于Swerling Ⅰ、Ⅲ型,這是因為Swerling Ⅱ、Ⅳ型起伏目標脈沖與脈沖間的目標幅度不相關(快起伏),無法保持相位相干性[17]。但根據文獻[18]的研究表明,對于一定起伏頻率范圍內的快起伏目標,長駐留信號能夠很好地起到平滑作用,減小目標起伏帶來的影響[19],相干檢測或可應用到快起伏模型中。

綜上,對抗隱身飛機等目標的低RCS特性需要相參積累,而對抗其RCS起伏特性需要高效的起伏目標檢測方法。若能將相干檢測應用于起伏目標,并對脈沖進行積累,就能充分利用回波信號的相位信息,提高雷達對這類目標的檢測概率,有效地應對隱身飛機等起伏目標的挑戰。本文在分析非相干檢測性能的基礎上,從檢測原理出發,提出將相干檢測應用于起伏目標以提高雷達檢測性能,基于重積分方程,推導了Swerling Ⅰ、Ⅲ型起伏目標相干檢測概率計算公式。并通過仿真驗證,對比分析了兩種檢測方式的檢測性能和起伏損失。除Swerling模型外,常用的目標起伏模型還有對數正態分布模型和萊斯分布模型等。文獻[20]和[21]還針對隱身飛機等目標的RCS起伏統計特性進行了專門研究。選用的起伏模型越精確,越能反映隱身飛機等目標的RCS特性[22],應用相干檢測后也越能提高雷達對于隱身飛機等目標的檢測概率。本文安排如下:第2部分,給出了起伏目標模型,第3部分,分析非相干檢測原理與性能,第4部分,分析相干檢測原理及其非起伏目標性能,第5部分給出起伏目標的相干檢測公式,第6部分通過仿真,驗證了起伏目標相干檢測公式的準確性,并對比非相干檢測進行分析,最后給出結論。

2 起伏目標模型

雷達目標通常由多個散射體組成,通過對來自單個散射體,既包含幅度信息又包含相位信息的接收信號采用相干求和的方法,就可以計算出目標的回波幅度。對于靜止目標,回波幅度通常是固定值,但運動目標相對雷達視線的姿態角會不斷變化,使散射體子矢量合成時,各自的相對相位隨機變化,從而產生回波幅度起伏。

通過特定的波形設計和對回波信號幅度與相位的處理、分析和變換,可以得到目標RCS及其起伏統計模型,這可以表征雷達目標的固有特征,所以雷達回波的起伏總是與雷達目標的RCS相聯系。由于雷達需要探測的目標十分復雜而且多種多樣,很難準確地得到各種目標截面積的概率分布和相關函數,通常是用一個接近而又合理的模型來估計目標起伏的影響并進行數學上的分析。

文獻[23]提出χ2分布模型且指出,若其參數設置得當,可用于近似幾乎其他所有模型。χ2分布模型為:

(1)

經典的Swerling模型把典型的目標起伏分為四種類型,且掃描至掃描間都是完全不相關的,其中Swerling Ⅰ、Ⅱ型和Ⅲ、Ⅳ型分別有相同的概率密度函數。χ2分布模型包含了傳統Swerling模型,但考慮到相干檢測只能用于Swerling模型中的Ⅰ、Ⅲ型,且為對比分析非相干檢測性能,本文只分析Swerling Ⅰ、Ⅲ型,即自由度為2、4的卡方分布模型[24]。

2.1 Swerling I型

根據文獻[25],Swerling Ⅰ型表示由多個均勻獨立散射子目標的組合,典型目標如前向觀察的小型噴氣飛機等。它的起伏特性為慢起伏,接收到的目標回波在任意一次掃描期間都是恒定的(完全相關),但是從一次掃描到下一次掃描完全是獨立的(不相關的)。假設不計天線波束形狀對回波信號幅度的影響,截面積σ服從自由度為2的卡方分布:

(2)

由于回波信號幅度A在一次掃描內有恒定的幅度值,掃描與掃描間則為瑞利分布,由于A2=σ,所以關于回波信號幅度A的概率密度函數為:

(3)

2.2 Swerling Ⅲ型

Swerling Ⅲ型表示一個占支配地位的大隨機散射體與其他均勻獨立散射體組合的目標,典型目標如螺旋槳飛機及直升機等,它的起伏特性也是慢起伏,其RCS服從自由度為4的卡方分布:

(4)

對應回波信號幅度A的概率密度函數為:

(5)

3 非相干檢測原理與性能

3.1 非相干檢測原理

假定噪聲是均值為0,功率為ψ2的加性高斯白噪聲,且和幅度為A的正弦信號是不相干的。正弦信號功率為A2/2,單個脈沖的信噪比可定義為SNR=A2/2ψ2。傳統雷達采用非相干檢測,根據雷達經典理論,其回波信號幅度的概率密度函數為[25]:

(6)

r為回波信號加噪聲的包絡,I0(β)是宗量為β的零階修正貝塞爾函數:

(7)

式(6)是一個萊斯分布函數,如果SNR=0時,即只有噪聲輸入的時候,變為瑞利分布:

(8)

當SNR非常大,即接收機輸入為回波信號和加性高斯白噪聲時,式(6)可近似成均值為A,方差為ψ2的高斯分布:

(9)

圖1 回波信號和噪聲概率密度分布圖Fig.1 Probability density distribution of echo signal and noise

以縱軸表示概率,橫軸表示相對電平,采用非相干檢測時,若SNR為10 dB,給定噪聲功率為1,回波信號和噪聲的概率密度分布如圖1(a)所示。

3.2 非起伏目標檢測性能

目標RCS的大小與雷達的檢測性能有直接關系,在工程計算中常把RCS視為常數,即非起伏目標。當噪聲和回波信號一起輸入到接收機時,對非起伏目標采用非相干檢測,根據文獻[17]和[25]可知,虛警概率為:

(10)

式(10)中,若已知虛警概率PF,非相干檢測門限VTN可表示為:

(11)

檢測概率為:

(12)

式(12)計算非常復雜,但如果給定的虛警概率較小,而要達到的檢測概率很大時,可用多種近似方法代替,目前最精確的近似值為[26]:

(13)

erfc(z)為補誤差函數:

(14)

在雷達檢測中,一般采用奈曼-皮爾遜(Neyman-Pearson)準則:先設定一個虛警概率,然后使得檢測概率最大。給定虛警概率PF時,式(13)是關于單次檢測概率PD的表達式,它是關于SNR的函數。

(15)

式(15)中,

(16)

上式為不完全伽馬函數的皮爾遜形式,VTNI表示非相參積累的檢測門限,即:

VTNI=Γ-1(1-PF,N)

(17)

式(17)中I-1表示不完全伽馬函數的反函數,檢測概率為:

(18)

3.3 起伏目標檢測性能3.3.1 Swerling Ⅰ型

Swerling Ⅰ型起伏目標的非相干檢測概率公式由Swerling提出:

PD=e-VTN/(1+SNR),N=1

(19)

(20)

3.3.2 Swerling Ⅲ型

Swerling Ⅲ型起伏目標的非相干檢測概率公式由Marcum提出,當脈沖積累數N=1或2時,檢測概率為:

(21)

(22)

式(21)中,N=1時VT=VTN,N=2時VT=VTNI;當N>2時,檢測概率為:

(23)

上式和式(20)都是對信號進行非相參積累后再進行非相干檢測,完全未利用回波信號的相位信息。

4 相干檢測原理及其非起伏目標檢測性能

4.1 相干檢測原理

在隱身飛機等起伏目標的挑戰下,傳統非相干檢測雷達的性能已經遠遠達不到作戰需求。通過增大發射功率等途徑提高檢測性能往往也受制于戰場環境,這就需要從雷達檢測原理出發,研究如何提高雷達的檢測性能。

雷達的檢測概率取決于檢測門限和檢測信噪比,要提高雷達的檢測概率,就需要降低檢測門限或者增加信噪比。雷達采用非相干檢測時,噪聲服從瑞利分布,回波信號服從高斯分布。降低檢測門限可提高檢測概率,但也會造成虛警概率的增加,往往得不償失。根據文獻[25],對于增加信噪比,如果依據雷達距離方程考慮雷達的“硬”對抗,方法僅包含兩類:增大有效發射功率、減小接收系統噪聲。由于電子對抗措施的存在,增大有效發射功率太易暴露且易遭受攻擊,而減小接收系統噪聲往往依賴于材料技術的進步,隨著高電子遷移率場效應管(High Electron Mobility Transistor, HEMT)的大量采用,接收系統噪聲水平難以再大幅降低。

在以上兩種增加檢測概率的方法均受限制的情況下,要提高雷達的檢測性能,就需要從檢測的基本原理出發,重新設計目標存在和目標不存在時的兩個概率密度函數。根據雷達原理可知,雷達的檢測方式有兩種,即非相干檢測和相干檢測。非相干檢測只利用了回波信號的振幅信息,而相干檢測不僅利用了振幅信息,還利用了其相位信息。對信號的先驗信息知道得越多,檢測性能越好,所以相干檢測的性能理應優于非相干檢測。

非相干檢測時,變量r表示回波信號加噪聲的包絡;假定回波信號完全相干,對于相干檢測,已知信號的相位信息,此時r表示回波信號和噪聲的綜合幅度。若雷達采用相干檢測,當有回波信號時,其幅度的概率密度函數是均值為A,方差為ψ2的高斯分布:

(24)

當只有噪聲時,是均值為0,方差為ψ2為高斯分布:

(25)

采用相干檢測時,仍舊給定SNR為10 dB,噪聲功率為1,回波信號和噪聲的概率密度曲線如圖1(b)所示。對比圖1(a),噪聲與回波信號的概率密度函數分得更開,其交點處概率值更低,也更容易選擇合適的檢測門限,達到更大的檢測概率。

4.2 非起伏目標相干檢測性能

對于非起伏目標,相干檢測的檢測概率要高于非相干檢測,文獻[29]已經給出了證明。虛警概率為:

(26)

式(26)中VTC代表相干檢測門限,若已知虛警概率時,可表示為:

(27)

式(27)中erfc-1(z)為補誤差函數的反函數,檢測概率為:

(28)

檢波器以前進行的積累稱為相參積累,檢波器不僅利用了回波信號的幅度信息,還利用了相位信息。在理想情況下,對N個回波脈沖進行相參積累,積累后的信噪比可以改善N倍,積累后的信噪比變為(SNR)CI=N×SNR,(SNR)CI為N個脈沖相參積累后的信噪比,將式(28)的SNR替換為(SNR)CI即為相參積累時的檢測概率。

需要指出的是,非相干檢測時,相參積累只是在短時內忽略了回波信號的相位變化,注定了積累時間不會太長。而相干檢測時,保證了回波信號的相位始終相關,可積累的時間更長。

對回波信號先進行相參積累,對積累后的信號進行相干檢測,完全利用了其相位信息,因而能達到的檢測概率也更大。

5 起伏目標相干檢測性能

對于起伏目標的相干檢測,回波信號可視為目標RCS變量σ和綜合幅度r的聯合概率密度函數。用f(σ)表示目標RCS的概率密度函數,f(r)表示綜合幅度的概率密度函數,所以聯合概率密度函數為f(r,σ)。根據簡單的概率論知識,得到f(r)需要用到以下關系式:

f(r,σ)=f(r/σ)f(σ)

(29)

(30)

將式(29)代入式(30)可得

(31)

要求得起伏目標相干檢測的檢測概率,實質上是求二重積分的解。先對f(σ)從0至無窮大積分,得到綜合幅度r的一維概率密度函數f(r),再對f(r)從檢測門限VTC到無窮大積分即可得到起伏目標的檢測概率。同式(13)一樣,在給定虛警概率時,所得到是單次檢測概率PD的值,是關于SNR的函數。

又因為目標RCS與回波功率成比例,若以回波信號幅度A為變量,檢測概率的求解方式相同,式(31)變為:

(32)

N個脈沖積累后進行相干檢測,單個脈沖信噪比為(SNR)CI,將A固定不變,根據信噪比定義可得脈沖積累后的ψ,帶入式(27)可得VTC,再根據式(31)和(32)可得脈沖積累時的相干檢測概率。

計算檢測概率前,首先要知道檢測門限值。非相干檢測的門限值計算,單個脈沖檢測時為式(11),脈沖積累時為式(17),需要計算不完全伽馬函數的皮爾遜形式的反函數等復雜公式,大大增加了計算的復雜度。而對于相干檢測,其檢測門限的計算公式為式(27),對比于非相干檢測,公式數量少,計算復雜度低。在檢測方式相同時,檢測門限值只與虛警概率和脈沖積累數有關,而與目標是否起伏無關。

基于重積分方程的相干檢測,實際計算時只需求得二重積分的解,且計算方法具有通用性。以RCS變量σ為例,根據不同的起伏模型,將對應的f(σ)帶入式(31),可得到不同起伏目標的檢測概率。相對于起伏目標的非相干檢測復雜的計算公式,公式簡潔且計算復雜度低。

綜上,對于起伏目標,相干檢測的檢測門限和檢測概率的計算難度均小于非相干檢測。下面以Swerling型起伏目標為例,對起伏目標相干檢測性能進行分析。

5.1 Swerling Ⅰ型

根據起伏目標相干檢測概率計算方法,對Swerling Ⅰ型起伏目標信號進行相干檢測,回波信號的概率密度函數按RCS形式可表示為:

(33)

將式(2)、(33)帶入式(31),可得:

(34)

根據式(32),按幅度形式可以表示為:

(35)

所以,對Swerling Ⅰ型起伏目標進行相干檢測,檢測概率為:

(36)

5.2 Swerling Ⅲ型

對Swerling Ⅲ型起伏目標進行相干檢測,回波信號的概率密度函數按RCS形式仍為式(33),將式(4)、(33)帶入式(31)可得:

(37)

根據式(32),按幅度形式,回波信號的概率密度函數可以表示為:

(38)

所以,對Swerling Ⅲ型起伏目標進行相干檢測,檢測概率為:

(39)

6 仿真驗證與分析

對非起伏目標進行非相干檢測時,式(13)為單個脈沖非相干檢測理論值,式(18)為脈沖積累時的理論值;相干檢測時,式(28)為理論值。對于Swerling Ⅰ型起伏目標,非相干檢測時,式(19)為單個脈沖非相干檢測的理論值,式(20)為脈沖積累時的理論值;相干檢測時,式(36)為理論值。對于Swerling Ⅲ型起伏目標,非相干檢測時,脈沖積累數N=1、2時,式(21)為理論值,N>2時,式(23)為理論值;相干檢測時,式(39)為理論值。

6.1 相干檢測仿真驗證

本次仿真主要通過MATLAB對推導的相干檢測公式進行仿真驗證及分析,仿真時假設目標回波完全相干。對于三種目標的非相干檢測和非起伏目標的相干檢測,MATLAB中已經有內置函數rocpfa,且其理論已經相當成熟,因此不再對非相干檢測進行仿真驗證。對于起伏目標相干檢測的理論值,主要用到的函數為erfcinv和integral2,前者對應erfc-1(z),可計算出檢測門限,后者為二重積分函數。已知檢測門限時,根據式(36)、(39)可計算出檢測概率。

下面對非起伏目標和Swerling Ⅰ、Ⅲ型起伏目標的相干檢測,通過MATLAB進行仿真驗證。對于非起伏目標,理想情況下,經相干檢測后,在接收機輸出處,回波信號的綜合幅度服從方差為ψ2,均值為A的高斯分布。對于Swerling Ⅰ、Ⅲ型起伏目標,回波信號的綜合幅度在掃描與掃描間分別服從自由度為2和4的卡方分布,相干檢測后仍舊服從高斯分布。

對于非起伏目標的相干檢測仿真,隨機產生一定長度的0、1信息序列,0代表無信號,1代表有回波信號且幅度為1,輸入至高斯信道模型,信道輸出序列是被干擾了的碼字序列。根據奈曼-皮爾遜準則,以SNR為參數,給定虛警概率,根據虛警概率求出檢測門限,分別對被干擾了的碼字序列進行門限判決,高于門限的信號判定為回波信號,統計出正確檢測數,最后除以初始回波信號數,多次仿真后求平均數,得到不同SNR對應的檢測概率。

對起伏目標相干檢測仿真時,生成的“1”碼元需要先輸入至卡方分布模型,得到自由度為2或4,平均幅度為1的碼字數列后,再輸入至高斯信道模型進行門限判決和統計檢測概率。

根據文獻[27],雷達通常工作時,虛警概率PF的值一般非常小,根據系統的種類不同,PF的典型值不會超過10-3,通常介于10-6和10-8之間,仿真時仍根據這個結論設定虛警概率值。

仿真分為兩部分:單個脈沖檢測時,分別給定虛警概率PF為10-2、10-4和10-8,以SNR為參數,研究虛警概率對檢測性能的影響;脈沖積累時,虛警概率固定為10-8,給定脈沖積累數N分別為10和50時,以SNR為參數,研究脈沖積累數N對檢測性能的影響。

對于非起伏目標和Swerling Ⅰ、Ⅲ型起伏目標的相干檢測,如圖2、圖3和圖4所示。單個脈沖檢測時,給定不同的虛警概率,以SNR為參數,在SNR相同時,相干檢測的檢測概率均大于非相干檢測。也可以這樣說,要達到相同檢測概率,相干檢測所需的SNR均小于非相干檢測。同時,對于兩種檢測方式,給定虛警概率越小,達到相同檢測概率所需的SNR越高。

圖2 單個脈沖信噪比對應的檢測概率(非起伏)Fig.2 Probability of detection versus single pulse SNR. Non-fluctuating targets

圖3 單個脈沖信噪比對應的檢測概率(Swerling Ⅰ型)Fig.3 Probability of detection versus single pulse SNR. Swerling Ⅰ

脈沖積累時,虛警概率固定為10-8,給定不同的脈沖積累數,以SNR為參數,SNR相同時,相干檢測的檢測概率均大于非相干檢測。換句話說,要達到相同檢測概率,相干檢測所需的SNR均小于非相干檢測。起伏目標的脈沖積累改善和非起伏目標相同。

圖4 單個脈沖信噪比對應的檢測概率(Swerling Ⅲ型)Fig.4 Probability of detection versus single pulse SNR. Swerling Ⅲ

圖5 單個脈沖信噪比對應的檢測概率對比Fig.5 Probability of detection versus single pulse SNR.Non-fluctuating versus Fluctuating

綜上,對于非起伏目標和起伏目標,分別給定相同的虛警概率和脈沖積累數,相干檢測的性能均優于非相干檢測,仿真結果驗證了理論推導。

6.2 起伏目標與非起伏目標對比

對于非起伏目標和Swerling Ⅰ、Ⅲ型起伏目標,給定虛警概率為10-8,分別給定脈沖數N為1、10,以SNR為參數,要達到的檢測概率為0.9時,所需的單個脈沖信噪比SNR如表1所示。

表1 檢測概率為0.9時所需的單個脈沖信噪比

N為1和10時,非相干檢測所需的SNR均大于相干檢測。在N由1增加為10,分別對Swerling Ⅰ、Ⅲ型起伏目標進行相干檢測,要達到相同檢測概率時,所需的SNR減小了10 dB,而非相干檢測所需的SNR的減少量小于10 dB。三種目標模型,脈沖積累數N為10時,相對于非相干檢測,相干檢測所需的SNR均減小了大約2 dB。

6.3 兩種檢測方式的起伏損失

目標的起伏引起目標RCS變小,進而造成SNR變小,檢測概率隨之大幅降低。起伏損失定義為給定相同虛警概率和檢測概率,起伏目標相比于非起伏目標SNR額外增加的值。

本文參考了文獻[17]的提出的計算方法,該方法的誤差可以控制在0.005 dB以內。對非起伏目標和Swerling Ⅰ、Ⅲ型起伏目標分別計算出對應的檢測概率PD的SNR,在檢測概率相同時,用起伏目標對應的SNR減去非起伏目標的SNR即為對應檢測概率PD的起伏損失。本文利用這個方法分析了相干檢測和非相干檢測的起伏損失。

給定虛警概率PF為10-8,單個脈沖檢測時,所得結果如圖6所示:在檢測概率較低時,非相干檢測的起伏損失值會略小于相干檢測;在檢測概率較大時,起伏損失趨于一致。

圖6 起伏損失與檢測概率的關系Fig.6 Fluctuation loss versus detection probability

對起伏損失的分析說明,相干檢測可以達到比非相干檢測更高的檢測概率,而兩者承受的起伏損失卻近似相等,所以其性能更好。

7 結論

本文給出了基于積分方程計算起伏目標相干檢測概率的方法,并以Swerling Ⅰ、Ⅲ型起伏目標為例,推導出了其相干檢測概率計算公式。通過仿真驗證,得到以下結論:對于同一起伏目標模型,給定相同虛警概率和脈沖積累數:a)以信噪比為參數,采用相干檢測的檢測概率要高于非相干檢測,且脈沖積累數越大,檢測概率更高;b)以檢測概率為參數,在檢測概率較大的情況下,相干檢測和非相干檢測的起伏損失是近似相等的。綜上,對于起伏目標,相干檢測的檢測性能要優于非相干檢測。雖然Swerling Ⅰ、Ⅲ模型已經不太適合于現代飛機,但其作為自由度分別為2、4的卡方分布,改變自由度即可擬合其他起伏目標模型,因而本文對于提高隱身飛機等目標的檢測概率有一定參考意義。實際應用過程中,應該根據目標起伏特性選擇最合適的起伏模型和脈沖積累數,以期達到最佳的檢測性能。

猜你喜歡
虛警信噪比脈沖
頻率步進連續波雷達電磁輻射2階互調虛警干擾效應規律
脈沖離散Ginzburg-Landau方程組的統計解及其極限行為
兩種64排GE CT冠脈成像信噪比與劑量對比分析研究
上下解反向的脈沖微分包含解的存在性
一種電阻式應變傳感器的數據處理方法
基于深度學習的無人機數據鏈信噪比估計算法
黃芩苷脈沖片的制備
低信噪比下基于Hough變換的前視陣列SAR稀疏三維成像
空管自動化系統二次代碼劫機虛警分析
民用飛機貨艙煙霧虛警淺析
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合