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

臨近空間風溫遙感干涉儀設計及正演

2020-09-05 12:02:52何微微武魁軍王后茂
光學精密工程 2020年8期
關鍵詞:大氣

何微微,武魁軍,傅 頔,王后茂,李 娟

(1.武漢科技大學 城市學院,湖北 武漢 430083;2.中國科學院 武漢物理與數學研究所,湖北 武漢 430071;3.中國科學院 西安光學精密機械研究所,陜西 西安 710119;4.中國科學院 國家空間科學中心,北京 100190)

1 引 言

臨近空間通常指高度介于20~100 km的大氣層區域[1]。臨近空間風溫探測對于大氣物理和空間科學的發展具有重要的學術意義。精確的風場與溫度場探測對于研究中層大氣的半年振蕩[2]、大氣潮汐結構的季節變化[3]、全球大氣動量和能量的傳播,以及球面大氣中行星波的激發具有重要的科學意義[4-5]。 此外,大氣風場溫度場的精準探測,對于提高衛星發射的成功率、保障航空航天安全性、提高中長期天氣預報的準確率具有十分重要的應用價值。

臨近空間的探測設備非常豐富,但由于技術限制,大氣風場、溫度場的探測主要局限在地基遙感探測技術,如流星雷達[6]、拉曼測溫激光雷達[7]、多普勒測風激光雷達[8]、氣輝成像干涉儀[9]等。隨著空間技術的不斷發展,提出了具有全球觀測尺度的星載風溫探測手段。星載氣輝成像干涉儀采用臨邊觀測模式,具有垂直分辨率高、光譜分辨率高、相位靈敏度高的優勢,已成為衛星遙感領域的前沿課題[10]。

星載成像干涉儀技術通常采用多普勒頻移法獲得風速,采用多普勒展寬法獲得溫度,即通過干涉圖的相位變化獲取氣輝的多普勒頻移信息以反演風速,通過干涉圖對比度的變化獲取氣輝多普勒展寬信息以反演溫度。由于干涉圖的相位變化與對比度變化相比更敏感,風場反演精度高于溫度場反演精度。然而,對處于熱力學平衡態的雙原子或多原子分子,可以采用雙線比率法(即兩條分離的輻射線的相對強度)確定大氣溫度。雙線比率法可提高溫度測量精度(約1~3 K)。1991年,搭載于上層大氣研究衛星UARS的風場探測干涉儀-風成像干涉儀(WIND Imaging Interferometer,WINDII)與高分辨力多普勒成像儀(High Resolution Doppler Imager,HRDI)分別采用了多普勒展寬法、雙線比率法兩種測溫方式。WINDII利用O原子的紅線(630.0 nm)和綠線(557.7 nm),測量了海拔75~320 km的風速(精度5~8 m/s)和溫度(精度20~75 K)[11]。HRDI利用O2分子電子態躍遷 的(0,0),(1,0),(2,0)3個振動帶的吸收及輻射光譜[12],獲得了平流層(10~40 km)到中間層與低熱層(50~120 km)的大氣風場(精度為5 m/s)和溫度場(精度為7 K)[13]。

