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

天然針闊混交林下穿透降雨特征及其模擬

2014-01-26 08:42:54康永祥劉建軍夏國威劉小林
水土保持通報 2014年3期
關鍵詞:模型

周 威,康永祥,劉建軍,夏國威,劉小林

(1.西北農林科技大學 林學院,陜西 楊凌712100;2.賀蘭山森林生態系統定位研究站,寧夏 銀川750021;3.甘肅小隴山森林生態系統定位研究站,甘肅 天水741022)

水在森林生態系統中的循環與分配整合了能量流動和養分循環等生態過程[1],因而水文功能在森林生態系統的作用中是非常重要的。張焜[2]、王波等[3]、孔維健等[4]、劉小林等[5]分別在不同地域對天然針闊混交林的水文功能進行了研究,為探求天然針闊混交林的水文作用提供了有價值的參考。但由于研究對象的林分特征、氣候條件、研究側重點等存在差異,要全面認識天然針闊混交林的水文效應還任重而道遠。

森林引起生態系統內降雨分配的時空差異,對林地水文過程產生顯著影響[6]。國內外針對林冠層分配降水研究較多,但在秦嶺西段小隴山的研究鮮見報道[7]。小隴山林區是兼有中國南北氣候特點的典型天然次生林區,天然次生針闊混交林是一個重要森林類型。研究小隴山林區森林冠層下穿透降雨規律及與大氣降雨之間的相關性,對揭示區域森林生態水文作用具有重要意義。

在森林生態系統中,不同冠層對降雨的截留效應響應不同,導致了林內穿透降雨的空間分布異質性[8]。地統計學是一種空間分析方法,用于描述區域化變量的空間分布特征,并提供了一種無偏最優的空間插值方法,可用來充分認識空間結構的特征,已廣泛應用于礦產、土壤、生態、氣象等因子分布圖的繪制[9]。本文試圖將地統計學的理論方法應用到穿透降雨的研究中,探討穿透降雨空間分布的特征。

本文通過對小隴山天然次生針闊混交林下穿透降雨特征及其空間分布格局與降雨量關系的分析研究,旨在進一步明確針闊混交天然林生態系統的水分運轉過程及功能的生態學機制,以求更深入地探究森林水文的過程及作用機理,以期為正確認識小隴山林區的森林水文效應及流域森林生態水文功能提供一定的觀測數據和經驗參數,對制定森林管理措施提供科學依據。

1 研究區概況

研究區位于甘肅省小隴山林業科學研究所沙壩實驗基地,地處小隴山林區中心地帶,行政區劃屬天水市秦州區娘娘壩鎮,東接觀音林場,南接高橋林場,西鄰麻沿林場,北連李子園林場。屬秦嶺南坡石質山地,海拔在1 550~2 100m,年平均氣溫1.2~8.8℃,最高氣溫30.3℃,最低氣溫-22.4℃,≥10℃的年積溫2 480.4℃,初霜期10月16日,終霜期5月4日,無霜期154d,年均降雨量757mm,平均年蒸發量1 012.2mm,平均相對濕度78%,屬大陸性季風氣候,為暖溫帶濕潤區。土壤以山地褐土和山地棕壤為主,成土母質為巖石,土層厚度30—55cm,土壤質地為砂壤,表土層腐殖質含量豐富,礦物質養分含量中等,呈微酸性到中性。

樣地設在沙壩實驗基地亂石窖溝,地理坐標為105°53′72″—105°53′73″E,34°8′56″—34°8′57″N,海拔1 630~1 660m,樣地坡度15°~25°,南坡中坡位。樣地為以銳齒櫟(Quercus aliena var.acuteserrata)為主的天然次生針闊混交林,幾乎無人為干擾,植被繁密。樣地內胸徑大于5cm的喬木總計120株,平均胸徑16.8cm,平均樹高19.9m,平均郁閉度0.93,喬木林群落的樹種組成為4銳齒櫟,2華山松(Pinus armandii),1四照花(Cronus japonica var.chinen-sis),1燈臺樹(Bothrocaryum controversum),2軟闊葉樹種。樣地內下木主要有白檀(Symplocos paniculata)、千金榆(Carpinus cordata)、鵝耳櫪(Carpinus turczaninowii)、葛棗獼猴桃(Actinidia polygama)、灰栒子(Cotoneaster acutifolius)等,平均高度2.5m,蓋度為60%。草本主要有菝葜(SmilaxL.)、羊胡子草(EriophorumL.)、懸鉤子(Rubus L.)等,蓋度15%。

