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

導航衛星星地雙基SAR時頻同步誤差估計方法

2023-05-11 12:51:02遆晶晶索志勇王婷婷趙秉吉張樂茹
西安電子科技大學學報 2023年2期
關鍵詞:信號

遆晶晶,索志勇,王婷婷,2,趙秉吉,張樂茹

(1.西安電子科技大學 雷達信號處理國家重點實驗室,陜西 西安 710071;2.四川航天電子設備研究所,四川 成都 610100;3.北京空間飛行器總體設計部,北京 100094)

1 引 言

地球同步軌道合成孔徑雷達(GEo synchronous Orbital Synthetic Aperture Radar,GEO SAR)具有寬覆蓋、高重訪、持續觀測能力強等優勢,是當前遙感領域的研究熱點,具有十分廣闊的應用前景[1-4]。與低軌合成孔徑雷達(Low Earth Orbital SAR,LEO SAR)衛星相比,GEO SAR成像最顯著的特點就是合成孔徑時間急劇增加,電離層影響、天線振動、星歷誤差等非理想因素[5]對于GEO SAR成像或定位的影響十分嚴重。近年來,研究人員從理論上研究了GEO SAR長合成孔徑時間對地成像機理的可行性。但由于目前國內外在軌的SAR基本上屬于低軌合成孔徑雷達衛星,對GEO SAR成像的可行性僅能進行部分驗證。隨著北斗導航系統的發展,研究人員發現北斗3號傾斜地球同步軌道(Inclined Geo Synchronous Orbit,IGSO)與GEO SAR的軌道高度均在3 500 km左右,屬于高軌衛星,因此利用北斗3號IGSO導航衛星與地面靜止放置直達波接收站和地面反射波接收站的星地雙基構型[6-7]可以進一步對GEO SAR的長合成孔徑時間聚焦成像處理進行驗證。

雙基合成孔徑雷達成像中最為關鍵的問題就是收發雙置引起的同步問題[8]。新體制下由于距離向和方位向二維耦合的存在,以及要求在很長的積累時間內保持系統的相干性等,使得系統對同步問題的要求更加嚴苛。雙基合成孔徑雷達系統研究初期是通過在收發站之間建立光纖鏈路達到同步的效果[9],隨著衛星導航系統的發展,研究人員提出了一系列利用導航衛星鏈路實現時間同步的方法,例如雙向授時[10-11],衛星共視[12]等;LI等[13]在實驗中使用GPS對系統的收發站進行授時來實現收發站之間的時間同步;中國科學院大學為LT-1衛星設計了利用雷達信號之間的時間間隙來交換同步信號[14-16]的同步方案;上述研究方法均為增加同步信號傳輸鏈路以實現雙基SAR的時頻同步,同時隨著導航衛星的快速發展,該方法結合了導航衛星晶振高精度的優點,使得收發信號之間時頻同步精度非常高,但是此類方法均需要增加系統的硬件和同步鏈路,大大增加了成本。2007年德國開展的TerraSAR-X雙基實驗首次利用直達波處理完成了星機雙基合成孔徑雷達同步處理[17],之后出現了許多基于TerraSAR的成像及同步研究[18-20]。TerraSAR-X衛星的運行高度為514 km,相比于高軌衛星來說具有傳輸延時短、路徑損耗小以及信噪比高的優點,因此地面接收到的高軌衛星信號相比于低軌衛星受到的干擾因素更多,高軌雙基SAR的時頻同步算法也更加復雜。謝輝等[21]提出了一種通過使用最小二乘計算得到系統時間誤差進而完成對雙/多基地雷達系統的時間進行校準的方法,但該方法僅對雙基合成孔徑雷達成像中的時間同步誤差進行了校正,并沒有考慮頻率同步誤差對成像質量的影響。北京理工大學的張天[22]在其博士論文中提出了基于北斗2號導航衛星的時頻同步方法,并通過實測數據驗證了該方法的有效性,但該方法運用在北斗3號IGSO導航衛星星地雙基SAR成像過程中存在時間同步后,距離脈壓信號包絡分布與理論推導模型失配的問題。