WINDII和HRDI取得成功之后,星載風溫探測吸引了越來越多的關注,其發展趨勢是將探測高度向較低海拔區域推進以獲得更廣泛的高度覆蓋范圍。加拿大 York大學研制的同溫層風場輸運干涉儀(Stratospheric Wind Interferometer for Transport studies,SWIFT) 采用Michelson干涉成像儀探測O3位于1 133.433 5 cm-1附近的振轉光譜目標譜線[14]。2015年,Shepherd團隊提出了SWIFT-DASH的改進方案,采用多普勒外差干涉儀(Doppler Asymmetric Spatial Heterodyne,DASH)替代Michelson干涉儀[15]。多普勒調制型分子濾光器(Doppler Wind and Temperature Sounder,DWTS)采用分子濾光技術,可同時獲取多普勒頻移和輻射線寬信息,從而實現25~250 km的連續風溫探測[16-17]。由于目標譜線均采用了中波紅外(Mid-wavelength Infared,MWIR)和長波紅外(Long-wavelength Infared,LWIR)輻射線的設計策略,SWIFT和DWTS的優勢具備全天時風溫探測的能力。然而,工作在MWIR或LWIR波段的星載紅外遙感器系統通常會受到來自光學元件熱背景的強烈干擾,測量不確定性較大。鑒于此,紅外遙感系統平臺需要較低的制冷溫度,因此增加了平臺成本,但SWIFT的溫度測量能力依舊相對較低。此外,DWTS的模擬結果表明,在50~100 km的高度內,DWTS無法獲得垂直軌道方向的風速分量;100~250 km高度內,雖然可以獲得垂直軌道方向分量的風速,但其不確定度是沿軌方向風速分量的10倍。

根據加拿大中層大氣模型[18],一項科學研究表明,為了改善天氣預報系統中的數據同化分析精度,40~60 km高度內,風速精度需達到5~10 m/s,溫度精度需達到5~10 K[19]。針對臨近空間的探測需求,可利用O2分子的紅外大氣帶(躍遷)1.27 μm波段的氣輝為目標譜線。鑒于O2(a1Δg)氣輝的輻射強度相對較強、自吸收能力相對較弱,適用于45~90 km高度內的風溫探測。加拿大StaSci計劃部署的中間層成像邁克爾遜干涉儀(Mesospheric Imaging Michelson Interferometer,MIMI)以及NASA的MIDEX計劃支持的波成像邁克爾遜干涉儀(Waves Michelson Interferometer,WAMI)均采用該譜帶的強、弱兩組氣輝譜線(每組3條)進行風溫探測[20]。中國科學院武漢物理與數學研究所的武魁軍等人利用O2分子1.27 μm近紅外氣輝對臨近空間風溫探測精度進行了研究[21-23],并分析了氫氧氣輝對1.27 μm氣輝風溫探測精度的影響[24]。

本文提出了臨近空間風溫探測干涉儀(Near-space Wind and Temperature Sensing Interferometer,NWTSI)儀器模型,通過廣角邁克爾遜干涉儀觀察1.27 μm波段的O2(a1Δg)氣輝,實現臨近空間風場、溫度場的星載遙感探測。與WAMI和MIMI類似, NWTSI重點關注低海拔區域(40~80 km)。NWTSI采用邁克爾遜干涉儀,由多普勒頻移法獲得風速,由雙線比率法獲得溫度,具有較高的風溫探測精度。本文介紹了NWTSI儀器模型,分析了大氣臨邊輻射光譜的特性并對NWTSI進行了正演仿真模擬,最后討論了風溫探測的不確定度。

2 儀器模型

NWTSI是一種基于廣角邁克爾遜干涉儀的星載風溫遙感探測載荷。它通過臨邊觀測的方式,在兩個相互正交的觀測方向上交替觀測。交替觀測由指向鏡將視場(Field of View,FOV)指向相應的觀測方向實現,視場角為1.5°,高度覆蓋20~120 km。NWTSI 前后兩個視場可先后觀測到相同的大氣區域,在410 km的高度上觀測間隔大約為9 min。

