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

被動微波遙感觀測資料干擾對地表參數反演的影響分析

2017-09-21 01:19:04吳瑩錢博王振會
自然資源遙感 2017年3期

吳瑩, 錢博, 王振會

(南京信息工程大學氣象災害教育部重點實驗室/氣候與環境變化國際合作聯合實驗室/氣象災害預報預警與評估協同創新中心/中國氣象局氣溶膠與云降水重點開放實驗室,南京 210044)

被動微波遙感觀測資料干擾對地表參數反演的影響分析

吳瑩, 錢博, 王振會

(南京信息工程大學氣象災害教育部重點實驗室/氣候與環境變化國際合作聯合實驗室/氣象災害預報預警與評估協同創新中心/中國氣象局氣溶膠與云降水重點開放實驗室,南京 210044)

以東亞陸地為研究區,利用2011年7月1—16日AMSR-E (advanced microwave scanning radiometer - earth observing system ) 二級亮溫數據,采用一維變分反演收斂度量識別法對研究區的無線電頻率干擾(radio-frequency interference,RFI)進行識別,并進一步分析了該干擾對星載被動微波數據反演地表參數的影響。研究表明,在東亞地區,AMSR-E C波段和X波段都普遍存在RFI,大多出現在工業區、科研中心、人口密集的大城市、機場和高速公路等區域; 一般情況下,C波段和X波段的干擾信號分布區域基本不重合; 水平和垂直極化通道都有較強的RFI存在,其時間上具有持續性,且不同極化方式下干擾信號的強度因地而異; 干擾信號的強度隨衛星地球方位角的變化而變化,當星載微波輻射計掃描到某一地球方位角度范圍內時這些視場才會受到干擾。研究還發現,用受RFI影響的微波數據反演得到的地表參數值誤差較大。因此,在采用被動微波觀測數據反演地表參數之前,有必要先有效地檢測和剔除RFI。

被動微波遙感; 地表參數; 無線電頻率干擾(RFI); 東亞陸地

0 引言

地表溫度[1-3]、土壤水分[4]、積雪[5]和植被等地表參數是水文模型、氣候及陸面過程模式中的重要參數[6-8]。傳統方法所測結果只代表觀測點的局部特征,而遙感可以快速同步獲取大面積區域的地表參數[6],提供二維陸面分布信息。因此利用衛星遙感數據反演地表參數,探討衛星資料的反演理論及其實際應用方法,已經成為遙感科學的一個重要研究領域[6]。

根據傳感器通道采用的波段不同,可將衛星傳感器分為3類: 可見光、紅外和微波傳感器。由于較可見光和紅外遙感而言,微波遙感不易受大氣影響,具有全天時、全天候的監測能力以及對云、雨和大氣較強的穿透能力,并且微波傳感器對于植被特性的變化、地表土壤水分和積雪參數十分敏感,微波數據已被廣泛應用于積雪、植被等地表參數的監測和反演應用之中[6]。近年來廣泛使用的星載微波傳感器有: 搭載于EOS/Aqua衛星上的先進微波掃描輻射計(advanced microwave scanning radiometer-earth observing system,AMSR-E)、搭載于美國國防部Coriolis衛星上的全極化微波輻射計(WindSat)、搭載于我國風云三號衛星上的微波成像儀(microwave radiation imager,MWRI)和搭載于日本GCOM-W1衛星上的AMSR-E的繼任者(AMSR-2)等。

然而,這些用于反演地表參數的低頻觀測資料均不同程度地受到地面無線電頻率的干擾。星載被動微波傳感器的無線電頻率干擾(radio-frequency interference,RFI)往往是由地面主動微波傳感器的發射信號或陸面反射輻射信號產生的,很容易覆蓋地表產生的相對較弱的熱發射輻射信號,使得星載被動微波傳感器接收的信息不能真實地反映地表狀況。如果不能準確地將其識別和剔除,往往會導致較大的反演誤差,從而顯著降低現有以及將來的被動微波資料的使用效率[9-13]。