針對上述問題,筆者提出了基于回波數據估計導航衛星測距碼延遲誤差,并進行誤差修正的星地雙基地合成孔徑雷達回波數據時頻同步誤差估計方法,并以北斗3號IGSO導航衛星為例,分析處理了雙基構型帶來的時頻同步誤差。該方法首先根據北斗3號IGSO導航衛星星歷和信號傳播延遲模型獲得直達波數據非理想采樣環境下正確的峰值位置序列來估計時間同步誤差和修正本地測距碼;然后通過修正后的本地測距碼對完成時間同步誤差校正的直達波數據進行脈沖壓縮來估計頻率同步誤差;最后利用估計的時頻誤差對反射波數據進行同步誤差校正和成像處理來實現目標的聚焦成像。實測數據的處理驗證了此方法的有效性。

2 星地雙基SAR系統空間觀測幾何模型

建立星地雙基合成孔徑雷達系統空間觀測幾何模型,如圖1所示。圖1中的坐標系為目標本體坐標系,坐標原點位于目標位置;Z軸為地心與目標的連線,方向指向目標;Y軸為IGSO導航衛星速度方向在Z軸垂直方向上的投影;X軸由右手螺旋定則確定。

圖1中PT(t)為IGSO衛星位置矢量;VT(t)為IGSO衛星速度矢量;AD為直達波接收天線位置矢量;AR為反射波接收天線位置矢量。實際外場試驗時,北斗3號IGSO衛星的運行高度約為35 786 km左右,而直達波接收天線與反射波接收天線間的距離為1 m左右,此時的空間大異構特性使得直達波通道和反射波通道之間信號的相位具有高相參性和時頻同步誤差的一致性[23],也即直達波信號和反射波信號的誤差特性可近似看成一樣的;利用此一致性可以消除實際成像處理時非理想因素對于成像的影響,最終得到一個聚焦良好的圖像。

圖1 星地雙基SAR系統空間觀測幾何模型

3 IGSO導航衛星星-地雙基SAR回波模型

3.1 發射信號模型

導航信號不同于常見的脈沖雷達合成孔徑雷達系統成像信號,如線性調頻(Linear Frequency Modulated,LFM)脈沖信號,導航信號為連續波信號,且導航信號中攜帶的導航信息會引入新的相位調制項。北斗3號導航系統現已發布的公開服務信號為B3I信號,它是由“測距碼+數據碼”采用二進制相移鍵控(Binary Phase Shift Keying,BPSK)調制在載波上構成的,其發射信號表達式為[24]

(1)

其中,j表示衛星編號;AB3I表示B3I信號振幅;CB3I表示B3I信號測距碼;DB3I表示B3I信號測距碼上調制的數據碼,數據碼是由“導航電文+Neumann-Hoffman”構成的;fc表示B3I信號載波頻率;φB3I表示B3I信號載波初始相位。筆者將針對北斗3號B3I信號進行分析,且各個編號的衛星的分析處理過程并無差別,故可將式(1)簡寫為

S(t)=AC(t)D(t)cos(2πfct+φ) 。

(2)

3.2 接收信號模型

回波信號是發射信號的時延,設τD、τR分別為直達波、反射波收發之間電磁波傳播路徑引起的時延,衛星星歷可通過導航電文解算得到,且電文中含有發射時間,因此接收信號的時延只與發射時間t有關。由圖1可知:

(3)

其中,c為光速;PT(t)表示衛星位置;AD表示直達波接收天線的天線相位中心位置;AR表示反射波接收天線的天線相位中心位置。

星地雙基合成孔徑雷達收發設備時間系統不一致使得發射通道與兩個接收通道之間存在時間同步誤差,地面接收端無法準確獲得發射信號的發射時間,即回波信號的起始采樣距離存在誤差;而收發本振不同則使得發射通道與兩個接收通道之間存在頻率同步誤差;發射通道與兩個接收通道之間的鐘差和頻差會造成后向投影(Back Projection,BP)成像結果散焦或者移位。假設Δt和Δfc分別為收發設備時間系統不一致和收發本振不同引起的信號誤差,則直達波、反射波基帶信號表達式為

(4)

(5)

其中,I(t)表示信號傳輸鏈路上非理想因素(電離層、對流層)引起的時延,可由電離層、對流層延遲改正預報模型[13]求得。

3.3 IGSO導航信號二維回波構建