圖1為NWTSI儀器模型的光學系統示意圖。該儀器參考了WAMI的儀器設計,由兩個望遠鏡、一個邁克爾遜干涉儀、一個窄帶濾光片、一個法布里-珀羅(F-P)標準具以及一個近紅外焦平面組成。視場由望遠鏡1(telescope 1)和視場光闌(FS1)決定。入射光由指向鏡進入望遠鏡1,然后入射到邁克爾遜干涉儀。望遠鏡1的角放大率為2,則邁克爾遜干涉儀的視場為3°×3°。邁克爾遜干涉儀長壁(LPA)的反射鏡分為厚度不同的四等份,每一分區較上一分區光程差(Optical Path Difference,OPD)增加λ/4,則四分區之間的階躍光程差為λ/4,λ/2和3λ/4。四分區的設計可實現4個干涉圖樣的同時采集,與WINDII的順序采集相比,該方法有效避免了測量過程中強度變化引起的誤差。邁克爾遜干涉儀短壁(SPA)的反射鏡安裝在電容式位置傳感器控制的壓電器件上,則干涉儀鏡子也可以步進,用于實現飛行期間的校準。平面鏡M3的作用是將光學元件折疊成緊湊的形狀。入射光經視場光闌FS2后進入望遠鏡2。望遠鏡2的角放大率為0.5,因此,濾光系統的視場再次變為1.5°×1.5°。邁克爾遜有視場展寬的能力,受視場角影響小,而F-P干涉儀透射函數卻受視場角的影響很大,該設計的主要意圖是有效解決F-P干涉儀視場角受限的問題。F-P標準具與干涉濾光片(F1與F2)共同作用實現目標輻射線的分離,同時阻止散射太陽光的進入。望遠鏡2聚焦在F-P標準具后面的金字塔棱鏡上。棱鏡的邊緣與邁克爾遜干涉儀長壁反射鏡的邊緣部分對齊,使得入射光投影到焦平面陣列(Focal Plane Array,FPA)的不同區域,每個區域對應不同的階躍相位。棱鏡后設置有分束鏡(BS),將光路分為兩個通道:強線通道與弱線通道。表1為NWTSI儀器的系統參數。

圖1 臨近空間風溫探測干涉儀的光學原理圖Fig.1 Optical schematic of near-space wind and temperature sensing interferometer (NWTSI)

表1 臨近空間風溫探測干涉儀的系統參數

3 仿 真

正演仿真的目的是模擬儀器的功能和效果,產生預期的干涉圖像。通過測量仿真可以實現儀器性能的獨立分析,并可對儀器性能做出評估。正演模型包含大氣輻射模塊,邁克爾遜干涉儀模塊,濾波器模塊,以及光學系統、傳感器陣列、紅外焦平面等設備的系統參數。

3.1 正演仿真

3.1.1 大氣輻射

1.27 μm氣輝具有輻射強度強、高度覆蓋范圍廣的優勢,是全球大氣風溫遙感的最佳目標譜線之一[25]。1.27 μm氣輝主要來源于O2第一激發態O2(a1Δg)到基態O2(X3Σg)的躍遷。實際電子躍遷過程中,振動態和轉動態的存在,使得1.27 μm氧氣輝成為包含約150條輻射線的光譜帶。

O2(a1Δg)氣輝臨邊輻射光譜由沿視向的路徑積分計算獲得。積分過程中大氣被分成許多離散的薄層,每一層都被視為是均勻的,逐層計算的輻射傳輸方程表示為:

L(v)l=L(v)l-1exp(-∑α(v)l,iul)+

(1)

其中:J(ν)l,i與α(ν)l,i分別是第i條臨邊觀測路徑、第l層大氣的輻射源函數與吸收系數;ul=nl(ZU-ZL)是第l層大氣的O2分子的分子柱密度;ZU與ZL是第l層大氣上、下邊界的臨邊切點高度。

圖2 O2紅外大氣波段在40,60 km切點高度上的輻射光譜與透射譜Fig.2 Limb radiance and transmittance spectra of O2 infrared atmospheric band at tangent heights of 40 km and 60 km

圖2描述了臨邊切點高度為40 km和60 km,不考慮大氣風場與衛星速度時O2(a1Δg)紅外大氣帶的臨邊輻射和透射光譜。圖中紅線代表臨邊輻射光譜,藍線代表臨邊透射光譜(彩圖見期刊電子版)。可見,處于大氣帶中心位置的輻射光譜具有較強的輻射強度,但兩翼位置處的輻射光譜自吸收效應較弱,并且輻射強度與吸收系數均隨切點高度的變化而變化。因此,兩翼位置處的輻射光譜適合低海拔區域的風溫探測,而中間區域的輻射光譜適合高海拔區域。

