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

基于HJ-1A星HSI數據的植被吸收光合有效輻射APAR估算

2017-05-11 11:22:12常雅軒于穎崔綠園侯吉宇曾德偉
森林工程 2017年2期
關鍵詞:模型研究

常雅軒,于穎,崔綠園,侯吉宇,曾德偉

(東北林業大學 林學院,哈爾濱 150040)

基于HJ-1A星HSI數據的植被吸收光合有效輻射APAR估算

常雅軒,于穎*,崔綠園,侯吉宇,曾德偉

(東北林業大學 林學院,哈爾濱 150040)

為了更準確地測量植被對太陽輻射的有效吸收量,本文對基于HJ-1A衛星的HSI數據進行植被吸收光合有效輻射APAR的研究,為陰陽葉的吸收有效光合輻射量的分離計算提供基礎,并為森林生態系統凈初級生產力計算提供準確的參數數據。研究的主要步驟有以下三步:首先對一景包含漠河氣象站的HSI數據進行預處理,其中包括HSI數據的格式轉換,去條帶處理以及大氣校正等;在此基礎上應用逐日太陽輻射模型計算日總太陽輻射,并通過NDVI計算得到光合有效輻射吸收比率FPAR;最后用FPAR與太陽輻射共同計算出APAR。最終得到處理后的圖像,并且驗證了逐日太陽輻射模型和研究區域氣象站數據的相關性,通過NDVI和SR的結合計算出FPAR,精度得到了提高;植被吸收光合有效輻射APAR計算完成后,將灰度圖像顯示為專題圖像。從專題圖中可以得出不同地區森林的APAR值,并最終推算出全年的APAR。

HSI數據處理;太陽輻射模型;FPAR;APAR

0 引言

吸收光合有效輻射(Absorbed Photosynthetically Active Radiation,APAR)指植被冠層吸收的參與光合生物量積累的光合有效輻射部分[1]。森林生態系統凈初級生產力在研究碳循環和全球氣候變化中具有重要地位,綠色植物冠層吸收光合有效輻射是CASA等碳循環模型中計算凈初級生產力的基礎參數,因此,準確估算區域APAR對研究區域森林生態系統的光合作用、凈初級生產力和碳循環具有十分重要的意義,并且研究APAR是進一步研究植被陰陽葉APAR分離與計算的基礎[2]。最初,由于輻射觀測點極少,獲得APAR的觀測資料較少,人們主要運用氣候學計算方法,包括理論和經驗方法,獲得任一地區的APAR分布情況,但具有依賴于地面的氣候觀測資料、計算方法參數多、計算過程復雜的缺陷,后來,一部分學者開始嘗試利用遙感信息來反演大范圍地區的APAR分布。

Hipps等研究了冬小麥抽穗前后APAR與LAI的函數關系[3];Meeree等研究得到植物的光合作用產物與APAR呈近似線性關系[4];Chen(1999)根據兩葉模型,將冠層植被分為陰葉和陽葉,分別計算陰陽葉的太陽輻射,可用于陰陽葉APAR的分離與計算[5]。在國內,張喜旺等利用MODISNDVI數據研究了2010年河南省APAR分布[6];李剛等基于MODIS數據研究了內蒙古地區APAR及其時空變化特征[7];陳新芳等研究對于大范圍APAR的監測和估算,通過對PAR和基于遙感數據計算的FPAR來實現[8]。本研究基于高光譜遙感數據分三步估算APAR。高光譜遙感數據具有高光譜分辨率,被廣泛應用于植被長勢檢測及植被生物化學信息提取,環境資源減災HJ-1A衛星搭載超光譜成像儀(Hyper-SpectralImager,簡稱HSI),其數據具有高光譜分辨率、較短的時間分辨率和較低的信噪比,且對廣大用戶免費下載,是本文研究大興安嶺漠河地區很好的遙感數據源。

1 研究區域概況及數據