進行星地雙基SAR成像處理之前,首先應將連續的導航信號按照一定的規則構建SAR成像所需的二維時域信號,也即對接收到的回波信號(4)和(5)進行二維分塊。(4)和(5)中測距碼周期為1 ms,為保留整個采集信號的完整性,可直接將連續的一維信號數據按照測距碼周期對應的脈沖重復頻率(Pulse Repetition Frequence,PRF)進行二維分塊,此時所有的回波數據都能被充分利用。接收信號(4)和(5)二維分塊后的時延和多普勒相位都是快慢時間的函數,但時延和多普勒相位隨時間的變化量存在明顯差異,時延可以近似認為是慢時間的矢量,而多普勒相位需要考慮隨快時間的變化情況[23]。

采用上述近似后,二維分塊后的直達波信號表達式為

(6)

其中,0≤η

(7)

其中,ΔTCA為接收端時間系統與發射端時間系統下測距碼的周期差,會造成一個隨慢時間η變化的時間徙動誤差。

由圖1可知,在目標本體坐標系下,北斗IGSO導航衛星的衛星高度遠遠大于AD的長度,衛星的運動速度也遠大于地表物體的運動速度,因此直達波和反射波二維時域信號中由衛星運動造成的多普勒頻偏近似可以認為一致,即

(8)

則二維分塊后的反射波信號表達式可寫為

(9)

4 IGSO導航衛星雙基地SAR系統的時頻誤差估計與補償方法

4.1 雙基地SAR系統時頻誤差估計與補償

直達波信號相比于反射波信號,具有斜距歷史簡單、信噪比高等優勢,且直達波和反射波同時接收的雙通道接收系統保證了兩個通道信號的相參性和時頻同步誤差的一致性。因此,筆者提出利用信噪比較高的直達波信號進行時頻同步誤差估計,并對反射波信號進行補償,然后再進行星地雙基合成孔徑雷達成像處理。北斗3號導航衛星測距碼模糊函數相比于線性調頻信號的模糊函數來說,對多普勒頻偏敏感度更高[25],因此成像處理時應當先進行距離向頻偏φD(η;τ)的移除,再進行距離向脈沖壓縮。在導航領域,一般使用信噪比較高的直達波信號通過捕獲跟蹤處理進行解擴來獲取距離向頻偏。實際信號采集環境中受到干擾信號或者硬件電路設備延遲的影響,去除距離向頻偏的直達波信號脈壓后的峰值位置序列仍然存在較多的干擾(如圖3去干擾前的峰值位置所示)。

因此,利用直達波信號進行時頻同步誤差估計首先應該獲取非理想采樣環境下正確的直達波脈壓后的峰值位置序列。由式(6)可知,非理想采樣環境下直達波距離脈壓后的峰值位置序列為

(10)

其中,mod(A,B)表示A除以B后的余數;Fs為信號的采樣率;Nr=Fs·T為距離向的點數;Pdis(η)為干擾導致的峰值位置偏移;從式(10)和式(7)可以看出,直達波的峰值位置Ppeak(η)與斜距歷史、常數項Δt、線性項η·PPRF·ΔTCA和傳播鏈路非理想延遲有關。在高軌SAR動輒幾百秒的合成孔徑時間里,衛星的運行軌跡不能再近似于直線運動使得衛星相對于地面的距離歷史表達式不再是簡單的線性關系,mod(A,B)運算和非理想因素誤差本身都是非線性的,因此線性估計方法無法獲取未受干擾的直達波距離徙動曲線。

為了解決上述的問題,提出了北斗3號雙基SAR基于直達波復信號的峰值位置序列恢復方法。以下將結合北斗3號導航衛星星地雙基試驗實測數據進行方法說明。

第1步 將地固系中的導航衛星星歷和直達波接收天線位置轉換至本體坐標系,計算斜距歷程和非理想因素引起的斜距誤差并對其進行多項式擬合后可表示為

RD(η)=|PT(η)-AD|+I(η)c≈Rlinear(η)+Rnon-linear(η) ,

(11)

其中,Rlinear(η)為多項式擬合后的線性分量,Rnon-linear(η)為2階及2階以上的非線性分量。

第2步 在Ppeak(η)項中移除Rnon-linear(η),從而得到一個線性分量如下:

(12)

(13)

設Na為方位向的采樣點數,N為常數;對該式進行NaN點傅里葉變換,可得

S(n)=FFFT(s(η),Na·N),n=0,1,2,…,NaN-1 。

(14)

將S(n)幅度的峰值位置記為nmax,則常數項和一次項系數估計值為

(15)

(16)

(a)

第5步 非理想采樣環境下的峰值位置序列Ppeak(η)去除非理想采樣環境干擾造成的位置偏移Pdis(η)得到直達波完整的峰值位置序列為

