?

EMVS互質面陣張量波束成形

2024-03-07 13:04周成偉史治國
信號處理 2024年2期
關鍵詞:波達面陣互質

鄭 航 周成偉*,2 王 勇 史治國,3

(1.浙江大學信息與電子工程學院,浙江杭州 310027;2.浙江省協同感知與自主無人系統重點實驗室,浙江杭州 310015;3.浙江大學國際聯合創新中心,浙江海寧 314400)

1 引言

在雷達、遙感、無線通信、射電天文等應用場景中,通過部署傳感器陣列對空間信源的波達方向進行精準測算,是實現目標測向、定位、跟蹤等諸多功能的基本前提[1-3]。波束成形作為一項基礎性陣列信號處理技術,能夠對傳感器陣列接收信號進行空域濾波形成指向性空間波束,通過波束掃描得到不同方向的空間響應,從而實現信源測向[4-7]。波束成形技術能夠有效抵抗干擾且易于部署,已被廣泛應用于各類傳感器陣列系統中,長期以來發揮著重要作用。隨著實際應用場景中對信源測向性能需求的不斷提高,現有系統所布設的傳感器陣列逐漸演變出新的特點。一方面,為了在系統軟硬件資源受限的條件下提升信號感知能力,稀疏化陣元排布的陣列架構開始得到關注[8-11],其中,互質陣列作為一種擁有系統化結構的典型稀疏陣列,能夠突破傳統均勻陣列信號處理在奈奎斯特采樣速率上的限制,在使用相同數量傳感器陣元的條件下擁有更大的陣列孔徑,從而有效提升波束分辨率[12-17];另一方面,為了滿足復雜信源測向場景感知空間極化信息的需求,電磁矢量傳感器(electromagnetic vector sensor,EMVS)陣列相比于傳統的標量傳感器陣列,能夠獲取入射信號的多極化狀態,通過利用更加全面的空間電磁信息與更強的干擾抵抗性,提供更為精準的波束控制與角度估計能力[18-19]。為此,如何發掘互質陣列和EMVS 陣元的獨特優勢,形成相匹配的波束成形算法設計框架,是當前一個國際前沿研究熱點。

為了發揮互質陣列的大孔徑優勢并降低波束成形器的軟硬件復雜度,面向互質線陣的自適應波束成形算法受到廣泛關注[20-23]。其中,文獻[23]結合互質線陣參數估計與稀疏協方差矩陣重建思想,提出了一種面向互質線陣的魯棒自適應波束成形算法。為了利用互質陣列波束成形實現信源測向,文獻[24]針對互質線陣的陣元稀疏排布特點,從互質線陣稀疏子陣波束譜的分布特性出發,提出了基于乘性準則和最小化功率準則的波束功率合成方法,通過虛峰抑制的波束掃描測算信源波達方向?;谠摷夹g框架,相關研究針對交叉項消除[25]、分辨率增強[26]和水聲探測[27]等具體問題及實際應用環境特點提出了相適應的改進方案,促進了互質線陣波束功率合成方法在信源測向上的推廣應用。另一方面,為了利用EMVS 陣列的極化域信息,文獻[28]提出了一種基于廣義旁瓣相消的EMVS 陣列波束成形方法,并證明了該方法相比于傳統Capon波束成形算法具有更小的波達方向估計誤差。文獻[29]考慮極化方向的魯棒約束,設計能夠應對陣元位置失配的EMVS 陣列魯棒波束成形器,進一步提升了波達方向的估計精度??紤]到互質陣列和EMVS 陣元均能夠有效促進波束成形輸出性能的提升,融合互質陣元排布特點與EMVS 陣元的EMVS 互質線陣應用研究開始興起。相關研究工作結合增廣虛擬域處理[30-31]與多重信號分類(multiple signal classification,MUSIC)方法,驗證了部署EMVS 互質線陣進行信源測向的可行性[32];文獻[33]基于EMVS 互質稀疏子陣信號的平行因子分析實現了信源測向的性能提升。為了滿足三維空間的信源測向需求,EMVS 互質線陣亟待向EMVS互質面陣演變,以便感知更大范圍、更高維度的空間信息。針對互質面陣信源測向需求,基于虛擬域稀疏恢復的方法能夠在克服相位模糊的前提下實現二維波達方向估計[34];文獻[35]利用一種互質雙平行陣列,在虛擬域上應用root-MUSIC方法以完成信源波達方向搜索。然而,這些方法將多快拍接收信號堆疊成矩陣并作后續統計處理。在EMVS互質面陣場景中,由于接收信號涵蓋了多維度信息,基于矩陣信號建模與處理的方法將破壞多維接收信號的原始結構,進而影響信源測向的性能。為此,如何設計匹配EMVS互質面陣多維信號結構的空域濾波算法,已成為當前一個亟待解決的研究難點問題。