本文采用一維變分反演(one dimensional variational retrieval,1-DVAR)收斂度量識別法分析了AMSR-E觀測資料在東亞陸地上RFI的時空分布特征、形成原因以及RFI對微波數據反演地表參數的影響,為提高星載微波資料在陸面過程模式及資料同化中的利用率提供依據。

1 數據源

本研究中,采用的被動微波遙感觀測資料是2011年7月1—16日AMSR-E的二級亮溫數據。WindSat,MWRI和AMSR-2等設備都具有和AMSR-E頻率相近的通道。AMSR-E搭載于2002年美國國家宇航局發射的地球觀測系統Aqua衛星上,提供6.925 GHz,10.65 GHz,18.7 GHz,23.8 GHz,36.5 GHz和89.0 GHz 6個頻率的水平和垂直2種極化方式,共12個通道的微波觀測值,主要用于觀測大氣、陸地、海洋和冰圈的參數。AMSR-E是圓錐型掃描輻射計,天線圓錐掃描角為47.4°,掃描幀幅寬度為1 445 km。星下點空間分辨率隨觀測頻率而變化,從56 km(6.925 GHz)到5.4 km(89.0 GHz)不等。AMSR-E資料每16 d為一個周期,每一個周期內覆蓋完全相同的區域,因此本文選用一個完整周期的AMSR-E二級亮溫數據。

在一維變分反演過程中,需要提取NCEP (NOAA’s national centers for environmental prediction) GDAS (global data assimilation system)客觀分析場中的相關參數作為代價函數迭代的初始場。GDAS系統每天生成4個時次(00UTC,06UTC,12UTC 和18UTC)、水平空間分辨率為1°×1°的各種大氣和地表參數場。

2 RFI識別算法

基于衛星微波和紅外觀測資料,采用1-DVAR不僅可以反演大氣參數(如大氣溫和濕垂直廓線)、云參數(云量和云頂高度)外,還可以反演某些地表參數(如地表溫度和地表發射率等)。變分反演的基礎是在觀測場與背景場誤差均服從高斯誤差分布的假定條件下,再對所定義的代價函數求最小化,得到最小誤差的分析場。其代價函數一般可以表示為[14]

(1)

式中:J(X)為代價函數;X為被反演的大氣(或地表)狀態變量;X0為大氣(或地表)狀態的先驗信息(稱作背景場向量);Ym為已獲得的觀測資料;B為背景場誤差協方差矩陣;E為觀測場誤差協方差矩陣;H為前向算子,代表一種具有某種復雜結構的從模式空間向觀測空間的映射。

通過對J(X)求導,并使其為0,求得目標函數式(1)達到最小值的解,即

(2)

可以得到

(X-X0)=△X=[(B-1+HTE-1H)-1HTE-1]×[Ym-H(X0)],

(3)

將式(3)用于迭代循環,一直到J(X)達到最小時,結束循環,從而可以得到使目標函數式(1)達到最小值的解,稱作一維變分產生的分析場(Xa),即

J(Xa)=minJ(X)。

(4)

對于衛星微波數據反演地球物理參數而言,這個前向算子H就是前向模式,即輻射傳輸模式。本文采用通用輻射傳輸模式(the community radiative transfer model,CRTM)作為1-DVAR中的前向模式。CRTM[15]由美國衛星資料同化聯合中心開發,適用于各種天氣條件,可以模擬所有微波頻率下由冰晶、雪晶、雨滴、霰粒和云中液態水等產生的散射,并生成所有大氣、地表參數相應的輻射值和輻射梯度(即雅可比矩陣)[14]。

在1-DVAR計算過程中,需要輸入GDAS 的大氣參數垂直廓線(溫度、濕度和云中液水含量)與地表參數(地表溫度、土壤濕度和植被覆蓋度等)作為反演的背景場。由于衛星微波觀測輻射對地面輻射率敏感性很高,因而將地表發射率作為預反演的大氣狀態變量的組成部分,在求代價函數極小值的過程中動態地調整地表發射率,有助于更好地反演地表參數。1-DVAR計算可得到地表和大氣參數2類,其中地表參數主要有地表溫度和地表發射率。