Ppeak_new(η)=Ppeak(η)-Pdis(η)≈

(17)

直達波實測數據去除非理想環境干擾前后的峰值位置序列如圖3所示。

圖3 去除非理想采樣環境干擾前后的峰值位置序列

將上述去除干擾后的直達波峰值位置矢量作為二維回波距離采樣起始位置,重新對一維直達波信號和反射波信號按照測距碼周期對應的PRF進行二維分塊、距離頻偏移除,可以保證直達波二維信號各距離向起點是與測距碼對齊的,直達波二維信號和反射波二維信號各距離向采樣起始距離是一致的,利用反射波和直達波之間的波程差可以實現反射波數據的時間同步處理。北斗3號導航衛星信號1個完整的測距碼上調制的數據碼D(η;τ)是定值,因此二維直達波和反射波信號表達式為

(18)

(19)

利用本地測距碼C(τ)對(18)、式(19)進行距離壓縮后可得

(20)

(21)

從式(20)可以看出重新劃分后的二維直達波信號利用本地測距碼C(τ)距離壓縮后的信號峰值位置均為1,由此保證同樣處理后的反射波信號式(21)逐方位向的起始采樣距離為導航衛星到直達波天線的斜距。反射波信號式(21)在進行BP成像時,可以按照各方位采樣時刻下,成像網格點到導航衛星和反射波接收天線的雙程斜距與導航衛星到直達波接收天線的單程斜距之間的斜距差獲得當前脈沖對成像網格點的貢獻,實現合成孔徑雷達聚焦成像處理。但在實測數據處理時,按照Ppeak_new(η)重新劃分后的直達波信號直接利用標準測距碼進行距離脈壓后信號的峰值位置不為1,如圖4(a)所示,這相當于利用反射波信號進行BP成像,距離脈壓后的信號的起始采樣距離未知,從而利用上述斜距差獲得當前脈沖對成像網格點的貢獻有誤,引起SAR圖像散焦。為了解決上述問題,進一步提出如下方法:

第1步 修正第二次進行距離壓縮的標準測距碼,修正后的表達式為

C′(η;τ)=C(τ-Pdis(η)) ;

(22)

第2步 利用修正后的測距碼C′(η;τ)對重新劃分的二維直達波信號進行距離壓縮后的表達式為

(23)

實測數據按照修正后的方法處理完之后直達波的峰值位置均為1,與式(20)達到的理論效果一致,如圖4(b)所示。

(a)

利用修正后的測距碼C′(η;τ)對重新劃分后的二維反射波信號進行距離壓縮后的表達式為

(24)

兩個接收通道共用一個本振信號保證了兩個通道信號相位的相參性,因此可以提取距離壓縮后的直達波信號峰值相位歷史用于反射波的頻率同步。距離壓縮后的直達波信號峰值相位歷史為

(25)

反射波信號相位去除直達波信號峰值相位可以實現頻率同步處理:

(26)

實現時頻同步后的反射波信號利用BiSAR實測數據處理最常用的BP算法進行聚焦成像。

4.2 雙基地SAR系統信號處理全流程

將對北斗3號IGSO導航衛星作為發射源,地面靜止接收反射波信號的星-地雙基地合成孔徑雷達系統進行聚焦成像處理,其成像處理流程如圖5所示。

圖5 基于時頻同步處理的星地雙基SAR成像處理流程

該流程包括以下步驟:

步驟1 (直達波峰值位置獲取) 利用直達波復數據和其解算得到的衛星星歷與多普勒頻偏,按照圖5左側處理流程得到正確的直達波峰值位置序列以及非理想采樣環境干擾引起的位置偏移Pdis(η)和衛星星歷等(如圖5虛線框所示);

步驟2 (時間同步) 根據正確的直達波峰值位置序列重新將導航接收信號(直達波、反射波)一維數據劃分為二維矩陣,實現時間同步;

步驟3 (頻偏移除) 重新劃分后的直達波、反射波信號進行頻偏移除;

步驟4 (距離壓縮) 根據干擾造成的位置偏移Pdis(η)重新更正本地測距碼后對直達波、反射波進行距離壓縮;

步驟5 (頻率同步) 提取直達波信號距離壓縮后的峰值相位矢量,并在反射波距離壓縮結果中進行移除,實現頻率同步;