2 研究方法

2.1 試驗布設與數據收集

林外大氣降雨量資料來源于小隴山森林生態系統定位研究站,位于甘肅省小隴山林業實驗局林業科學研究所沙壩實驗基地,地理坐標為105°54′E,34°07′N。在沙壩實驗基地天然次生林中選取一塊投影28.5m×27m的樣地,借助測繩和羅盤儀沿等高線精確圈出樣地。在樣地內網格狀機械布設100個雨量桶(直徑20cm),保持雨量筒高度距地面60cm,按品字形水平布設,行與行、列與列之間投影距離均為3m,用于收集樣地的穿透降雨量。每次降雨后及時測定各個雨量桶內的水量,計算雨量時根據收集面積將雨量桶內每次降雨量換算為mm單位。

2.2 數據分析

采用Microsoft Excel 2007和SPSS 18.0對數據進行一般處理分析,利用穿透降雨的測定數據,計算得到每次降雨事件每個雨量筒的穿透率(相同時間下林內穿透降雨量與林外大氣降雨量的比值)。完成基本統計分析后,用SPSS 18.0對數據通過直方圖和正態QQ圖進行正態分布檢驗,確定各穿透率數據均符合正態分布,不需要進行數據轉換。

然后利用地統計軟件GS+9.0對穿透率數據進行半方差函數分析和參數計算,所謂半方差函數,又叫半變異函數,就是兩點間插值的方差的1/2,其表達式為:

式中:γ(h)——間距h的半方差函數值;N(h)——以h為間距的穿透雨觀測點對的數目;Z(xi)——樣點在xi處的穿透率;Z(xi+h)——與xi相隔間距h處樣點的穿透率。下同。半變異函數值隨樣點間距的增加而增大,并在一定的間距(稱為變程,range)達到一個基本穩定的常數(稱為基臺,sill)[10]。

采用GS+9.0中不同類型半變異模型進行擬合,因相比于圓形模型、指數模型、高斯模型和線性模型而言,球面模型擬合效果最好,故本文所用半變異函數模型均為球面模型[11]。

式中:C0——塊金,即間距為0時的半方差函數值,是由實驗誤差和隨機因素引起的變異;C——結構方差,指空間自相關部分引起的變異;(C0+C)——基臺,表示系統內的總變異,一般情況下,基臺值越大表示總的空間異質性程度越高,反之越小;a——變程,指空間自相關距,表示觀測值之間的距離大于該值時是相互獨立的,小于該值時存在空間自相關性[12]。

采用ArcGIS 10.0這一GIS軟件平臺,利用地統計學中常用的數據插值空間變異分析方法——Kriging插值法,進行穿透降雨空間分布的預測模擬分析。

Kriging法是建立在半變異函數理論和結構分析的基礎上,在有限的區域內對區域化變量的取值進行無偏最優估計的一種方法,可對周圍的實測值進行加權以得出未測位置的預測,公式由數據的加權總和組成:

式中:Z(si)——第i個位置處的測量值;λi——第i個位置處的測量值的未知權重;s0——預測位置;N——測量值數量。權重λi取決于測量點、預測位置的距離和預測位置周圍的測量值之間空間關系的擬合模型。

使用ArcGIS 10.0的地統計分析模塊中普通Kriging法的球面模型對穿透率數據進行插值,預測出樣地中300多萬個點的穿透率值,繪制成柵格圖像數據,直觀地模擬樣地的降雨分布。

通過計算得到各降雨事件的穿透率平均值、標準差(S)和變異系數(Cv),作為評價不同降雨量時穿透率的平均狀況和總體不均勻度的指標。

3 結果與分析

3.1 各降雨事件穿透降雨特征

在本實驗觀測期2011年8月至2012年7月共收集了29次降雨事件數據,總降雨量為947.3mm,總的穿透雨量和穿透率分別為742.6mm和0.783 9。通過對這29次降雨事件實測數據的分析,發現當大氣降雨量小于1.9mm時,林冠截留了幾乎全部降水,沒有穿透雨產生,當大氣降雨量大于1.9mm時,樣地中100個樣點平均穿透降雨量與大氣降雨量緊密相關,穿透降雨量隨大氣降雨量的增加而增大,呈顯著的線性正相關關系(圖1a)。擬合的方程如下:

式中:T——穿透降雨量(mm);P——大氣降雨量(mm)。