張量作為一種多維度的數據結構,被用于表征三維及三維以上的數據。豐富的張量分解工具為張量數據的高維特征提取和分析提供了技術支撐[36]。近年來,張量模型被廣泛應用于陣列信號處理領域,用于表征多維信號并實現張量化參數估計[37-40]。通過引入張量進行EMVS互質面陣接收信號建模,能夠有效保留其所涵蓋的復雜空間/時間/電磁特征,并利用張量代數工具,在不破壞信號原始結構的前提下實現高性能的角度信息提取??紤]到張量模型的優越性,面向多維信號的張量波束成形算法成為一個新興的研究熱點。具體而言,文獻[41-43]分別提出了面向均勻面陣的最小均方誤差(minimum mean square error,MMSE),線性約束最小方差(linear constrained minimum variance,LCMV)以及最小方差無畸變響應(minimum variance distortionless response,MVDR)等傳統波束成形設計準則的張量表征方法,并利用張量分解對所設計張量權重優化問題進行求解,驗證了對張量信號進行結構化空域濾波的可行性。進一步地,上述張量波束成形算法亦被推廣應用于極化敏感陣列及大規模分布式陣列等復雜場景中[44-45]。然而,現有方法僅適用于滿足奈奎斯特采樣速率的均勻面陣場景,而在EMVS 互質面陣場景下,陣元的稀疏排布將導致模型失配問題,由此引入的虛峰將嚴重干擾張量波束成形器對信號源波達方向的精準刻畫。文獻[46-47]提出了面向互質面陣的稀疏張量模型,分別從虛擬域張量信源辨識度優化和非連續虛擬域張量填充的角度,實現了高自由度、高精度的二維波達方向估計。然而,這些方法沒有考慮部署EMVS互質面陣。因此,如何匹配EMVS 稀疏陣元排布特點進行虛峰抑制,是利用EMVS 互質面陣張量波束成形進行有效信源測向的關鍵。

本文提出了一種面向EMVS 互質面陣的張量波束成形算法,實現了基于EMVS 互質面陣接收信號的多維空域濾波與信源精準測向。首先,將構成EMVS 互質面陣的稀疏均勻子面陣接收信號表示為一對涵蓋多維度空間電磁信息的張量。在張量信號建模的基礎上,構建面向稀疏均勻子面陣張量信號的結構化空域濾波準則,并設計相應的張量化最小方差無畸變響應優化問題。為了有效求解該優化問題,對張量權重進行canonical polyadic分解,將原始張量化最小方差無畸變響應優化問題轉換為對應波達方向信息維度與極化狀態信息維度的子問題,并根據這些子問題所對應空間維度信息的耦合關系設計交替迭代求解方法,以獲得張量權重的全局最優解。進而,基于稀疏均勻子面陣的陣元互質排布特點分析,理論證明各子面陣所對應虛峰的互不重疊特性。從抑制虛峰的角度出發,設計波束功率的互質合成方法,以此構造EMVS 互質面陣張量波束功率圖,在波達方向平面上精準刻畫信源空間方位。所提算法既能夠充分利用EMVS 互質面陣接收信號的多維結構化信息,又能夠匹配其陣元稀疏排布特點實現虛峰抵消。仿真結果表明,所提算法具有相較于EMVS 均勻面陣張量波束成形算法更好的信源測向效果。

2 EMVS互質面陣的張量信號建模