步驟6 (BP成像) 對同步后的反射波數據在時域完成場景回波信號的BP成像。

5 北斗3號IGSO導航衛星實測數據成像

5.1 IGSO導航衛星星地雙基試驗配置

為了驗證時頻同步處理方法的可行性,將北斗3號IGSO導航衛星作為發射源獲取的BiSAR試驗實測數據進行成像處理。北斗導航信號到達地球表面時的信號強度非常小,約為-163 dBW[23],因此無法直接利用地面散射后的反射波信號進行處理得到合成孔徑雷達圖像。為實現北斗導航信號對地成像,在試驗場景中放置了有源轉發器將轉發器接收天線接收到的直達波信號進行提純放大后再由發射天線發出至反射波接收天線,轉發器所在位置可以相當于一個強散射點用于地面成像和分辨率的評估。提純性轉發器雖然解決了導航信號弱的問題,但是其處理過程也增加了反射信號的復雜度,例如碼延遲和鐘差,對此在外場試驗前將轉發器對回波錄取的影響做了一些考慮和處理,以使實驗符合實際的雙基散射情況。試驗中轉發器和直達/反射接收設備處各設置了一個相同的高穩定銣鐘,通過采用高精度原子鐘的方法為系統提供時間基準;成像系統工作前首先對兩個銣鐘完成了鐘差的修正;其次,對于導航反射接收信號進行了碼延遲測試,利用測試結果在信號成像處理之前對其碼延遲進行了校正。星地雙基試驗配置如圖6所示。

圖6 星地BiSAR試驗幾何構型與屋頂試驗光學照片

5.2 試驗結果

北斗3號IGSO導航衛星星歷參數與電離層延遲改正模型參數在導航電文中進行播發,因此,北斗3號IGSO導航衛星星-地雙基SAR直達波接收信號應首先進行預處理獲得星歷和參數。預處理包含捕獲、追蹤和解碼定位3個步驟。圖7 (a)為北斗3號IGSO導航衛星雙基SAR試驗捕獲的結果,由圖可知,捕獲到IGSO衛星編號(PRN)為6,8,9,13,16,38,39;圖7(b)為數據采集時間內可見衛星的天空圖。由天空圖可知,數據采集時間段內編號為6、9、16的IGSO衛星位于試驗場地的西南方位,編號為8、13的導航衛星位于東南方位。

(a)

考慮到雙基觀測幾何以及信噪比等因素,以編號為8的北斗3號IGSO導航衛星為例分析星地雙基SAR時頻同步誤差估計與補償后的成像結果,以驗證筆者所提方法的有效性。由于本次試驗中直達波接收天線與轉發器之間的距離較近,反射波與直達波信號中由電離層、對流層延遲引起的相位誤差可以不考慮,但是在成像場景較大的情況下,此誤差對于成像的影響不可忽略。8號星的數據所利用的時長為628 s,能夠完全驗證GEO SAR長合成孔徑時間成像方法的有效性。

圖8(a)為8號導航衛星反射波未進行時頻同步誤差補償的BP成像結果圖,該處理未能完成對于強散點的聚焦成像;圖8(b)為采用文獻[21]方法進行時頻同步誤差估計與補償后的BP成像結果圖,該方法在對時間同步處理后的信號進行脈沖壓縮時,未對距離向匹配濾波器(測距碼)進行修正,導致直達波距離脈壓后信號的包絡分布與理論推導模型失配,如圖4(a)所示,因此反射波信號進行BP成像時距離向起始采樣距離存在誤差,將會引起SAR圖像的散焦。圖8(c)為采用文中所研究方法對測距碼延遲誤差進行修正后再對回波信號進行頻率同步誤差估計與補償的BP成像結果圖,由圖可知采用此方法實現了強散點的聚焦成像;由圖8(a)(c)補償前后的成像結果可知雙基地SAR回波進行時域成像前需完成時頻同步誤差的估計與補償。由圖8(b)(c)的成像結果可知,導航信號的時頻同步在利用直達波峰值位置進行時間同步處理之后,需要將距離向匹配濾波器(測距碼)按照測距碼延遲誤差進行相應的修正以確保目標的聚焦成像。

圖8 時頻同步誤差校正前后8號衛星星地雙基合成孔徑雷達成像結果

表1 8號衛星強點目標成像指標