隨著降雨量的增加,穿透率也呈增加的趨勢。對穿透率與降雨量的關系進行回歸分析,比較后得出對數函數具有較好的擬合性(圖1b),其關系式為:

式中:Tr——穿透率;P——大氣降雨量(mm)。

圖1 穿透雨量、穿透率及穿透雨量變異系數與大氣降雨量的關系

可見,穿透率隨降雨量的變化基本分為3個階段:在小雨雨量級(單場降雨量10mm以下),穿透率隨著降雨量的增加而迅速增大,為快變期;在中雨雨量級(單場降雨量10~25mm),穿透率增加速率明顯變緩,為漸變期;在大雨雨量級(單場降雨量25~60mm)和暴雨雨量級(單場降雨量60mm以上),穿透率逐漸趨于平穩,波動范圍在0.793 9~0.907 5,為穩定期。表1的穿透率平均數也可以印證這一結論。

穿透降雨變異系數(降雨事件的穿透降雨標準差與穿透降雨量的比值)揭示穿透降雨的空間分布不均勻程度,其隨著降雨量的增加而呈下降的趨勢。對穿透降雨變異系數與降雨量的關系進行回歸分析,經過比較得出Mitscherlich模型具有較好的擬合性(圖1c),Mitscherlich模型是Richards模型的一個特例,最初是描述生物生長量的模型。其關系式為:

式中:Tv——穿透降雨變異系數;P——大氣降雨量(mm)。

穿透降雨變異系數隨降雨量的變化也可分為快變期、漸變期和穩定期3個階段,但變化規律與穿透率隨降雨量的變化相反。在小雨雨量級內,穿透降雨變異系數隨著降雨量的增加而迅速減小,在中雨雨量級之后逐漸趨于穩定,波動范圍在0.193 0~0.193 2。

表1 不同雨量級穿透率的不同指標特征

3.2 不同雨量級下穿透降雨分布格局

穿透率觀測值的半方差函數分析(表2)表明:基臺值(C0+C)隨雨量級的增大而減小,說明降雨量增大使得降雨穿透率在空間上趨向于均勻化;各穿透率半方差函數模型空間自相關范圍(變程a)大致介于5.5~6.2m;模型殘差值RSS是選擇模型的主要依據,所有擬合的半變異函數模型殘差值均極小,RSS越小則擬合效果越好,說明各模型擬合效果都很好;決定系數R2表示穿透率半方差模型的解釋效率,即模型擬合精確度的數字表示,可見暴雨的模擬精度最高,小雨的模擬精度最低,模擬精度隨雨量級的增大而升高,這可能與穿透降雨不確定度隨降雨量的增加而減小有關。

表2 不同雨量級下穿透率半變異函數理論模型及其參數

表2中總體指各樣點29次穿透降雨率平均值的集合,小雨、中雨、大雨、暴雨分別指各樣點在不同雨量級下的穿透降雨率平均值的集合。

塊金與基臺的比值C0/(C0+C)表示實驗誤差和隨機因素引起的空間變異占系統總變異的比例。從表2可見,每組穿透率數據C0/(C0+C)值都很小,說明實驗誤差和隨機因素引起的空間變異很小,主要是空間自相關部分引起的變異。但C0/(C0+C)值隨雨量級的增大而增大,說明隨著降雨量的增加,隨機因素引起的空間異質性逐漸增加,由結構性的空間因素引起的空間異質性逐漸下降,即與樹冠結構的相關性逐漸降低。

從空間自相關的角度來看,C0/(C0+C)表示變量的空間相關性程度,如果C0/(C0+C)<25%,說明系統具有強烈的空間相關性,如果25%≤C0/(C0+C)≤75%,說明系統具有中等的空間相關性,如果C0/(C0+C)>75%,則說明系統空間相關性很弱[13]。而表2中每組數據C0/(C0+C)都顯著小于25%,因此均表現出強烈的空間相關性。其中,小雨的空間自相關最顯著。