研究區域是黑龍江省大興安嶺地區,該地區位于東經121°12′~127°00′、北緯50°00′~53°33′。東部緊挨小興安嶺,西部和內蒙古緊密相連,南部與松嫩平原交界,北部與俄羅斯只有一江之隔。大興安嶺地區平均海拔573 m,年平均氣溫-2.6℃,年平均降水量428.5~526.8 mm,該地區與全國其他地區均不相同,是典型的寒溫帶大陸性季風氣候,全年無霜期80~110 d,冰封期180~200 d。太陽輻射數以百萬計共4 500耳,年日照2 600 h,年累積溫度2 100℃。

HJ-1A星HSI數據根據其校正程度共分為0~5級產品,分別為原始數據產品、輻射校正產品、系統幾何校正產品、幾何精校正產品、高程校正產品和標準鑲嵌圖像產品[9]。普通用戶通過訪問中國資源衛星應用中心的網站就可以免費獲得HSI的2級產品,本文采用HSI的2級產品。研究選取了包含122°22′E,53°28′N漠河氣象站的黑龍江省大興安嶺地區的一景圖像,該景圖像的產品號是1049632,圖像的詳細信息見表1。

表1 實驗數據的詳細參數

2 研究方法

本研究的主要流程包含三步:首先對一景包含漠河氣象站的HSI數據進行預處理,包括HSI數據的格式轉換,去條帶處理以及大氣校正等;第二,在此基礎上應用逐日太陽輻射模型計算日總太陽輻射,并通過歸一化植被指數(NDVI)和比值植被指數(SR)計算得到光合有效輻射吸收比率FPAR;最后用FPAR與太陽輻射共同計算出APAR。流程如圖1所示。

圖1 技術路線圖Fig.1 Technology roadmap

2.1 HSI數據預處理

HSI數據的預處理包括圖像旋轉、去條紋、大氣校正和MNF噪聲分離變換等過程。HJ-1A星 HSI數據 2 級產品存在像元灰度值比較小、與周圍有明顯差異的垂直條帶。條帶噪聲的存在嚴重影響了圖像的質量和應用,使有用波段信息不能完全被人們獲得和利用。因此,本文選取了全局去條紋法對條紋進行去除(圖2)[10]。HSI數據去條紋處理時認為圖像是按列垂直排列的,所以在進行去條紋處理前要對圖像先進行旋轉。之后,對進行條帶處理后的HSI數據采用ENVI的FLAASH模塊進行大氣校正[11]。

圖2 去條紋處理前后圖像對比Fig.2 Image comparison before and after stripe dealing

對比大氣校正前后的光譜曲線可知(圖3),去條紋的噪聲處理和大氣校正后原始圖像在近紅外波段的反射峰值有所提高,突出了植被特征,在藍光處出現了吸收谷,消除了大氣的多次散射影響。最后,對處理前后的圖像做MNF噪聲變換,圖像處理后圖像信息主要集中的波段噪聲較少。

2.2 逐日太陽輻射的模擬

太陽輻射是地球表層上的物理、生物和化學過程的主要能源,也是生態系統過程模型、水文模擬模型和生物物理模型研究中的必要參數[12]。本文選取大興安嶺研究區域中的漠河氣象站,將該氣象站的氣象數據和輻射數據作為逐日太陽輻射模型中緯度、日照時數和實測日太陽輻射的輸入數據,對逐日太陽輻射進行模擬。

太陽輻射是太陽直接輻射和散射輻射之和,時間尺度不同太陽總輻射值不同。利用逐日太陽輻射模型可以計算出一天之中的太陽總輻射,一天當中太陽總輻射與當天日長和日照時數有關,逐日太陽輻射利用Angtrom-Prescott方程計算[13]:

H=HL×(a+b×S/SL)。

(1)