圖8(d) (e)為強散射點距離向和方位向的剖面圖,可見兩維壓縮結果均近似為sinc函數的形狀。由于采用文獻[21]方法進行時頻同步誤差估計與補償后的BP成像結果圖未能完成對于實測數據強散射點的聚焦成像,故無法將兩種算法的目標成像指標作對比說明。因此文中8號衛星強點目標成像指標僅對于文中所提算法實測數據成像結果的分辨率與雙基地幾何模型下的理論分辨率進行了對比,計算了分辨率在方位向與地距向的展寬比例,如表1所示。由表1可知實測數據的處理結果較接近理論值,有誤差存在的原因主要是北斗衛星的星歷數據為非精密軌道數據,且硬件存在一定的延遲且無法精確校正。8號導航衛星反射波成像結果驗證了星地雙基地幾何模型下長合成孔徑時間對地成像的可行性和本文時頻同步誤差估計與補償方法的正確性。

6 結束語

筆者研究了北斗3號IGSO導航衛星長合成孔徑時間下星地雙基合成孔徑雷達系統的時頻誤差估計方法。首先針對時間同步后的二維直達波距離脈壓信號包絡分布與理論推導模型失配的問題,提出了基于回波數據的測距碼時延誤差估計與修正方法;其次利用修正后的直達波脈壓信號峰值相位可用于反射波信號進行頻率同步,從而為反射波進行BP成像創造了條件;最后通過北斗3號IGSO導航衛星進行了成像試驗。根據時頻同步誤差補償前后以及測距碼修正前后的實測數據成像結果,驗證了所提方法的有效性。

猜你喜歡
信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
7個信號,警惕寶寶要感冒
媽媽寶寶(2019年10期)2019-10-26 02:45:34
孩子停止長個的信號
《鐵道通信信號》訂閱單
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
基于Arduino的聯鎖信號控制接口研究
《鐵道通信信號》訂閱單
基于LabVIEW的力加載信號采集與PID控制
Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: 国产丝袜无码一区二区视频| 一级一毛片a级毛片| 99在线视频免费| 亚洲日韩AV无码一区二区三区人| 日本手机在线视频| 国产三区二区| 欧美专区日韩专区| 国产丝袜第一页| 精品无码人妻一区二区| 在线观看亚洲人成网站| 国产欧美精品午夜在线播放| 中文字幕色在线| 久久99这里精品8国产| 欧美黑人欧美精品刺激| 久久99精品久久久久久不卡| 欧美视频在线第一页| 成人夜夜嗨| 伊人激情综合网| 国产jizz| 国产成人综合日韩精品无码首页| 成人伊人色一区二区三区| 性69交片免费看| 狼友视频一区二区三区| 亚洲激情99| 91久久国产综合精品女同我| 色婷婷在线播放| 日韩在线观看网站| 中文字幕av一区二区三区欲色| jizz在线观看| 人妻中文久热无码丝袜| 亚洲av成人无码网站在线观看| 伊人久综合| 国产9191精品免费观看| 一区二区三区在线不卡免费| 99视频有精品视频免费观看| 亚洲天堂网在线视频| 亚洲欧美一区二区三区图片| 成年人国产视频| 国产乱论视频| 成年午夜精品久久精品| yjizz视频最新网站在线| 真人免费一级毛片一区二区 | 精品国产成人高清在线| 激情成人综合网| 日韩A级毛片一区二区三区| 国产91麻豆免费观看| 国产资源免费观看| 91无码网站| 黄色三级网站免费| 拍国产真实乱人偷精品| 久久美女精品| 国产区成人精品视频| 亚洲国产理论片在线播放| 亚洲精品中文字幕无乱码| 色亚洲成人| 亚洲欧美日韩中文字幕在线| 亚洲AV无码乱码在线观看代蜜桃| 粉嫩国产白浆在线观看| 日本午夜精品一本在线观看| 亚洲天堂伊人| 国产在线观看成人91| 色视频国产| 一级看片免费视频| 成人国产免费| 日本伊人色综合网| 视频二区亚洲精品| 日韩欧美中文字幕在线韩免费| 凹凸国产熟女精品视频| 伊人久综合| 国产精品香蕉在线| 老司国产精品视频| 红杏AV在线无码| 日本免费精品| 真人免费一级毛片一区二区| 九色在线视频导航91| 日本亚洲欧美在线| 日韩第八页| 国产91小视频| 国产9191精品免费观看| 99热国产这里只有精品无卡顿"| 欧美黄色网站在线看| 国产精品视频导航|