反演計算中,收斂度量是用前向算子最后一次模擬的亮溫值和測量值之間所有殘差的均方根。計算公式如下

χ2=[Ym-H(X)]T×E-1×[Ym-H(X)] 。

(5)

χ2作為是否可以達到收斂的判據,也用來衡量前向模式的優劣。只有選入并在1-DVAR計算中有效使用的那些通道的測量值才用于計算這個度量。通常,當χ2≤1時,認為可以達到收斂。然而,可以根據實際情況把這個標準放寬到10,即認為χ2>10時,反演結果不可靠。反演過程中發現,1-DVAR算法中的收斂度量值和RFI信號的強弱有著極強的相關性,χ2的值越大,意味著該處的RFI信號越強[16]。

3 結果與分析

AMSR-E 2011年7月10日一維變分法反演過程中計算χ2的分布如圖1所示。

圖1 2011年7月10日東亞陸地1-DVAR收斂度量分布(升軌)

從圖1中可以看出,在東亞地區χ2值較大的區域主要分布在中國的京津冀地區、長江中下游地區,日本的北海道部分地區、太平洋沿岸和瀨戶內海沿岸的狹長地帶以及沿海平原地區等。RFI主要來源于城區及工業區的無線電通信、雷達以及有線電轉播信號等,和該地區人口密集的城市和工業區分布基本一致。

大多數情況下,陸地上(非冰雪覆蓋區)的發射和散射主要取決于土壤含水量、地表粗糙度、地形地貌、地表溫度和植被覆蓋等。由于土壤中水分和植被的吸收作用隨頻率增加而加大,地表亮溫隨頻率增加而呈上升趨勢。在AMSR-E的低頻率通道,大氣相對透明,輻射計接收到的微波輻射的光譜特性主要由地表的發射和散射決定。隨著頻率的增加,地表和植被的散射效應也增加,成為亮溫降低的一個因素。而當頻率低于30 GHz時,散射效應通常是有限的。因此,AMSR-E低頻率通道的RFI是最可能造成C波段通道亮溫值遠遠高于X波段、X波段遠遠高于K波段的原因,從而在低頻波段上導致負頻譜梯度的產生。

可以通過各個通道數據之間的聯系(頻譜差異),用RFI指數來量化RFI的強度。RFI指數的定義[9]如下

RFIf1,p=TBf1,p-TBf2,p,

(6)