如圖1 所示,EMVS 互質面陣P 由兩個稀疏均勻子面陣P1和P2組成,P1和P2分別包含M1×N1和M2×N2個傳感器陣元,且(M1,M2)和(N1,N2)為兩對互質整數。稀疏均勻子面陣P1的傳感器陣元在x軸和y軸方向上的間隔分別為M2d和N2d,其中d=λ/2 為單位陣元間距,λ表示信號波長。類似地,稀疏均勻子面陣P2的傳感器陣元在x軸和y軸方向上的間隔分別為M1d和N1d。由此,P1中第(m1,n1)個傳感器陣元在x軸和y軸方向上的位置分別為x1(m1)=(m1-1)M2d和y1(n1)=(n1-1)N2d,其 中m1=1,2,…,M1,n1=1,2,…,N1;P2中第(m2,n2)個傳感器陣元在x軸和y軸方向上的位置分別為x2(m2)=(m2-1)M1d和y2(n2)=(n2-1)N1d,其中m2=1,2,…,M2,n2=1,2,…,N2。由于 兩個稀疏 均勻子面陣的傳感器陣元排布滿足互質條件,P1和P2只有在位于原點處的傳感器陣元(m1=n1=m2=n2=1)相互重合,而其余傳感器陣元在xoy平面內均不重合;因此,EMVS互質面陣P中傳感器陣元的數量為M1N1+M2N2-1。同時,與傳統標量傳感器不同的是,EMVS 互質面陣中各傳感器陣元利用三個相互正交的電偶極子和三個相互正交的磁偶極子來實現極化狀態的感知,形成六路輸出以描述完備的電磁場信息。

圖1 EMVS互質面陣結構示意圖(以M1=N1=5,M2=N2=4為例)Fig.1 EMVS coprime planar array geometry(M1=N1=5,M2=N2=4)

假設K個遠場窄帶非相關信號源從{(θk,?k),k=1,2,…,K}方向入射至EMVS 互質面陣,其中θk和?k分別表示第k個信號源的方位角和俯仰角,且θk?[-π/2,π/2],?k?[-π/2,π/2]。EMVS 互質面陣中各傳感器陣元的六路輸出同時包含了波達方向信息和極化狀態信息,可表示為一個電磁響應矢量:

分別表示對應第k個信號源的空域響應矩陣和極化狀態矢量,γk?[0,2π]和ηk?[-π,π]分別表示對應第k個信號源的極化輔助角和極化相位差,且

為了保留稀疏均勻子面陣Pi(i=1,2)在t時刻接收信號的多維電磁信息,即x軸方向、y軸方向的波達方向信息以及極化狀態信息,將其建模為一個三維張量:

分別表示EMVS互質面陣在x軸和y軸方向上的導引矢量,對應第k個信號源,且μk=sin?kcosθk,νk=sin?ksinθk,sk(t)為t時刻的采樣信號,t=1,2,…,T,T為采樣快拍數,?表示外積操作,(?)Τ表示轉置操作,Ni(t)為t時刻獨立同分布的加性高斯白噪聲張量。進一步地,將所采集T個采樣快拍的三維張量信號{Xi(t),t=1,2,…,T}在第四維度(即時間維度)上進行疊加,構造出稀疏均勻子面陣Pi的四維張量接收信號:

其中,sk=[sk(1),sk(2),…,sk(T)]Τ?CT表示第k個信號源的多快拍采樣信號矢量,Ni表示四維高斯白噪聲張量。

3 面向EMVS 互質面陣的張量波束成形器設計

本節介紹所提面向EMVS互質面陣的張量波束成形器,包括稀疏均勻子面陣張量信號的空域濾波準則,張量化最小方差無畸變響應優化問題分解策略及對應的交替迭代求解方法,基于張量空間質數分解唯一性定理的稀疏均勻子面陣波束分布特性分析,以及基于波束功率互質合成的張量波束成形輸出。

3.1 張量信號空域濾波準則設計

由于構成EMVS 互質面陣的兩個稀疏均勻子面陣中的陣元排布不滿足奈奎斯特采樣定理,由此帶來的相位模糊將導致波束功率圖上產生虛峰,進而直接影響信源測向性能。為此,本小節從面向稀疏均勻子面陣張量信號的空域濾波準則入手,提出張量化最小方差無畸變響應優化問題,為后續剖析虛峰在波達方向平面上的分布特性并設計張量波束功率的互質合成方法提供理論基礎。

具體地,針對稀疏均勻子面陣Pi在t時刻的張量接收信號Xi(t),設計一個具有相同三維結構的張量權重Wi?CMi×Ni×6對其進行加權,如圖2 所示,輸出信號yi(t)可以表示為如下形式:

圖2 稀疏均勻子面陣張量信號的張量權重加權示意圖Fig.2 Illustration of spatial filtering of the sparse subarray tensor signals with a tensor weight

圖3 張量波束功率互質合成處理框架示意圖Fig.3 Illustration of coprime synthesis of tensor beam power

為了獲得兩個稀疏均勻子面陣所對應的張量權重Wi,需要最小化張量波束成形的平均輸出功率,并保證波束掃描方向信號無失真,具體優化問題可表示為:

表示輸出信號功率,

表示張量信號Xi(t)的六維協方差張量,

表示稀疏均勻子面陣Pi對應第k個信號源的三維導引張量,| · |表示復數的求模操作,E[·]表示取期望操作。這里,分別表示波束掃描的方位角和俯仰角,。在實際應用中,六維協方差張量Ri可通過計算Xi(t)的采樣協方差張量近似得到,即:

式(9)所構造張量化最小方差無畸變響應優化問題以張量信號加權輸出功率為優化目標,在極化輔助角γk和極化相位差ηk理想已知的條件下,當波束掃描至信號源波達方向上,即時,所對應張量波束功率呈現為最大響應,即在張量波束功率圖對應位置上形成主瓣;由此,通過遍歷波束掃描方位角和俯仰角的取值范圍,能夠在波達方向平面上有效反映信源的方位。

3.2 基于張量權重分解的交替迭代求解方法

由于優化問題(9)的優化目標和約束條件中包含張量外積和內積項,因此該問題是一個非凸優化問題。傳統優化方法以矩陣信號處理和矩陣代數為基本前提,無法對所提張量化最小方差無畸變響應優化問題進行直接求解。同時,優化問題(9)的非凸屬性也導致其無法提供一個閉式的全局最優解,嚴重影響了所提算法的實用性。為此,本文提出對張量權重進行分解,將原始優化問題(9)轉換為對應波達方向信息維度和極化狀態信息維度的子問題,并設計相應的交替迭代求解方法。具體而言,將三維張量權重Wi與Xi(t)的各空間維度信息逐一對應,可將Wi用canonical polyadic 分解的方式表示為對應于x軸波達方向信息、y軸波達方向信息和極化狀態信息的波束成形權重矢量,具體形式如下:

將式(14)代入式(8)中,則采用張量權重Wi對張量信號Xi(t)的加權可等價表示為三個權重矢量對其進行結構化加權,如圖2所示,則輸出信號yi(t)可等價表示為:

其中,×r表示沿著第r維度的張量-矩陣內積操作。

基于式(15)的張量信號結構化加權形式,利用對應稀疏均勻子面陣Pi任意兩個維度的權重矢量對Xi(t)進行加權,則相應的剩余維度,即x軸波達方向信息維度、y軸波達方向信息維度和極化狀態信息維度的輸出信號可分別表示為:

相應地,稀疏均勻子面陣Pi的張量波束成形輸出信號yi(t)可等價表示為:

其中()·H表示共軛轉置操作。將式(19)中輸出信號yi(t)的三種等價表示形式代入張量化最小方差無畸變響應優化問題(9)中,可將其分解為以下三個子問題:

分別表示對應x軸波達方向信息維度、y軸波達方向信息維度和極化狀態信息維度的協方差矩陣。由于這些子問題的優化目標和約束條件中均不包含張量,因此它們是凸優化問題,能夠通過拉格朗日乘子法有效求得它們的局部最優解。

進而,根據這些子問題所對應空間維度信息的耦合關系,可以對它們進行交替迭代求解得到張量權重的全局最優解。具體而言,子問題(20)優化目標中的Rxi包含了其他兩個子問題(21)、(22)中的優化變量wyi,wpi,子問題(21)優化目標中的Ryi包含了其他兩個子問題(20)、(23)中的優化變量wxi,wpi,而子問題(22)優化目標中的Rpi包含了其他兩個子問題(20)、(21)中的優化變量wxi,wyi,因此基于交替迭代優化的思想對上述三個子問題進行聯合求解。首先,將波束成形權重矢量wxi,wyi和wpi初始化為,根據式(23)計算得到對應x軸波達方向信息維度的協方差矩陣Rxi;隨后,利用拉格朗日乘子法求得子問題(20)的局部最優解:

將式(26)中的局部最優解wˉxi代入式(24)中計算協方差矩陣Ryi,并通過拉格朗日乘子法求得子問題(21)的局部最優解:

重復式(26)~(28)的優化過程交替迭代更新三個權重矢量,并在每次迭代中根據式(15)計算輸出信號yi(t)的功率σ=|yi(t)|2。在第η次迭代中(η≥2),對比輸出信號功率ση與前一次迭代所對應輸出信號功率ση-1的差值,當Ξη=ση-ση-1≤δ時,上述三個子問題的交替迭代優化收斂,其中δ>0 表示一個預定義的收斂閾值。由此,通過序貫求解子問題(20)~(22),可以得到面向稀疏均勻子面陣Pi的全局最優張量權重以此為基礎,下一小節通過研究兩個稀疏均勻子面陣的波束分布特性,設計匹配陣元稀疏排布特點的EMVS 互質面陣張量波束功率處理方法,從而實現虛峰抵消的張量波束成形。

3.3 稀疏均勻子面陣的波束分布特性分析

定理1P1和P2所對應張量波束功率在信號源方向上形成主瓣,而在其他位置形成的虛峰互不重疊。

根據式(5)、(6)中導引矢量ai(μk),ai(νk)的定義,結合式(30)、(31)的等價關系,從ai(μk),ai(νk)的指數項推導可得[48]:

其中,P1,P2,Q1,Q2為非零整數。從式(32)、(33)中可得出以下關系:

然而,由于互質整數對(M1,M2)和(N1,N2)的最大公約數均為1,并不存在能夠滿足式(34)等價條件的非零整數P1,P2,Q1,Q2。因此,除了信號源方向(θk,?k)以外,并不存在其他模糊方向能夠使得張量波束功率在對應位置上同時輸出最大響應。由此可知,僅在對應信號源波達方向(θk,?k)的位置上形成主瓣,而在其他位置形成的虛峰互不重疊。證畢。

定理1 基于張量權重的可分解性,將張量化無畸變響應約束投影至不同維度上,利用張量空間的質數分解唯一性證明了EMVS 互質面陣中兩個稀疏均勻子面陣所對應虛峰的互不重疊特性。該結論為后續通過稀疏均勻子面陣波束功率的互質合成實現虛峰抵消的張量波束成形提供了基礎。

3.4 基于波束功率互質合成的張量波束成形

4 仿真實驗

首先,通過繪制張量波束功率圖開展仿真實驗對比,以驗證所提EMVS 互質面陣張量波束成形算法的有效性??紤]如圖1 所示共包含40 個陣元的EMVS 互質面陣結構,假定兩個信號源的入射方向分別為(θ1,?1)=(15.5°,20.5°)和(θ2,?2)=(45.5°,60.5°),對應已知的極化輔助角和極化相位差為(γk,ηk)=(40°,20°),k=1,2。所提算法的收斂閾值設置為δ=10-3,波束掃描的角度間隔設置為Δ=0.1°。在信源信噪比SNR=5 dB,采樣快拍數T=300的條件下,在圖4(a)、(b)中繪制稀疏均勻子面陣P1和P2的張量波束功率圖,在 圖4(c)中繪制基于乘積最小化處理的互質合成張量波束功率圖,并將其分別投影至方位角維度和俯仰角維度上,如圖4(d)、(e)所示。

圖4 所提EMVS互質面陣張量波束成形器的張量波束功率圖Fig.4 Τensor beam power patterns for the proposed EMVS coprime planar array tensor beamformer

從圖4(a)、(b)可以看出,得益于EMVS 六路輸出中所涵蓋的豐富空間電磁信息,以及所提張量信號的結構化空域濾波準則,兩個稀疏均勻子面陣的張量波束功率均在對應信號源波達方向的位置上呈現為主瓣;然而,由于陣元的稀疏排布特性,在非信號源波達方向的位置上同時呈現出多個虛峰。根據3.3 小節所分析的稀疏均勻子面陣波束分布特性,所對應的虛峰位置互不重疊。通過對兩個稀疏均勻子面陣張量波束功率進行互質合成處理,如圖4(c)、(d)、(e)所示,EMVS 互質面陣的張量波束功率在對應信號源波達方向的位置上呈現出明顯的主瓣,而在其他位置上的虛峰被有效抵消。