用Kriging插值法繪制穿透降雨率在各雨量級下的空間分布圖(圖2),由圖2中可直觀地看出穿透降雨在不同雨量級的空間分布格局及其差異,穿透率的空間分布受雨量級的影響呈現出不同的特點。在小雨條件下,穿透率普遍偏小,其分布東北部及個別區域穿透率相對略高,其他各處較為平均,穿透率處于0.75以下(0.6~0.75)的占到樣地總面積的80%以上。中雨條件下的穿透率略有增大,穿透率在整個樣地中的分布比較平滑勻稱,其分布更具有代表性。大雨條件下,穿透率整體較高,但從表2來看,其分布的不均勻性并不很強,變化相對比較平緩,85%以上的區域穿透率在0.8~0.95。相對來說,在暴雨條件下的穿透降雨分布呈現出一些差異,其分布的不均勻性更明顯更突出。但在不同雨量級下,穿透降雨率分布形態也有不同程度的相似,如東北部均存在一個相對高穿透率區域,中部偏西及西南角的同一位置均存在相對低穿透率區域。因為除降雨量大小外,林冠下穿透降雨的空間異質性同時受多個其他因素綜合影響,包括葉面積指數、冠層厚度、離樹干距離、冠層郁閉度、冠形、雨強、風速、風向、坡度等。

圖2 基于100個樣點在不同雨量級下的穿透降雨率平均值預測的穿透降雨率在樣地內分布

3.3 穿透率實測數據與預測數據的比較

由表1可知,實測與Kriging預測的穿透率的平均數差別很小,可以看成是一致的,但是實測數據與預測數據的標準差和變異系數差別很大,預測數據的標準差和變異系數明顯小于實測數據。且根據表1實測數據可知,雨量級從小雨增大時,穿透降雨率的變異系數減小。由此可知,穿透降雨分布的不均勻性隨雨量級的增加而下降。但預測數據不符合上述規律,在預測數據中,大雨的標準差和變異系數最小,其次是小雨以及中雨,暴雨的標準差和變異系數最大,且其遠高于小雨、中雨和大雨,小雨、中雨和大雨的標準差和變異系數差別不大。由此推斷,除大雨外,穿透降雨分布的不均勻性隨雨量級的增加而增大。實測數據與預測數據產生矛盾。

雖然圖3中預測值和實測值頻數分布的基本趨勢是一致的,但預測值的最大值比實測最大值小,而預測值的最小值比實測最小值大,這意味著預測值比實測值的極差縮小。因此,實測數據頻數分布圖相對矮寬,預測數據的頻數分布圖相對高窄,Kriging預測數據的離散度下降,分布更為集中,于是預測值的分布不能完全再現樣本分布,即存在“平滑效應”。李超等[14]對平滑效應的產生原因做了詳盡闡釋,平滑效應的存在會導致預測值的空間連續程度增強,變異程度下降。

盡管Kriging預測值存在條件偏差,但由于在Kriging方程組中對無偏性進行了強制約束,使其仍然是全局最優的。對于平滑效應的認識也要全面、客觀,要針對具體的研究目的而定。

雖然平滑效應導致了以上的問題,但正是Kriging預測值存在平滑效應才使得Kriging法能夠用來繪制穿透率的等值線和趨勢面,對預測結果的離散程度沒有任何約束的插值方法是無法繪制等值線或等值面的。

因此,Kriging預測數據在反映總體分布規律時數據不一定精確,但可以表示相對大小的趨勢,統計分析應以實測數據為主。同時,對插值結果的可靠性可以用Kriging方差來指示,Kriging方差越大說明結果的可靠性越差,否則結果就越可靠[15]。可見大雨條件下的穿透率插值結果最可靠,暴雨條件下的最不可靠。

圖3 不同雨量級穿透率頻數分布

4 討論

本試驗在觀測期內,降水總量為947.3mm,總穿透雨量和穿透率分別為742.6mm和0.7839。穿透降雨量隨大氣降雨量的增加而增大,呈顯著線性正相關關系。很多相似研究從不同尺度、不同林分、不同氣候等角度均認為線形方程能較好地擬合穿透雨量和降雨量之間的關系[16-18]。

穿透率隨著大氣降雨量的增加而逐漸增大,最后趨于穩定值,對數函數可以較好地描述穿透率隨降雨量的變化。這一結果與翟杰等[19]對紫金山麻櫟的研究,鮑文等[20]對岷江上游油松的研究結果相一致,而劉章文等[21]對祁連山灌叢的研究認為指數函數擬合效果更好。實際上劉章文等所采用的是Mitscherlich模型,采用Mitscherlich模型擬合本實驗數據的穿透率與降雨量的方程為:

可見Mitscherlich模型與對數函數擬合精度完全相同(R2=0.439)。將 Mitscherlich模型式(7)與對數函數公式(5)中的Tr用T/P代替,將其分別經過轉換得到穿透雨量與降雨量的關系公式:

以上兩式與線性方程(4)擬合精度完全相同(R2=0.994),表明其與線性方程均可描述穿透雨量與降雨量的關系。因而可知,本實驗采用對數函數與Mitscherlich模型也可較好擬合穿透率與降雨量之間的關系。李振新等[22]認為穿透雨率同林外降雨量之間的關系可以用logistic曲線方程模擬的結論不適用于本研究,模擬精度相對稍差(R2=0.428),但其穿透率存在一個最大值的特點與Mitscherlich模型相似。

穿透降雨變異系數與降雨量呈負相關關系,Mitscherlich模型可以較好地描述穿透降雨變異系數隨降雨量的變化。說明穿透降雨空間分布的不均勻度隨大氣降雨量的增加有下降的趨勢,最終逐漸趨于穩定值。這與戰 偉 慶 等[23]、Rodrigo A 等[24]的 研 究觀點一致,但所用擬合模型不同,他們分別認為指數模型和對數模型較為合適。戰偉慶[23]認為造成這一現象的原因是,當林冠達到飽和持水量時,林外降雨幾乎全部轉化為穿透雨。而我們認為此時林外降雨并未全部轉化為穿透雨,仍然有樹體吸收和蒸發等過程發生,只是林外降雨轉化為穿透雨的比例趨于穩定。

5 結論

穿透降雨率的半變異函數分析表明,天然針闊混交林冠下穿透降雨率的不均勻度隨雨量級的增大而減弱,空間分布趨向于均勻化。且隨著雨量級的增大,隨機因素對穿透率的影響相對增加,而林冠結構等結構性因素對穿透降雨率產生的影響有逐漸變小的趨勢,但結構性因素始終是最主要的影響因素。從穿透率的插值分布圖可以直觀地看到穿透率的空間分布在不同雨量級存在著明顯的差別。

由于Kriging預測數據存在“平滑效應”,使得預測值的分布不能完全體現實際分布,故Kriging預測數據只用來反映總體分布規律,數據不精確,統計分析應以實測數據為主。

天然次生針闊混交林下穿透雨空間分布的差異對林地的土壤水分分布和養分循環利用及生物多樣性必然會產生重要影響,有待于進一步研究。

[1] Anna A,Anselm R.Trace metal fluxes in bulk deposition,through-fall and stem-flow at two evergreen oak stands in NE Spain subject to different exposure to the industrial environment[J].Atmospheric Environment,2004(38):171-180.

[2] 張焜.重慶四面山4種類型天然林水文功能研究[D].北京:北京林業大學,2012.

[3] 王波,張洪江,杜士才,等.三峽庫區天然次生林凋落物森林水文效應研究[J].水土保持通報,2009,29(3):83-87.

[4] 孔維健,周本智,安艷飛,等.天然次生林和人工毛竹林水文生態特征比較[J].水土保持研究,2010,17(1):113-116.

[5] 劉小林,鄭子龍,藺巖雄,等.甘肅小隴山林區主要林分類型土壤水分物理性質研究[J].西北林學院學報,2013,28(1):7-11.

[6] Sun Ge,Mcnulty S G,Amaty A D M,et al.A comparison of the watershed hydrology of coastal forested wetlands and the mountainous uplands in the Southern US[J].Journal of Hydrology,2002,263(1/4):92-104.

[7] 馬熙淵,馬立琦.小隴山林區主要林分對降水的再分配[J].甘肅林業科技,2003,28(2):1-4.

[8] 譚俊磊,馬明國,車濤,等.基于不同郁閉度的青海云杉冠層截留特征研究[J].地球科學進展,2009,24(7):825-833.

[9] 時忠杰,王彥輝,熊偉,等.單株華北落葉松樹冠穿透降雨的空間異質性[J].生態學報,2006,26(9):2877-2886.

[10] 史利江.基于GIS和地統計學的土壤養分空間變異特征研究[D].上海:上海師范大學,2006.

[11] Esr.ArcGIS Resources[OL].[2013-03-10].http:∥resources.arcgis.com/en/help/previous-help/index.html.

[12] 姜城.不同經營體制下土壤養分空間變異規律及管理技術的研究[D].北京:中國農業科學院研究生院,2000.

[13] 王麗霞,段文標,陳立新,等.紅松闊葉混交林林隙大小對土壤水分空間異質性的影響[J].應用生態學報,2013,24(1):17-24.