式中:TB表示亮溫度; 下標p表示水平h或垂直v極化方式;f1和f2表示2個相鄰的頻率(f10,RFIf1,p數值越大意味著f1頻率通道的RFI強度越大。

由于χ2值和RFI信號的強弱有著高度相關性[16],可以利用線性擬合的方法來擬合出χ2值和各通道RFI強度的關系曲線。在本文研究區域內,獲得線性關系如下

χ2=2.422 81+0.321 787RFI6,v+0.061 115 8RFI6,h,

(7)

χ2=3.569 67+0.089 525 1RFI10,v+0.285 240RFI10,h,

(8)

式(7)和式(8)線性擬合的標準差分別為1.366 7和1.246 5。顯而易見,RFI的強度越大,對應的χ2值也越大。

相對于中國東部而言,RFI在日本分布尤為顯著。日本的工業區主要分布在太平洋沿岸和瀨戶內海沿岸的狹長地帶,這里分布著京濱工業區、名古屋工業區、阪神工業區、瀨戶內海沿岸工業區和北九州工業區等5大日本傳統工業區。新興工業的布局主要取決于人才、交通和環境因素,日本的九州島(“硅島”)、筑波(“科學城”)等新興工業區主要分布在九州南部和本州島的東北部,依托于現代飛機場和四通八達的高速公路網。在工業區內,日本的RFI除了受無線電通信的影響,還要考慮到衛星的影響,如阪神工業區周圍就分布著10多個衛星城,也就是說這里的RFI受衛星以及無線電通信的綜合影響。

由于AMSR-E是圓錐型掃描儀器,其數據中提供了衛星的地球方位角數值,定義為衛星掃描方向相對于觀測視場正北面的方位,即觀測視場和衛星的連線在地球上的投影與正北方向的夾角,取值范圍為[-180°,180°],順時針方位為正值,反之為負值。根據Aqua衛星的軌道路徑,選取圖1中2個區域(B1: N 35°~37°,E 139°~142°,B2: N 38.5°~41°,E 140.5°~143°)進一步討論RFI出現位置和強度與儀器的地球方位角之間的關系,如圖2和圖3所示。

(a) 6.925 GHz水平極化 (b) 6.925 GHz垂直極化 (c) 10.65 GHz水平極化(d) 10.65 GHz垂直極化

圖22011年7月1—16日B1區域亮溫的觀測方向和強度間的關系

Fig.2ObservationdirectionandthemagnitudeofbrightnesstemperatureofAMSR-EinB1duringJuly1-16,2011

(a) 6.925 GHz水平極化 (b) 6.925 GHz垂直極化 (c) 10.65 GHz水平極化(d) 10.65 GHz垂直極化

圖32011年7月1—16日B2區域亮溫的觀測方向和強度間的關系

Fig.3ObservationdirectionandthemagnitudeofbrightnesstemperatureofAMSR-EinB2duringJuly1-16,2011

對于B1區域,6.925 GHz和10.65 GHz的4個通道均不同程度地出現了異常高值的亮溫點,說明該地區在C波段和X波段都存在RFI。C波段垂直極化時的RFI比水平極化時的RFI分布更廣,強度也更大。而且,6.925 GHz的2個通道幾乎在所有地球方位角范圍內均出現RFI,而10.65 GHz的2個通道只在一定地球方位角范圍內([-180°,-150°],[110°,180°])才出現RFI,且10.65 GHz比6.925 GHz的 RFI強度更大。對于B2區域,6.925 GHz和10.65 GHz 2個頻率也均受RFI影響,且水平極化時的RFI較垂直極化時強。6.925 GHz的2個通道出現RFI對應的衛星地球方位角范圍([-60°,-30°],[110°,150°])比10.65 GHz通道的RFI分布更廣,說明6.925 GHz的RFI的方向性沒有10.65 GHz的RFI方向性顯著。

以NCEP GDAS的客觀分析場作為1-DVAR算法的背景場,利用2011年7月10日AMSR-E的亮度溫度數據進行反演計算地表參數。反演結果如圖4所示,其中(a)—(d)分別是6.925 GHz 和10.65 GHz水平、垂直極化時的亮溫值分布。圖4(e)和(f) 分別是1-DVAR計算得出的地表溫度和地表發射率分布。

(a) 6.925 GHz水平極化 (b) 6.925 GHz垂直極化(c) 10.65 GHz水平極化

(d) 10.65 GHz垂直極化(e) AMSR-E反演的地表溫度 (f) AMSR-E 反演的地表發射率(6.925 GHz垂直極化)

圖42011年7月10日東亞陸地AMSR-E6.925GHz和10.65GHz亮溫和反演的地表溫度、地表發射率

Fig.4AMSR-Ebrightnesstemperatureat6.925GHz,10.65GHzandAMSR-EretrievedlandsurfacetemperatureandlandsurfaceemissivityovereasternAsiaonJuly10,2011

從圖4中可以發現,亮溫值分布中顯示存在大量高異常亮溫值的點(紅色區域),而在地表溫度和地表發射率分布中顯示有斷斷續續的白色區域及一些零散分布的地表溫度偏高的紅色區域。主要由于該區域存在強RFI使得變分反演的收斂度量值過大,從而得不到反演值,造成反演失敗(白色區域); 較弱的RFI使得反演值偏差較大,呈現許多地表溫度反演異常偏高的孤立點(紅色區域),與周圍地表溫度的分布不連續。而通常情況下,自然地表發射輻射特性往往呈現連續、平滑的分布特征。因此,如果反演之前,不對使用的微波資料進行RFI檢測,將會大大降低反演精度,特別是強RFI的存在大大降低了微波資料的使用效率(即得不到反演值)。

4 結論

本文基于2011年7月16 d的AMSR-E二級亮溫數據,用1-DVAR收斂度量識別法檢測出東亞陸地上RFI信號的時空分布和方向特性,并分析了其產生的原因及其對地表參數反演的影響,得出如下結論:

1)以東亞陸地為例,用1-DVAR收斂度量識別法可檢測RFI的位置和強度。