式中:H為日實測太陽總輻射;HL為晴天狀態下的地面總輻射;S為日照時數;SL為日長;a和b為兩個經驗參數,本文引用童成立等利用實測總輻射和日照百分率月平均值等資料計算得到的適合我國大部分地區的經驗參數a為0.248,b為0.752[13]。

圖3 圖像處理前后植被光譜曲線對比Fig.3 Vegetation spectral curve comparison beforeand after image processing

理想狀態下(晴天狀態)大氣上空太陽輻射通過大氣層到達地面的輻射量可以表示為:

HL=0.8×HO。

(2)

式中:HL為理想狀態下(晴天狀態)的地面總輻射;HO為大氣上空的太陽總輻射。

大氣上空的太陽總輻射HO計算公式:

HO=(1/π)·Gsc·EO·(cosΦ·cosδ·sinWs+(π/180)·sinΦ·sinδ·Ws。

(3)

式中:Gsc為太陽常數,其值一般為1 367 W/m2(相當于118.108 MJ/(m2·d));EO為地球軌道偏心率校正因子;Φ為緯度;δ為太陽赤緯;Ws為時角(角度)。

2.3 光合有效輻射

到達地球表面的全部太陽輻射只有一部分特定波段的太陽輻射能被綠色植物用來進行光合作用,這部分能量被稱作光合有效輻射,其中光合有效輻射量與到達地表的太陽總輻射有關,利用公式計算:

PAR=H×ηPAR。

(4)

式中:PAR為光合有效輻射;H為日實測太陽總輻射;ηPAR為光合有效輻射占太陽輻射的比值,不同方法得到的ηPAR有差異,國內學者估算PAR時多采用ηPAR為0.49~0.50,本文采用Szeice測定的 為0.5[14]。

2.4 植物吸收光合有效輻射

植物吸收光合有效輻射與植被類型和植被覆蓋狀況有關,用植物光合有效輻射吸收率來表示植物對光合輻射的吸收效率,研究表明光合有效輻射吸收率與NDVI和植被類型有關,本文APAR的估算是基于CASA模型的,計算APAR是以像元為單位逐一計算的,植物吸收光合有效輻射可以表示為[15]:

APAR=PAR×FPAR。

(5)

式中:APAR為植物吸收光合有效輻射;PAR為光合有效輻射;FPAR為光合有效輻射吸收率,利用植物吸收光合有效輻射計算公式就可以計算出每一個像元的植物吸收光合有效輻射,從而確定整個研究區域的總的植物吸收光合有效輻射。

2.5 光合有效輻射吸收比例FPAR的計算

2.5.1 歸一化植被指數NDVI計算FPAR

植被對太陽有效輻射的吸收比例FPAR取決于植被類型和植被覆蓋情況,由遙感數據獲得的歸一化植被指數FPAR能很好地反映植物覆蓋狀況[16]。在一定范圍內,FPAR與NDVI之間存在著線性關系[17],這一個關系可以根據FPAR的最大值、最小值以及該種植物類型對應的NDVI的最大值、最小值來確定。由于FPAR的最大值最小值與植被類型無關,分別取0.95和0.001[18],所以具體公式如下:

(6)

2.5.2 比值植被指數SR計算FPAR

進一步研究表明,FPAR與比值植被指數SR也存在著較好的線性關系[19],可由以下公式表示:

(7)

其中,FPAR的最大值和最小值分別取0.95和0.001,SR(x,t)由公式(8)確定:

(8)

由于本研究區域為大興安嶺區域,該區域的植被類型主要為落葉針葉林,所以根據朱文泉[20]對中國各陸地植被類型NDVI和SR的最大值最小值列出的分布表(表2)查得本式中NDVImax為0.738,NDVImin為0.023,SRmax為6.63,SRmin為1.05。

表2 中國陸地植被類型NDVI和SR的最大值與最小值

3 結果與分析

3.1 逐日太陽輻射的結果驗證

選取了大興安嶺地區的漠河氣象站進行比較分析,其輻射觀測臺站地理信息見表3。

表3 輻射觀測臺站地理信息

根據逐日太陽輻射模型模擬計算了大氣上界和地面晴空條件下的太陽輻射日變化,逐日太陽輻射值選取了2013年漠河氣象站全年記錄的日太陽實測輻射值。結果表明模擬計算得到的晴空條件下太陽輻射值與實際觀測到的逐日太陽輻射最大值基本一致(圖4),少數奇異點(327.66 MJ/m2)為缺測數據,由于大氣上界的太陽輻射穿過大氣層時,大氣中的各種成分會對太陽輻射有吸收、散射作用,最終導致太陽輻射的衰減,使得到達地球表面的太陽輻射值小于大氣上界的太陽輻射值。

圖4 漠河氣象站地區大氣上界、地面晴空條件下和實際觀測太陽輻射日變化Fig.4 Observed diurnal variation of solar radiation under atmospheric boundary and ground clear-sky conditions in Mohe Weather Station area.

利用逐日太陽輻射模型計算得到漠河氣象站一年中每日的太陽輻射模擬值與動態監測的太陽輻射真實值,進行回歸分析(圖5),模型模擬結果與漠河觀測站點結果有很高的相關性(表4),決定系數(R2)達到了0.95,平均偏差MBE為-1.06 MJ/(m2·d),表明逐日太陽輻射模擬值偏小,并且大都小于實際觀測的太陽輻射值,均方根誤差RMSE反映了模擬值與實際觀測值間的偏離程度,為1.94 MJ/(m2·d),其結果小于童成立等對黑龍江漠河的研究結果[21],以太陽輻射預測值為橫坐標,殘差為縱坐標,得到殘差圖(圖6),圖中可以看出描繪的點圍繞殘差等于0的直線上下隨機散布,表明變量實測太陽輻射與模擬太陽輻射有顯著的線性相關關系,進而表明逐日太陽輻射模型可用于模擬漠河氣象站地區的日太陽輻射值,實現了由漠河氣象站的點輻射值,驗證太陽輻射模型并應用到一景遙感圖像的面輻射值,可用于APAR的計算,但是日照時數和日序數作為逐日太陽輻射模型的輸入參數,在時間步長為旬、月和年的情況下,模擬太陽輻射值計算的工作量勢必增大。

圖5 太陽日輻射模擬值與真實值比較Fig.5 Comparison between simulated and actual values of solar radiation

圖6 模擬太陽輻射的殘差圖Fig.6 Simulated residual plots of solar radiation

站點R2MBE(MJ/m-2·d-1)RMSE(MJ/m-2·d-1)黑龍江漠河095282-10568193812

3.2 基于NDVI和SR估算的FPAR結果

以往研究表明,通過NDVI估算出來的FPAR值高于實測值,通過SR估算出來的FPAR值低于實測值,但是SR估算出來的誤差小于直接由NDVI估算出來的結果,所以基于這種情況,Los SO[22]將這兩種方法結合起來,取兩者的平均值作為FPAR的估算值,此時估算的FPAR與實測值之間的誤差達到最小。本研究最終將基于NDVI估算的FPAR和基于SR估算的FPAR的兩個結果加和取平均值作為FPAR的最終估算值[23]。

3.3 植物吸收光合有效輻射APAR的估算結果

通過ENVI軟件計算研究區遙感圖像的APAR,并利用密度分割工具將APAR的大小分割成13等級,分不同的顏色顯示,運用ArcMap制作成一景遙感圖像中APAR分布圖,具體結果如圖7所示。

從大興安嶺漠河地區APAR分布圖可以看出,APAR值為5~6 MJ/(m2·d)和6~7 MJ/(m2·d)的分布最多,占據了整個區域的65%以上,通過與HSI數據的標準假彩色影像圖進行對比,APAR值為5~6 MJ/(m2·d)所分布的地方為大片平地上的森林,APAR值為4~5 MJ/(m2·d)所分布的地方為坡地上的森林,而APAR值為3~4 MJ/(m2·d)所分布 的地方為平地上的草地和低矮的植被,裸地上沒有APAR值。由于本研究的是大興安嶺漠河地區一天的APAR值,用一天的APAR值估計年APAR值約為1 800 MJ/(m2·a),在對我國植被凈第一生產力的估算中[24-26],東北長白山、小興安嶺、大興安嶺等林區的APAR值一般在1 300~2 200 MJ/(m2·a),所以本研究基本符合要求。

圖7 大興安嶺漠河地區一景遙感圖像的APAR分布圖Fig.7 APAR distribution of remote sensing images of a scene in Daxinganling Mohe area

4 結論

植被吸收光合有效輻射是計算森林生態系統凈初級生產力的一個非常重要的參數,本研究基于HJ-1A星HSI數據的APAR研究方法,以一景遙感圖像作為數據源,運用ENVI對HSI數據進行預處理,根據漠河氣象站的氣象和輻射數據,應用逐日太陽模型模擬出太陽總輻射,利用NDVI和SR估算FPAR,最后基于CASA模型以像元為單位逐一計算出APAR。主要研究結果如下:

(1)本文選取了全局去條紋法對HSI數據旋轉后垂直分布的條紋進行去除,可以在保證信息質量的前提下,最大程度的去除條紋的影響,提高圖像質量。

(2)本文對逐日太陽輻射模擬結果進行統計分析,得出模擬值與實測值相比普遍偏低,但偏差不大,RMSE為1.94 MJ/m2·d,以日照時數和日序數作為模型輸入參數時,如果在時間步長為旬、月、年時,模擬太陽輻射值的工作量增大,模型方法將在之后研究中進一步改進。與已有的太陽輻射模型相比,該模型具有易操作的特點,所需的輸入資料簡單(緯度,日照時數,日序數),實測太陽輻射與模擬太陽輻射有顯著的線性相關關系,模擬值與實測值之間達到了顯著相關水平,決定系數(R2)為0.95,模擬漠河氣象站地區的每日太陽輻射精度高,實現了由漠河氣象站的點輻射值,到一景遙感圖像的面輻射值,可作為APAR計算的重要參數。

對CASA模型的算法進行了改進和簡化,改進了其中的FPAR算法,由于由NDVI計算出的FPAR值偏高,由SR估計的FPAR偏小,其中SR由NDVI計算得到,所以將NDVI和SR的影響結合起來考慮,使得最后得到的FPAR結果更接近于真實值。

本研究獲得的NDVI數據在無植被地區有時會出現較大的值,因此這需要將遙感數據和地面植被數據在計算前后進行比對,盡量消除這種由于數據本身原因而導致結果誤差的特殊情況。

[1]王保林,王晶杰,楊勇,等.植被光合有效輻射吸收分量及最大光能利用率算法的改進[J].草業學報,2013,22(5):220-228.

[2]何明珠.遙感數據與兩葉光能利用率模型同化的總初級生產力模擬[D].南京:南京大學,2013.

[3]Hipps L E,Asrar G,Kanemasu E T.Assessing the interception of photosynthetically active radiation in winter wheat[J].Agricultural Meteorology,1983,28:253-259.

[4]Meere K,McCree J.A solarineter for measuring photosynthetically active radiation[J].Agric Meteorol,1966,3(5-6):353-366.

[5]Chen J M,Liu J,Cihlar J,et al.Daily canopy photosynthesis model through temporal and spatial scalingforremote sensing applications[J].Ecological Modelling,1999,124(2-3):99-119.

[6]張喜旺,于寧,秦奮,等.基于遙感的植被吸收光合有效輻射估算[J].生態科學,2013,32(5):604-608.

[7]李剛,辛曉平,王道龍,等.基于MODIS數據的草地植被光合過程參數估算——以內蒙古自治區為例[J].中國草地學報,2008,33(3):1-7.

[8]陳新芳,安樹青,陳鏡明,等.森林生態系統生物物理參數遙感反演研究進展[J].生態學雜志,2005,24(9):1074-1079.

[9]鈕立明,蒙繼華,吳炳方,等.HJ-1A星HSI數據2級產品處理流程研究[J].國土資源遙感,2011(1):77-82.

[10]高海亮,顧行發,余濤,等.Hyperion遙感影像噪聲去除方法研究[J].遙感信息,2014,29(3):3-7.

[11]Gao B,Montes M J,Davis C O.Atmospheric correction algorithms for hyperspectral remote sensing data of land and ocean[J].Remote Sensing of Environment,2009,113(1):17-24.

[12]齊月,房世波,周文佐.近50年來中國地面太陽輻射變化及其空間分布[J].生態學報,2014,34(24):7444-7453.

[13]左大康,王懿賢,陳建綏.中國地區太陽總輻射的空間分布特征[J].氣象學報,1963(1):78-96.

[14]Szeicz G.Solar radiation for plant growth[J].Journal of Applied Ecology,1974,11(11):617-636.

[15]董丹,倪健.利用CASA模型模擬西南喀斯特植被凈第一性生產力[J].生態學報,2011,31(7):1855-1866.

[16]董泰鋒,蒙繼華,吳炳方.基于遙感的光合有效輻射吸收比率(FPAR)估算方法綜述[J].生態學報,2012,32(22):7190-7201.

[17]Ruimy A,Saugier B.Methodology for the estimation of terrestrial net primary production from remotely sensed data[J].Journal of Geophysical Research,1994,97:18515-18521.

[18]高艷妮,于貴瑞,張黎,等.中國陸地生態系統凈初級生產力變化特征——基于過程模型和遙感模型的評估結果[J].地理科學進展,2012,31(1):109-117.

[19]Los S O,Justice C O,Tucker C J.A global 1by 1NDVI dataset for climate studies derived from the GIMMS continental NDVI data[J].International Journal of Remote Sensing,1994,15:3493-3518.

[20]朱文泉,潘耀忠,張錦水.中國陸地植被凈初級生產力遙感估算[J].植物生態學報,2007,31(3):413-424.

[21]童成立,張文菊,湯陽,等.逐日太陽輻射的模擬計算[J].中國農業氣象,2005,26(3):165-169.

[22]Los S O.linkages between global vegetation and climate:an analysis based on NOAA advanced very high resolution radiometer data[D].Oregon:National Oceanic and Atmospheric Administration,1998.

[23]Peng D,Zhang B,Liu L.Comparing spatiotemporal patterns in Eurasian FPAR derived from two NDVI-based methods[J].International Journal of Digital Earth,2012,5(4):1-16.

[24]陳正華,麻清源,王建,等.利用CASA模型估算黑河流域凈第一性生產力[J].自然資源學報,2008,23(2):263-273.

[25]呂丹桔,施心陵,吳云鵬.基于Labview的植物紅光吸收動態監測[J].林業機械與木工設備,2015,43(1):35-36+42.

[26]孫存舉,雷建華,唐果山,等.基于SPOT NDVI影像的阿壩州植被變化遙感監測[J].四川林業科技,2014,35(6):56-59.

Vegetation Absorbed Photosynthetically Active Radiation Estimates Based on HJ-1A Satellite HSI Data

Chang Yaxuan1,Yu Ying*,Cui Lvyuan2,Hou Jiyu2,Zeng Dewei2

(College of Forest,Northeast Forestry University,Harbin 150040)

In order to measure the effective absorption of solar radiation by vegetation,HSI data obtained by HJ-1A satellite was used to estimate the absorbed photosynthetic active radiation.The basis was provided for the separate calculation of sunlit and shadedAPAR.The accurate parameter data was also provided for the net primary productivity calculation of forest ecosystem.This study includes three main steps:Firstly,the HSI data containing Mohe weather station were pre-processed,which includes format conversion,stripe-removing and atmospheric correction and so on.Secondly,daily solar radiation model was used to calculate the daily total solar radiation,and then theFPARwas calculated byNDVI.FPARand solar radiation were used to obtain theAPAR.Finally the processed image was obtained and the correlation between daily solar radiation model and regional meteorological station data was verified.The precision of FPAR calculated by combination ofNDVIandSRwas improved.The grayscale image was displayed as a thematic image when theAPARcalculation of vegetation was completed.TheAPARvalues of different regions was obtained from the thematic maps,and the annualAPARwas finally calculated.

HSI data processing;solar radiation model;FPAR;APAR

2016-08-30

大學生創新項目(20150225068)

常雅軒,本科生。研究方向:遙感與地理信息系統。

*通信作者:于穎,博士,講師。研究方向:定量遙感與地理信息系統。E-mail:yuying4458@163.com

常雅軒,于穎,崔綠園,等.基于HJ-1A星HSI數據的植被吸收光合有效輻射APAR估算[J].森林工程,2017,33(2):22-27.

S 771

A

1001-005X(2017)02-0022-06

猜你喜歡
模型研究
一半模型
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品黄色片| 一级毛片免费的| 国产超碰一区二区三区| 久久99热这里只有精品免费看| 国产午夜在线观看视频| 国产成人无码综合亚洲日韩不卡| 国产成人亚洲无吗淙合青草| 日韩av资源在线| 亚洲美女AV免费一区| 国产91小视频| 中文无码精品A∨在线观看不卡 | 综合色在线| 国产成人午夜福利免费无码r| 亚洲欧洲自拍拍偷午夜色| 亚洲人成网站在线播放2019| 日本人真淫视频一区二区三区| 欧美色视频在线| 国产成人精品综合| 在线中文字幕日韩| 黄色免费在线网址| 99爱视频精品免视看| 欧美成人在线免费| 黄色污网站在线观看| 国产免费久久精品44| 成人中文字幕在线| 无码又爽又刺激的高潮视频| 99草精品视频| 国产欧美中文字幕| 久久久久久久久久国产精品| 欧美成人亚洲综合精品欧美激情| 在线不卡免费视频| 亚洲国产成人超福利久久精品| a级毛片一区二区免费视频| 亚洲区一区| 亚洲第一极品精品无码| 国产日韩丝袜一二三区| 亚洲午夜福利精品无码| 人妻无码一区二区视频| 色综合中文综合网| 久久青草热| 日韩欧美网址| 国产性猛交XXXX免费看| 亚洲伦理一区二区| 国产欧美日韩精品综合在线| 好久久免费视频高清| 99国产精品国产高清一区二区| 婷婷综合在线观看丁香| 久久综合亚洲鲁鲁九月天| 波多野结衣无码中文字幕在线观看一区二区 | 国产精品一老牛影视频| 国产69精品久久久久妇女| 国产成人区在线观看视频| 国国产a国产片免费麻豆| 国产成人高清精品免费软件| 91系列在线观看| www.99精品视频在线播放| 国产精品亚洲五月天高清| 91精品国产福利| 国产日本欧美亚洲精品视| V一区无码内射国产| 91小视频在线观看免费版高清| 亚洲91精品视频| 国产人妖视频一区在线观看| 色婷婷电影网| 91久久国产成人免费观看| 欧美一区中文字幕| 久久青草精品一区二区三区| 中字无码精油按摩中出视频| 国产成人无码播放| 欧美精品1区| 国产精品自在拍首页视频8| 蜜臀AV在线播放| 日韩在线视频网站| 亚洲AV无码乱码在线观看裸奔 | 伊人丁香五月天久久综合 | 亚洲精品午夜无码电影网| 国产丰满大乳无码免费播放| 免费可以看的无遮挡av无码| 91无码人妻精品一区| 国产成人久久综合一区| 黄色三级网站免费| 久久综合五月婷婷|