[14] 李超.土壤水分的空間變異特性分析和農田干旱評估[D].北京:清華大學,2009.

[15] Yamamoto J K.An alternative measure of the reliability of ordinary Kriging estimates[J].Mathematical Geology,2000,32(4):489-509.

[16] 胡珊珊,于靜潔,胡堃,等.華北石質山區油松林對降水再分配過程的影響[J].生態學報,2010,30(7):1751-1757.

[17] 劉世海,余新曉.北京密云水庫庫區水源涵養林冠層水文特征研究[J].林業科學,2005,41(1):194-197.

[18] 薛建輝,郝奇林,吳永波,等.3種亞高山森林群落林冠截留量及穿透雨量與降雨量的關系[J].南京林業大學學報:自然科學版,2008,32(3):9-13.

[19] 翟杰,張志民,胡海波,等.紫金山麻櫟林降水分配格局研究[J].林業科技開發,2011,25(4):63-67.

[20] 鮑文,包維楷,何丙輝,等.岷江上游油松人工林對降水的截留分配效應[J].北京林業大學學報,2004,26(5):10-16.

[21] 劉章文,陳仁升,宋耀選,等.祁連山典型灌叢降雨截留特征[J].生態學報,2012,32(4):1337-1346.

[22] 李振新,鄭華,歐陽志云,等.岷江冷杉針葉林下穿透雨空間分布特征[J].生態學報,2004,24(5):1015-1021.

[23] 戰偉慶,張志強,武軍,等.華北油松人工林冠層穿透雨空間變異性研究[J].中國水土保持科學,2006,4(3):26-30.

[24] Rodrigo A,Avila A.Influence of sampling size in the estimation of mean throughfall in two Mediterranean holm oak forests[J].Journal of Hydrology,2001,243(3/4):216-227.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 97在线碰| 国产永久免费视频m3u8| 午夜福利无码一区二区| 国产精品所毛片视频| 欧美在线黄| 97视频在线精品国自产拍| 欧美激情视频一区| 国产福利在线免费观看| 亚洲第一天堂无码专区| 99精品在线视频观看| 国产无码网站在线观看| 99热亚洲精品6码| 亚洲精品大秀视频| 成人毛片免费在线观看| 日本免费一级视频| 成人日韩精品| 国产91成人| 一级毛片在线直接观看| 亚洲精品少妇熟女| 亚洲无码在线午夜电影| 成人另类稀缺在线观看| 欧美国产日韩另类| 狼友视频国产精品首页| 午夜日韩久久影院| 国产特一级毛片| 狠狠色丁香婷婷| av一区二区三区在线观看| 久久中文字幕2021精品| 日韩欧美国产成人| 国产国模一区二区三区四区| 国产成人永久免费视频| 亚洲成aⅴ人片在线影院八| 视频二区中文无码| 国产伦片中文免费观看| 亚洲码在线中文在线观看| 久久久久亚洲精品无码网站| 白丝美女办公室高潮喷水视频 | 9999在线视频| 视频二区亚洲精品| 亚洲,国产,日韩,综合一区 | 91在线无码精品秘九色APP| 小说区 亚洲 自拍 另类| 欧美一区二区三区香蕉视| 毛片久久久| 亚洲熟女中文字幕男人总站 | 欧美、日韩、国产综合一区| 久久96热在精品国产高清| 真人免费一级毛片一区二区 | 国产激情无码一区二区免费| 国产亚洲视频免费播放| 成人a免费α片在线视频网站| 91香蕉国产亚洲一二三区| 97在线公开视频| www.亚洲国产| 日韩无码视频播放| 午夜高清国产拍精品| 思思99思思久久最新精品| 丁香五月婷婷激情基地| 99偷拍视频精品一区二区| 亚洲中文字幕国产av| аv天堂最新中文在线| 欧美国产日产一区二区| 日韩福利视频导航| 亚洲国产欧美国产综合久久 | 真实国产乱子伦视频| 亚洲欧美日韩高清综合678| 亚洲精品成人7777在线观看| 69av免费视频| 国产在线日本| 2022国产91精品久久久久久| 国内精品九九久久久精品| 在线看AV天堂| 麻豆国产原创视频在线播放| 国产九九精品视频| 亚洲无码视频一区二区三区| 日韩激情成人| 亚洲一区二区在线无码| 欧美怡红院视频一区二区三区| 亚洲第一av网站| 四虎精品国产永久在线观看| 在线观看网站国产| 91毛片网|