2)該研究區內,AMSR-E的RFI主要分布在工業區、科研中心、人口密集的大城市、機場和高速公路周邊。

3)研究區除個別地區同時存在6.925 GHz和10.65 GHz的RFI信號,AMSR-E的RFI在C波段和X波段的分布區域重合較少,大部分位于不同的地理位置,其時間上具有持續性。

4)總體上,RFI呈現出較強的方向性,RFI的強度隨衛星的地球方位角變化而變化,當星載微波輻射計掃描到某一地球方位角度范圍內時這些視場才受RFI的影響。

5)RFI會使受影響區域地表參數的反演值偏高,嚴重影響參數的反演精度。故有效地進行RFI 檢測并進行校正方法的開發是提高星載被動微波遙感資料利用率和地表參數反演精度的必要進程。

[1] 周芳成,宋小寧,李召良.地表溫度的被動微波遙感反演研究進展[J].國土資源遙感,2014,26(1):1-7.doi:10.6046/gtzyyg.2014.01.01. Zhou F C,Song X N,Li Z L.Progress of land surface temperature retrieval based on passive microwave remote sensing[J].Remote Sensing for Land and Resources,2014,26(1):1-7.doi:10.6046/gtzyyg.2014.01.01.

[2] 劉晶,馬紅章,楊樂,等.基于被動微波的地表溫度反演研究綜述[J].遙感技術與應用,2012,27(6):812-821. Liu J,Ma H Z,Yang L,et al.A survey of surface temperature retrieval by passive microwave remote sensing[J].Remote Sensing Technology and Application,2012,27(6):812-821.

[3] 毛克彪,施建成,李召良,等.用被動微波AMSR數據反演地表溫度及發射率的方法研究[J].國土資源遙感,2005,17(3):14-17.doi:10.6046/gtzyyg.2005.03.04. Mao K B,Shi J C,Li Z L,et al.The land surface temperature and emissivity retrieved from the AMSR passive microwave data[J].Remote Sensing for Land and Resources,2005,17(3):14-17.doi:10.6046/gtzyyg.2005.03.04.

[4] 鮑艷松,毛飛,閔錦忠,等.基于FY-3B/MWRI數據的裸土區土壤濕度反演[J].國土資源遙感,2014,26(4):131-137.doi:10.6046/gtzyyg.2014.04.21. Bao Y S,Mao F,Min J Z,et al.Retrieval of bare soil moisture from FY-3B/MWRI data[J].Remote Sensing for Land and Resources,2014,26(4):131-137.doi:10.6046/gtzyyg.2014.04.21.

[5] 孫知文,于鵬珊,夏浪,等.被動微波遙感積雪參數反演方法進展[J].國土資源遙感,2015,27(1):9-15.doi:10.6046/gtzyyg.2015.01.02. Sun Z W,Yu P S,Xia L,et al.Progress in study of snow parameter inversion by passive microwave remote sensing[J].Remote Sensing for Land and Resources,2015,27(1):9-15.doi:10.6046/gtzyyg.2015.01.02.

[6] 施建成,杜陽,杜今陽,等.微波遙感地表參數反演進展[J].中國科學(地球科學),2012,42(6):814-842. Shi J C,Du Y,Du J Y,et al.Progresses on microwave remote sensing of land surface parameters[J].Science China Earth Science,2012,55(7):1052-1078.