進一步地,對比所提EMVS 互質面陣張量波束成形器與EMVS 均勻面陣張量波束成形器[44]的波束成形輸出性能。EMVS 均勻面陣按照5 行8 列的結構排布40個陣元,以保證陣元個數與EMVS互質面陣相同。假定兩個信號源的入射方向分別為(θ1,?1)=(20.5°,25.5°) 和(θ2,?2)=(30.5°,45.5°)。在信噪比SNR=5 dB,采樣快拍數T=300的條件下,繪制兩種方法所對應的張量波束功率圖,如圖5(a)、(b)所示;在信噪比SNR=-10 dB,采樣快拍數T=300 的條件下,繪制兩種方法所對應的張量波束功率圖,如圖5(c)、(d)所示。

從圖5中可以看出,在不同信噪比場景下,所提EMVS 互質面陣張量波束成形器的輸出性能均優于采用相同陣元數目的EMVS 均勻面陣張量波束成形器。得益于EMVS 互質面陣的大孔徑優勢,所提EMVS 互質面陣張量波束成形器相較于EMVS均勻面陣能夠輸出更加精尖且能量集中的主瓣,在對應信號源波達方向的位置具備更強的波束指向性。與此同時,由于所提張量波束功率乘積最小化處理準則充分利用稀疏均勻子面陣的陣元互質排布特點削弱虛峰響應,EMVS 互質面陣張量波束成形器對應的張量波束功率圖上無明顯虛峰干擾,而EMVS 均勻面陣張量波束成形器所對應的張量波束功率圖卻仍存在較強的旁瓣,將對信源空間方位的精準刻畫造成嚴重干擾。

最后,為了驗證所提算法在信源測向精度上的優越性,將其與典型的EMVS 互質面陣虛擬域MUSIC方法[32],以及EMVS均勻面陣張量波束成形方法進行性能對比。假定兩個信號源的入射方向分別 為(θ1,?1)=(30°,25°)和(θ2,?2)=(45°,52°),在圖6中展示對這兩個信號源進行測向的50次重復實驗結果。在不同的仿真環境下,所提算法都能夠精確的估計兩個信號源的來波方向,而相比之下,虛擬域MUSIC方法存在較大的信源測向偏差。同時,在采用張量波束成形準則的條件下,EMVS均勻面陣由于孔徑受限,也無法提供準確的測向結果。綜合上述仿真結果,本文所提EMVS 互質面陣張量波束成形器充分利用了多維稀疏陣列的大孔徑優勢和陣元結構化稀疏排布特點,能夠獲得較EMVS 均勻面陣更加精尖的波束輸出,以及較傳統EMVS 互質面陣測向方法更為精確、穩定的波達方向估計結果,為實際應用場景中的高性能信源測向提供了保障。

圖6 信源測向結果對比圖Fig.6 Comparison of source direction estimations

5 結論

本文綜合考慮了EMVS 互質面陣的接收信號多維空間電磁信息結構和陣元稀疏排布特點,提出了一種EMVS 互質面陣張量波束成形算法。該算法構建了EMVS 互質面陣張量信號的結構化空域濾波準則,基于張量權重分解設計對應波達方向信息維度和極化狀態信息維度的子問題,并提出針對子問題的交替迭代求解方法獲得張量權重的全局最優解。在此基礎上,分析了由兩個稀疏均勻子面陣的陣元互質排布特點所帶來的虛峰互不重疊特性,從抑制虛峰的角度設計了稀疏均勻子面陣張量波束功率的互質合成方法,從而在對應信號源波達方向的位置形成指向性波束。仿真結果表明,所提算法獲得了優于現有EMVS 均勻面陣張量波束成形器的輸出性能。

猜你喜歡
波達面陣互質
基于互質陣列的信號波達方向估計算法
一種改進面陣紅外相機畸變校正方法研究
基于STC12C5A的雙足機器人設計
基于積分球數據的面陣航測相機影像輻射校正
一種有色噪聲背景下混合信號的波達方向估計算法
衛星導航干擾信號波達方向估計分析*
Short-range Radar Detection with(M,N)-Coprime Array Configurations
基于分離式電磁矢量傳感器陣列的相干信號波達方向估計
基于級聯MUSIC的面陣中的二維DOA估計算法
基于超聲陣列傳感器與遺傳MUSIC的局放源波達方向估計
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合