在眾多的輻射光譜中選擇最佳的目標譜線是星載風溫成像干涉儀研制的關鍵之一。Ward等人提出了系統的選線標準[21],即目標譜線的選取需遵循三個原則:溫度靈敏度高,譜線分離性好和高度覆蓋范圍廣。NWTSI目標譜線的選取參考了WAMI與MIMI的選線方案,選用強線組與弱線組共6條譜線,每組包含3條譜線。強線組譜線波數分別為:7 908.973 6 cm-1,7 909.653 5 cm-1,7 911.013 3 cm-1;弱線組譜線為:7 821.110 5 cm-1,7 822.221 6 cm-1,7 822.952 5 cm-1。20~60 km的低空區域存在較強的自吸收現象,選用弱線來減少自吸收對反演的影響;60~120 km的高空區域體輻射率相對較低,信號也相對較弱,但高空的自吸收效應幾乎不存在,可選用強線探測,增大信噪比,提高測風精度。強、弱線的聯合使用能夠覆蓋20~120 km的臨近空間區域,既避免了低空自吸收效應對風溫探測帶來的影響,又保證了高空風溫探測的精度。

3.1.2 濾波器

O2分子的1.27 μm輻射帶包含數量眾多并且分布十分密集的譜線。為了實現目標譜線的分離,濾光系統需要有大約0.1 nm的帶寬。因此,需將F-P標準具與窄帶濾波器聯合使用。此外,窄帶濾波器還具備降低背景光亮度、提高條紋對比度的作用。

NWTSI的濾波系統由F-P標準具和窄帶濾波器(F1和F2)組成。系統的總濾波函數由二者共同決定。F-P標準具的帶通濾波函數采用Airy 函數模型,如式(2)所示,設計參數如表1所示。

(2)

其中:σ是入射光的波數,r是標準具的反射率,netal是標準具材料的折射率,d是標準具材料的厚度,θ是標準具腔內的折射角。

圖3 弱線組與強線組的最佳透射函數Fig.3 Optimal transmittance functions of weak group and strong group

采用雙線比率法獲得溫度信息,NWTSI要在同一視場中同時觀察3條譜線。為了降低成本,只使用一個標準具,其后利用分束鏡將光線分離為弱線、強線兩個通道。圖3給出了強線組與弱線組的最佳透射率函數。圖中S1: 7 908.97 cm-1,S2: 7 909.65 cm-1,S3: 7 911.01 cm-1;W1: 7 821.11 cm-1,W2: 7 822.22 cm-1,W3: 7 822.95 cm-1。

3.1.3 邁克爾遜干涉儀

邁克爾遜干涉儀的光程差隨離軸角的變化如下[25]:

(3)

其中:nL與nS分別是邁克爾遜干涉儀長臂和短臂的折射率;i是像素離軸角;tL與tS分別是干涉儀長臂和短臂的長度。邁克爾遜干涉儀的光學參數如表1所示。將式(3)應用于每個像素,可得光程差的像素分布如圖4所示。

圖4 邁克爾遜干涉儀的光程差像素分布Fig.4 Pixel distribution of optical path difference of Michelson interferometer

步進相位(k-1)/2π,(k=1,2,3,4),通常指的是視場幾何中心位置處的步進相位。每個像素第k次步進的相位可表示為:

(4)

其中d是邁克爾遜干涉儀的步進尺寸。

利用邁克爾遜干涉儀的相位,以及濾波器的帶通濾波函數,結合圖2中的臨邊光譜輻射特性,即可獲得NWTSI的干涉圖像。

3.2 測量仿真

圖像的像素值由大氣輻射傳輸方程(式(1))干涉圖決定。給定像素的干涉圖可表示為[25]:

[1+Uljcos(2πvΔlj+φklj]dv,

(5)

其中:Iklj為第l行第j列第k次步進的干涉圖,flj(v)是相對總濾波函數,Ulj是儀器可見度,Δlj是光程差,φklj是邁克爾遜干涉儀第k次的步進相位,ν是波數,Rlj是儀器響應度。Rlj可定義為[22]:

(6)

其中:AΩ是光學系統的光學擴展量,t是干涉時間,q是探測器的量子效率,h是普朗克常量,c是真空光速,τ是濾波器與光學系統的透過率,v0是O2輻射線的中心波長。

大氣風場的精準測量需要較大的視場,從而獲得較大的響應度。增加入射光立體角Ω的尺寸可以增大響應度,但這會使得光程差增加,引起條紋對比度的降低。因此,傳統儀器在立體角的擴大上有一定的局限性。NWTSI儀器臨邊觀測視場為1.5°×1.5°,第一臺望遠鏡的角放大率為2,則邁克爾遜干涉儀的視場為3°×3°。根據視場展寬原理,邁克爾遜干涉儀可以保證在大光學擴展量的情況下實現高分辨率。然而,F-P標準具的立體角則受限于:

ΩFP=2πn2/RFP,

(7)

其中RFP=σ0/δσFP,為F-P標準具的分辨率。因此,望遠鏡2的角度放大率設計為0.5,則F-P標準具的視場再次變為1.5°×1.5°。

WINDII的測量方式采用逐像素測量的四點采樣法,即通過比較四步進所得的干涉圖像獲得多普勒風引起的干涉圖相移,最終反演得到風速。NWTSI的四分區設計則可以實現4個干涉圖樣的同時采集,并通過金字塔形棱鏡同時成像在FPA的4個不同區域。通過比較FPA不同區域干涉圖的相移即可實現大氣風場的探測。

圖5所示為強線組與弱線組4個階躍相位的干涉圖像的正演仿真。如圖5所示,干涉圖樣與像素位置、階躍相位、濾波透過函數在視場上的變化、以及臨邊輻射光譜隨切點高度的變化均有相關性。干涉圖樣中的每個像素代表對應空間區域臨邊輻射光譜的積分強度。NWTSI對應的空間區域尺寸約為1 km × 1 km。像素間干涉圖樣的差異主要來自于臨邊觀測氣輝強度隨切點高度的變化,以及光學系統對干涉圖樣的調制。正演仿真過程中,氣輝光譜輻射只考慮垂直方向的變化,忽略視場相對于地平線的傾斜以及地球的曲率,并假設儀器可見度和響應度對所有像素均相同。

圖5 四個階躍相位的正演仿真干涉圖像Fig.5 Forward simulation of interference images for four phase steps

探測器第l行第j列像素的干涉圖總相位可以通過四點采樣法得到。總相位減去地球自轉相位、衛星運動相位以及儀器相位即可得到風場引起的相位。與地基觀測技術相比,臨邊觀測儀器的光學擴展量要小的多,然而,由于臨邊觀測具有較長的積分路徑,使得臨邊輻射譜的強度遠大于地基觀測的強度。

大氣風場、溫度場測量的精度受到信噪比(Signal-to-Noise Ratio,SNR)和臨邊觀測權重的共同影響。在干涉圖像的模擬中考慮了3種主要的噪聲源,即散粒噪聲、讀出噪聲以及探測器暗噪聲(相關數據如表1所示)。強線組與弱線組的SNR和臨邊觀測權重隨高度的變化如圖6所示。由圖可見,強線組的SNR峰值位于40~45 km處,在該區域外緩慢下降,且在小于40 km區域下降得更快。強線組的SNR在較高的海拔區域(約60 km)達到峰值。低海拔區域O2(a1Δg)態的碰撞淬滅作用,導致態密度衰減,從而使得臨邊觀測權重隨切點高度的減小而減小。

圖6 信噪比與臨邊觀測權重隨高度的變化Fig.6 Signal-to-noise ratio (SNR) and limb-view weight as functions of altitude

4 誤差分析

依據光程差的變化,像素干涉圖可以表示為截斷傅里葉級數。傅里葉系數J1,J2和J3(也稱為表觀量)與干涉圖相關[26]。風和溫度的測量不確定度可以通過誤差傳遞來確定。風和溫度可以由表觀量J1,J2和J3導出。視向風可由J2和J3計算得到,大氣溫度由J1確定。

(8)

4.1 風速誤差

風速vw可以由干涉圖的相移δφ獲得[25]:

(9)

其中:c代表真空光速,相位φ可以式(10)計算[22]:

(10)

每一個像素都可以計算得到視向風以及隨機誤差。視向風的隨機方差為[26]:

(11)

圖7 強線組(藍線)、弱線組(綠線)輻射線以及聯合譜線(紅線)多普勒風測量的隨機誤差標準差Fig.7 Random error standard deviation of Doppler wind by strong(blue line) and weak groups(green line) and combined lines (red line)

圖7所示為強線組與弱線組兩組輻射線風速測量的隨機誤差標準差(彩圖見期刊電子版)。由圖可見,對于低海拔區域(低于42 km)的風速探測,使用弱輻射線更加精確,而高海拔區域使用強線更有優勢。這是由于強輻射線在低空具有強烈的自吸收效應(如圖1所示)。

4.2 溫度誤差

大氣溫度由溫度敏感性不同的兩條獨立輻射線的吸收率的比值決定[27],即有:

(12)

式中:E″是躍遷的低態能級的能量;S0(T)為參考溫度T0=296 K時的線強度,A是輻射線的吸收率積分,h與kB分別是普朗克常量與玻爾茲曼常量。

溫度的不確定度由線強度和吸收率積分的不確定度決定,可表示為[28]:

(13)

吸收率積分的不確定度主要來自干涉圖平均值J1的信噪比。線強度不確定度是由固有不確定度ΔS(T0)與溫度相關的不確定度ΔST(T)的傳遞引起的。ΔS(T0)可以由數據庫獲得[28],ΔST(T)與線強度成正比:

(14)

式中Q(T)是分子配分函數。

圖8所示為強線組與弱線組溫度探測的隨機誤差標準偏差分布(彩圖見期刊電子版)。由圖可見,弱線組在切點高度40~75 km區域的隨機誤差水平為1.5~2 K,強線組的誤差約為2~3 K。高海拔區域由于信號較弱,誤差水平較大。測溫精度除了受信噪比影響以外,還與躍遷的低態能級的能量差有關,弱線組測溫精度高于強線組的原因在于它的低態能級的能量差比強線組大。

圖8 強線組(藍線)與弱線組以及聯合譜線(綠線)以及聯合譜線(紅線)溫度探測的隨機誤差標準偏差分布Fig.8 Temperature error profiles by using three emission lines in strong(blue line) and weak groups(green line)and combined lines (red line)

5 結 論

本文設計了臨近空間風溫遙感干涉儀,以1.27 μm附近的O2(a1Δg)氣輝為目標源,采用廣角邁克爾遜干涉儀,實現了臨近空間區域風場、溫度場的同時探測。首先通過逐行積分算法,并結合最新光化學反應模型(包含光譜參數、速率常數和太陽通量),開發了O2(a1Δg)氣輝臨邊輻射光譜的輻射傳輸模型。與WINDII和HRDI等其他風溫干涉儀不同,NWTSI利用兩組強度不同的目標譜線,僅對臨近空間區域進行風溫探測。按照譜線選取的原則選擇了7 820~7 824 cm-1內的3條譜線組成弱線組,7 908~7 912 cm-1內的3條譜線組成強線組。弱線用于低空探測,以避免自吸收效應對測量結果的影響;強線用于高空探測,以實現高的測量精度。此外,建立了由大氣輻射傳輸模塊,邁克爾遜干涉儀模塊,濾波器模塊,以及光學系統、傳感器陣列、紅外焦平面等設備的系統參數組成的正演模型,并計算得到了強線組與弱線組的正演仿真干涉圖。通過正演模型獲得了臨邊觀測圖像,并對風速及溫度的測量不確定度進行了計算分析。風溫探測的精度受到了信噪比和臨邊觀測權重的共同作用。由于低海拔區域O2(a1Δg)態的碰撞淬滅效應,臨邊觀測權重隨高度的降低而降低,導致低海拔區域探測精確度的降低。模擬結果表明,風速測量精度為1~3 m/s,溫度測量精度為1~3 K,水平分辨率沿視線約350 km,在40~80 km高度內沿軌道約170 km,滿足臨近空間風溫探測精度的要求。

猜你喜歡
大氣
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
首次發現系外行星大氣中存在CO2
科學(2022年5期)2022-12-29 09:48:56
宏偉大氣,氣勢與細膩兼備 Vivid Audio Giya G3 S2
太赫茲大氣臨邊探測儀遙感中高層大氣風仿真
有“心氣”才大氣
如何“看清”大氣中的二氧化碳
學生天地(2020年18期)2020-08-25 09:29:24
大氣穩健的美式之風Polk Audio Signature系列
稚拙率真 圓融大氣
中國篆刻(2017年3期)2017-05-17 06:20:46
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
主站蜘蛛池模板: 国产成人亚洲毛片| 亚洲一区毛片| 亚欧成人无码AV在线播放| 精品久久久久久成人AV| 91视频首页| 国产精品亚洲精品爽爽| 中文字幕乱码二三区免费| 国产成人一区| 久久综合伊人77777| 久久精品亚洲中文字幕乱码| 亚洲成人动漫在线| 国产97视频在线观看| 国产一在线观看| 天堂av高清一区二区三区| 亚洲天堂色色人体| 亚洲成a人片7777| 一级片免费网站| 亚洲第一黄片大全| 97se亚洲综合不卡| 国产 在线视频无码| 97se综合| 成人欧美日韩| 中文字幕佐山爱一区二区免费| 波多野结衣AV无码久久一区| 国产成人a在线观看视频| 久草网视频在线| 亚洲香蕉伊综合在人在线| 国产又粗又爽视频| 四虎永久在线视频| 亚洲人网站| 色噜噜综合网| 草逼视频国产| 无码中文字幕乱码免费2| 欧美一级片在线| 欧美日韩一区二区三区在线视频| 国产Av无码精品色午夜| 国产在线八区| 日韩 欧美 小说 综合网 另类| 日韩专区第一页| 国产午夜一级淫片| 尤物视频一区| 在线色国产| 国产麻豆福利av在线播放| 黑人巨大精品欧美一区二区区| 久久久精品国产亚洲AV日韩| 日本草草视频在线观看| 91精品免费久久久| 欧美人与牲动交a欧美精品 | 青草娱乐极品免费视频| 全裸无码专区| 毛片免费观看视频| 国产青榴视频在线观看网站| 国产亚洲精品91| 无码精油按摩潮喷在线播放| 日本亚洲成高清一区二区三区| 成人综合网址| 亚洲无码高清一区二区| 91国内视频在线观看| 欧美区一区二区三| 欧美色视频在线| 91精品啪在线观看国产60岁| 这里只有精品在线| 热99精品视频| 欧美翘臀一区二区三区 | 在线观看无码av五月花| 日韩视频免费| 玖玖精品视频在线观看| 国产视频大全| 亚洲av无码人妻| 国产v精品成人免费视频71pao| 动漫精品啪啪一区二区三区| 免费毛片网站在线观看| 色综合久久综合网| 综1合AV在线播放| 久久这里只有精品66| 国产又粗又猛又爽| 欧美成人手机在线观看网址| 无码专区国产精品一区| 国产欧美视频一区二区三区| 中文字幕丝袜一区二区| 国产黄在线观看| 国产精品不卡片视频免费观看|