[7] 張廷軍,晉銳,高峰.凍土遙感研究進展:被動微波遙感[J].地球科學進展,2009,24(10):1073-1083. Zhang T J,Jin R,Gao F.Overview of the satellite remote sensing of frozen ground:Passive microwave sensors[J].Advances in Earth Science,2009,24(10):1073-1083.

[8] 李芹,鐘若飛.模擬AMSR-E數據的地表多參數反演[J].國土資源遙感,2011,23(1):42-47.doi:10.6046/gtzyyg.2011.01.08. Li Q,Zhong R F.Multiple surface parameters retrieval of simulated AMSR-E data[J].Remote Sensing for Land and Resources,2011,23(1):42-47.doi:10.6046/gtzyyg.2011.01.08.

[9] Li L,Njoku E G,Im E,et al.A preliminary survey of radio-frequency interference over the U.S. in Aqua AMSR-E data[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(2):380-390.

[10]Njoku E G,Ashcroft P,Chan T K,et al.Global survey and statistics of radio-frequency interference in AMSR-E land observations[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(5):938-947.

[11]Lacava T,Coviello I,Faruolo M,et al.A multitemporal investigation of AMSR-E C-band radio-frequency interference[J].IEEE Transactions on Geoscience and Remote Sensing,2013,51(4):2007-2015.

[12]Zou X L,Zhao J,Weng F Z,et al.Detection of radio-frequency interference signal over land from FY-3B microwave radiation imager(MWRI)[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(12):4994-5003.

[13]官莉,張思勃.星載微波輻射計歐洲大陸無線電頻率干擾分析[J].光學學報,2014,34(7):0728004. Guan L,Zhang S B.Source analysis of radio-frequency interference over Europe land from advanced microwave scanning radiometer-E[J].Acta Optica Sinica,2014,34(7):0728004.

[14]Boukabara S A,Garrett K,Chen W C,et al.MiRS:An all-weather 1DVAR satellite data assimilation and retrieval system[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(9):3249-3272.

[15]Weng F Z,Han Y,van Delst P,et al.JCSDA Community radiative transfer model(CRTM)[C]//Proc.14th Int.ATOVS Study Conf.2005:217-222.

[16]Adams I S,Bettenhausen M H,Gaiser P W,et al.Identification of ocean-reflected radio-frequency interference using WindSat retrieval chi-square probability[J].IEEE Geoscience and Remote Sensing Letters,2010,7(2):406-410.

(責任編輯:陳理)

Effectofradio-frequencyinterferenceonthelandsurfaceparametersretrievalfrompassivemicrowaveremotesensingdata

WU Ying, QIAN Bo, WANG Zhenhui

(KeyLaboratoryofMeteorologicalDisaster,MinistryofEducation(KLME)/JointInternationalResearchLaboratoryofClimateandEnvironmentChange(ILCEC)/CollaborativeInnovationCenteronForecastandEvaluationofMeteorologicalDisasters(CIC-FEMD)/KeyLaboratoryforAerosol-Cloud-PrecipitationofChinaMeteorologicalAdministration,NanjingUniversityofInformationScienceandTechnology,Nanjing210044,China)

Radio-frequency interference (RFI) over eastern Asia land was detected and analyzed using one dimensional variational retrieval (1-DVAR) convergence metric method from AMSR-E (the advanced microwave scanning radiometer - earth observing system) Leval 2A measurements during July 1-16, 2011. And then its influence on the retrieval of surface parameters was studied. It is found that the RFI signals are detected both at C and X band channels of AMSR-E over eastern Asia, and the signals are most densely concentrated in industrial zones, scientific research centers, metropolises, airports and highways. Moreover, RFI signals at C and X band normally do not coincide with the same distribution area. AMSR-E RFI over eastern Asia land exists along both horizontal and vertical polarization channels. Furthermore, the intensity of AMSR-E RFI varies with the earth azimuth angle of the satellite; measurements are contaminated by RFI only when the spaceborne microwave radiometer is within some earth azimuth angle range. Lastly, it is also found that retrieved land parameters have large deviations from RFI contaminated microwave measurements. Therefore, it is expected to detect even weakened RFI effectively prior to retrieving land surface parameters from passive microwave remote sensing measurements.

passive microwave remote sensing; land surface parameter; radio-frequency interference (RFI); eastern Asia land

10.6046/gtzyyg.2017.03.26

吳瑩,錢博,王振會.被動微波遙感觀測資料干擾對地表參數反演的影響分析[J].國土資源遙感,2017,29(3):176-181.(Wu Y,Qian B,Wang Z H.Effect of radio-frequency interference on the land surface parameters retrieval from passive microwave remote sensing data[J].Remote Sensing for Land and Resources,2017,29(3):176-181.)

2015-12-29;

2016-03-02

國家自然科學基金項目“FY-3微波數據RFI訂正及我國典型地區地表微波發射率反演研究”(編號: 41305033)、江蘇省基礎研究計劃—青年基金項目“微波地表溫度計算及其對中國典型地區地表發射率反演改進研究”(編號: BK20150911)和江蘇高校優勢學科建設工程資助項目(PAPD)共同資助。

吳瑩(1980-),女,講師,博士,主要從事大氣探測與大氣遙感方面的教學和研究工作。Email: wuying_nuist@163.com。

P 422.2; TP 722.6

: A

: 1001-070X(2017)03-0176-06

主站蜘蛛池模板: 亚洲色图在线观看| 精品欧美一区二区三区在线| 成人午夜视频免费看欧美| 2022国产无码在线| 国产精品999在线| 国产玖玖视频| 天天摸夜夜操| 91久久夜色精品| 国产91精选在线观看| 色天堂无毒不卡| 香蕉蕉亚亚洲aav综合| 国产成人亚洲无吗淙合青草| 亚洲第一区在线| 久久精品aⅴ无码中文字幕| 欧美国产在线精品17p| 亚洲国产精品日韩欧美一区| 国产99视频在线| 亚州AV秘 一区二区三区| 色综合久久88| 国产成人一级| 国产精品深爱在线| 欧美日本中文| 最新亚洲人成网站在线观看| 国产人成在线观看| 色爽网免费视频| 51国产偷自视频区视频手机观看| 国产在线无码av完整版在线观看| 国产精品欧美激情| 91偷拍一区| 久久人体视频| 亚洲综合天堂网| 亚洲国产中文精品va在线播放 | 一本色道久久88综合日韩精品| 亚洲婷婷六月| 一区二区欧美日韩高清免费 | 日韩一区精品视频一区二区| 99精品伊人久久久大香线蕉| 久久婷婷人人澡人人爱91| 日韩成人免费网站| 黄色片中文字幕| 亚洲最新地址| 日韩 欧美 国产 精品 综合| 欧美激情视频一区| 国产在线麻豆波多野结衣| 国产色爱av资源综合区| 激情乱人伦| 91亚瑟视频| 中文字幕在线观| 免费在线看黄网址| 狼友视频国产精品首页| 波多野结衣在线se| 国产在线精品人成导航| 欧美成人综合视频| 色偷偷一区| 91啦中文字幕| 欧美综合中文字幕久久| 国产一级做美女做受视频| 91成人在线免费观看| 国产精品无码AV片在线观看播放| 午夜欧美在线| 亚洲成av人无码综合在线观看| 日韩精品视频久久| 日韩视频福利| 伊人久热这里只有精品视频99| 亚洲欧美自拍中文| 夜夜操狠狠操| 91探花国产综合在线精品| 亚洲三级色| 99热这里只有精品在线观看| 91亚洲精选| 亚洲香蕉伊综合在人在线| 亚洲成人网在线播放| 中文无码伦av中文字幕| 91娇喘视频| 欧美人与性动交a欧美精品| 亚洲天堂在线视频| www亚洲精品| 亚洲AV无码乱码在线观看裸奔| 国产国产人在线成免费视频狼人色| 成人蜜桃网| 国产偷国产偷在线高清| 日韩成人高清无码|