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

基于地球河流扇預測模型的火星地表水徑流量預測

2024-02-26 10:02:12張元福張森黃云英孫世坦袁曉冬王敏張曉晗陳冬
沉積學報 2024年1期

張元福,張森,黃云英,孫世坦,袁曉冬,王敏,張曉晗,陳冬

1.中國地質大學(北京)能源學院,北京 100083

2.中國地質大學(北京)科學研究院,北京 100083

3.中國石油勘探開發研究院,北京 100083

4.中國石油化工股份有限公司石油勘探開發研究院,北京 100083

0 引言

自20世紀60年代起,人類對火星的探測從未間斷,并在火星沉積學、地形地貌、物質成分、水環境等方面有所進展[1],隨著天問一號在烏托邦平原的成功著陸,我國的火星探測也進入了新的歷史階段。

水活動歷史和沉積體系演化一直是火星研究的中心主題,而且二者通常相互關聯。許多研究表明火星上具有豐富的沉積歷史、水環境歷史、古生命跡象等[2-5]。在火星的后期演化中,其水分和大氣逐漸逸失,導致火星氣候干燥、溫度下降、大氣稀薄[6],氧化作用、水活動和化學風化隨之減弱,最終使火星早期沉積演化的產物更容易保存下來[7]。這為通過沉積學研究推測火星水活動提供了可能[8-9]。

目前,大多數研究者通過遙感探測、光譜成像、衛星圖像等技術對火星水活動歷史和沉積體系進行研究[10-14]。來自火星隕石的同位素證據和來自大氣測繪數據的解釋表明,火星在Noachian 早期(大約在4 100 Ma 之前)可能存在一個全球地表水水庫[15]。McCubbinet al.[16]和Hurowitzet al.[17]對火星上的含水巖漿作用探究發現,含水巖漿活動可能是古代地表水的重要來源。有大量的地貌學和礦物學證據表明早期火星上存在液態流體,但其性質可能與地球上的水不同。Fukushiet al.[18]通過對Gale隕石坑內環形山湖泊沉積物進行分析,證實了早期Gale 湖水的高鹽度性,并認為火星表面流體偏酸性,碳酸鹽礦物的選擇性溶解廣泛存在。同時,沉積體系研究也是水活動歷史研究的一部分[19-20]。許多研究者推測火星上大量溝壑地形,特別是火星南部高地的數百個溝谷網絡,是由液態水所沖蝕而成[21-25]。而火星上不斷被發現的類似地球的三角洲、沖積扇和河道沉積也是地表水曾經存在的重要依據[26]。美國航空航天局的“機遇號”曾在Meridiani 平原發現了層狀沉積巖,Squyreset al.[14]的研究表明這種層狀沉積巖是由水侵蝕產生的硅酸鹽、硫酸鹽、赤鐵礦等物質構成。Morganet al.[12]在研究火星Saheki 隕石坑沖積扇時發現,Saheki扇中的大部分層狀沉積物由漫灘沉積和淺水湖泊沉積的細沉積物組成。

綜上,眾多學者通過對火星上沉積物和沉積巖的研究,均推測火星在地質歷史中存在較為活躍的水環境,并且改變了火星上的大量巖石[11]。但是,對火星上地表水的規模、流量、降水量等定量參數缺乏有效的研究手段。

2016 年,Birchet al.[27]通過研究土星最大的衛星——土衛六上沖積扇和河流扇的形成與分布,推斷了河道徑流量對扇體演變的影響,并且將扇形沉積體的沉積特征與地球上相似的沉積體進行了對比。與地球相比,火星的質量、重力加速度較小(火星上的重力加速度為3.71 m/s2,僅為地球上的37.9%)。火星與地球在重力加速度方面存在差異,而重力對水流作用和沉積物運移均有影響[28]。在低重力環境下,火星的水流作用變弱,顆粒運動和懸浮所需的剪切應力值也降低。因此,火星和地球在沉積物運移能力上差異不大[29]。這為通過沉積體系研究火星地表水活動提供了一個思路,也可以進一步論證地外行星存在水的證據以及恢復其古沉積環境等[21]。河流扇的概念是在沖積扇的研究中逐漸演化和分離出來的。Schumm[30]于1977年在The Fluvial System一書中將河流作用引入具有穩定流體供給的扇形沉積體中,是河流扇概念的基本雛形。2020年,河流扇被定義為一種發育在山口或平原地帶,內部以河流沉積作用為主,在地質時間尺度上形成沿上游頂點向下游區域發散的扇狀沉積體[31]。相較于觀測火星上的沖積扇,河流扇具有更大的面積和更穩定的河道供給。這些特征使其比沖積扇更容易觀測。而且,沉積過程的相對穩定也使得通過河流扇研究水活動歷史成為了可能。

通過建立地球河流扇數據集和扇體預測模型,對火星河流扇進行識別和分析,進而定量推測火星局部區域的地表水徑流量及其來源。該研究將為進一步深入研究火星地表水活動提供一個新的方法。

1 地球河流扇數據集的建立

1.1 地球河流扇數據采集

現代河流扇沉積的識別和數據采集主要利用美國谷歌公司的衛星圖片。通過對世界范圍內的河流扇進行識別和全方位測量,建立包含河流扇形成和發育信息的數據集。

對全球范圍的河流扇識別,主要分三個步驟:首先,對全球河流進行檢索,標定山區季節性河流和常年穩定性河流,篩選出河流扇最有可能發育的區域。其次,在篩選出的區域進行河流扇的初次判別。判別標準如下:(1)上游方向有穩定的水道供給,在某一點(頂點)向下游方向呈放射狀展布,形成扇形沉積體;(2)海拔高度自頂點向河道下游方向下降;(3)從頂點往下游方向沒有其他支流匯入;(4)終止于陸地(鹽湖、濕地、末端干涸),終止于軸向河流,或終止于穩定水體(湖泊、海洋)。最后,需要對識別出的扇體進行復核確認,并根據河流扇的發育背景和沉積特點,把沖積扇、三角洲等相類似的沉積體排除[31]。

采集的數據主要包含平面發育要素數據和控制因素數據兩大類。

平面發育要素數據主要包括河流扇長度、寬度、面積、內部河道形態及末端終止類型。長度是指扇體的頂端至扇端邊緣的距離;寬度是二分之一長度處的扇體寬度中值;面積包含活動的扇體面積和廢棄部分的面積;扇體類型是根據扇體是否發育受限劃分;河道樣式是指扇體內部的河道形態;終止類型主要是終止于陸地、軸向河流或穩定水體。

控制因素數據是指不依賴于扇體本身的各類環境參數。我們錄入了7 種控制因素數據。盆地類型數據是通過美國地質調查局發布的全球盆地分布圖確定;氣候類型數據和降雨量數據通過美國氣象局公布的全球氣候分布帶和全球降雨量分布確定;經緯度、坡度、集水區面積和扇體至山前距離數據則是通過谷歌地圖內的定位、測距、測高和測量面積這四個功能采集,數據類型及單位見表1。

表1 河流扇數據類型采集表Table 1 Fluvial fan types

1.2 地球河流扇的分布和特征

通過采集到的數據,發現全球河流扇分布具有明顯的緯度分帶特征。大多數扇體發育在中緯度地區。在北半球,河流扇集中在30°~50° N 緯度帶(49.5%),在南半球集中在15°~35° S 緯度帶(21.9%),這兩個緯度帶發育了全球(71.4%)的河流扇(圖1)。

圖1 河流扇全球分布情況條形圖為各緯度河流扇的發現數量,紅框為4個主要發育區Fig.1 Global distribution of fluvial fans the bar chart shows the number of fluvial fans found in different latitudes,the red box shows four main development areas

此外,河流扇在不同的氣候、盆地類型下的分布也存在明顯差異。河流扇主要分布在半干旱和干旱氣候條件下,分別占統計樣本總量的68.4% 和14.4%。熱帶季風氣候下也較為常見,為8.4%,其他氣候下發育的河流扇相對較少(圖2a)。如圖2b 所示,河流扇主要發育在前陸盆地內,占統計樣本總量的50.9%,克拉通盆地和大陸裂谷盆地占比分別為23.3%和12.3%。河流扇在被動陸緣盆地、弧前盆地、弧后盆地發育較少,在走滑盆地中河流扇發育最少,僅為1.3%。通過對河流扇形成和發育的控制因素分析,發現干旱和半干旱氣候下的前陸、克拉通和大陸裂谷盆地都有利于河流扇發育。其中,氣候條件具有明顯的規律性,三種有利的盆地構造類型分別屬于擠壓性、拉張性和穩定性的克拉通等盆地。這說明相對于構造條件,氣候是河流扇形成的關鍵條件。

圖2 河流扇在不同氣候和盆地類型下的發育情況(a)不同氣候下的河流扇發育情況;(b)不同盆地類型下的河流扇發育情況Fig.2 Development of fluvial fans in different climates and basin types(a) development of fluvial fans in different climates;(b) development of fluvial fans in different basin types

河流扇形狀大多為自頂點向下游發散的扇體,其中部分扇體受附近地形地貌的影響,向兩側發散受限,最終發育成狹長的帶狀扇體。據此,將河流扇分為限制型河流扇和非限制型河流扇。非限制型扇體是指不受周圍地形明顯限制,扇體形態完整,圓心角大于90°。限制型扇體是指受到高山或峽谷影響,呈條帶狀,圓心角通常小于90°。統計結果表明,限制型河流扇占比為46%,非限制型河流扇占比為54%(圖3a)。限制型扇體長寬比的最大值為6.98,最小值為0.54,平均值為2.07;而非限制型扇體的長寬比的最大值為5.12,最小值為0.37,平均值為1.60(圖3b)。總體上,發育狹長的限制型河流扇,長寬比相對較大。

圖3 河流扇外部形態比例(a)外部形態比例;(b)外部形態與長寬比關系Fig.3 Proportion of fluvial fan shapes(a) external morphological proportion;(b) relationship between external morphology and aspect ratio

2 地球河流扇發育面積預測模型

2.1 河流扇形成和發育面積的控制因素

基于數據集,我們發現構造、氣候是控制河流扇發育的兩大環境因素。在此基礎上,分析了影響河流扇發育面積可量化的直接因素,包括緯度、區域年平均徑流量、集水區面積、距山前距離等。

河流扇面積主要集中在100~10 000 km2,其中100~1 000 km2占 比41.4%,1 000~10 000 km2占 比35.1%(圖4a)。通過對河流扇發育面積與各類控制因素作相關分析可以發現,河流扇發育面積與緯度、地形坡度、集水區面積、距山前距離、年平均徑流量都有正相關關系。

圖4 河流扇面積分布及控制因素(a)河流扇面積分布;(b)面積與緯度(北緯);(c)面積與年徑流量;(d)面積與坡度;(e)面積與集水區面積;(f)面積與山前距離Fig.4 Distribution of fluvial fan area and controlling factors(a) fluvial fan area distribution;(b) area and latitude (north latitude);(c) area and annual runoff;(d) area and slope;(e) area and catchment area;(f) area and piedmont distance

首先,河流扇的面積與所在緯度有一定關系,在北半球,扇體的面積隨著緯度的增加而減小(圖4b)。緯度與氣溫直接相關,這反映了溫度對扇體面積的強烈作用。在中低緯度地區,溫度較高,且由于地勢平坦,河流扇內部河道的流向受地形的限制小,更容易發散,從而形成面積更大的河流扇體。由于南半球陸地面積小,河流扇發育主要集中在澳大利亞的干旱地區(圖1),扇體的面積和緯度之間并無明顯規律。

其次,河流扇的面積與坡度呈明顯的負相關關系,扇體的面積隨著坡度的增加而減小(圖4c)。面積大于1 000 km2的扇體坡度普遍小于1°。面積大于10 000 km2的扇體坡度普遍小于0.1°。

扇體面積和山前距離也具有很好的線性關系(圖4d),河流扇距離山口越遠,形成的河流扇面積越大,反映了地形地貌對扇體面積的強烈控制作用。河流扇距離山口遠,坡度降低,同時也表明河流扇上游河道的輸移能力較強,更容易形成面積較大的河流扇。

最后,河流扇面積和集水區面積呈明顯的正相關(圖4e),相關系數達0.445。集水區面積越大,形成的河流扇面積越大。同時,河流扇面積與集水區年徑流量也具有一定的線性關系(圖4f),反映了物源供給對扇體面積的控制作用。集水區年徑流量大的河流扇,其物源供給更加充分,更易形成面積更大的河流扇。

2.2 河流扇發育面積的神經網絡預測模型

為了保證預測模型的有效性,預測模型的參數應與河流扇體發育的控制因素有關,且易于獲取。基于控制因素分析,氣候、構造、物源條件是河流扇形成的主要控制因素。從河流扇數據集分析結果來看,扇體面積與緯度、降雨量、坡度、集水區面積和山前距離關系最密切(圖5a)。研究采用神經網絡建立扇體的非線性預測模型。通過對345 個河流扇樣本的神經網絡訓練,建立了神經網絡預測模型,并對5%的獨立樣本進行檢驗。發現河流扇的面積與年徑流量、緯度、距離山口距離、地形坡度及集水區面積具有很好的相關性。訓練樣本的擬合相關性達0.889(圖5b),測試樣本的擬合相關性達0.912(圖5c)。

圖5 地球上河流扇面積的神經網絡預測(a)河流扇面積神經網絡訓練示意圖;(b)訓練樣本預測結果與實測結果對比;(c)測試樣本預測結果與實測結果對比;(d)非限制型訓練樣本預測結果與實測結果對比;(e)非限制型測試樣本預測結果與實測結果對比Fig.5 Neural network prediction of the areas of fluvial fans on Earth(a) schematic diagram of fluvial area neural network training;(b) comparison between training sample prediction results and measured results;(c) comparison between the predicted results of test samples and the measured results;(d) comparison between the predicted results of unrestricted training samples and the measured results;(e) comparison between the predicted results of unrestricted test samples and the measured results

由于河流扇可以根據外部形態分為限制型和非限制型兩種類型。這兩種類型由于受控制因素的影響不同,其預測模型也必然不同。相比較而言,非限制型扇體由于扇體形態的發育受限制比較小,和控制因素之間的相關性更高。基于上述分析,通過對186個非限制型河流扇樣本進行神經網絡訓練,并隨機抽取10個獨立樣本進行檢驗。發現非限制型河流扇的面積與年徑流量、緯度、距離山口距離、地形坡度及集水區面積具有更好的相關性。訓練樣本的擬合相關性達0.998(圖5d),測試樣本的擬合相關性達0.949(圖5e)。

3 火星河流扇的識別和及對比

3.1 火星河流扇的識別

在火星南部高地廣泛分布著形成于諾亞紀和西方紀時期的網狀河谷和沖溝[25,32-33]。網狀河谷單個寬度一般小于5 km,長度可達數百至數千千米。其中,大量的網狀河谷屬于沖積扇,也有部分網狀河谷從特征上判斷屬于河流扇。隨著火星表面遙感影像的分辨率提高,部分河流扇可以被準確地識別出來。如圖6a是位于火星南部霍爾頓(Holden)隕石坑(26°S,34° W)東北方向一個河流扇發育點(24.3° S,33.5° W),衛星照片(1.4~6.0 米/像素)由美國航空航天局的火星軌道飛行器攝像機(MOC)拍攝。以下簡稱為霍爾頓河流扇。Malinet al.[34]把這個網狀沉積體籠統地解釋為扇形沉積體,并以此推測了火星上地表水的持續流動。根據最新的沉積學研究[31-32],該扇體上游具有穩定的河道供給,自頂點向下游發散,更符合河流扇沉積的特點。

圖6 火星霍爾頓河流扇識別示意圖(a)火星霍爾頓河流扇的位置及集水區范圍;(b)地球安德胡伊河流扇的扇體形態和期次;(c)火星霍爾頓河流扇的扇體形態和期次Fig.6 Schematic diagram of fluvial fan in the Holden crater on Mars(a) location and catchment area of Holden fluvial fan in Mars;(b) fan shape and stage of Andkhoy fluvial fan on Earth;(c) fan shape and stage of Holden fluvial fan on Mars

根據地球河流扇數據集,阿富汗的安德胡伊(Andkhoy)河流扇(65.09° E,35.88° N)與火星的霍爾頓河流扇非常類似。安德胡伊河流扇(圖6b)屬于非限制型河流扇,面積為491.35 km2,長31.22 km,寬25.5 km,坡度0.101°,其山前距離為43.2 km,集水區面積為6 316.8 km2。而霍爾頓河流扇(圖6c)面積為84.35 km2,長10.19 km,寬13.43 km,坡度為1.174°,距山前距離2.68 km。該河流扇上游存在兩級分水嶺,分別對應兩級集水區,其中二級集水區面積為698.4 km2,一級集水區面積為5 262.3 km2。

地球和火星上的兩個河流扇形態類似,扇體發育受側向限制較小,形態發育完整,各自發育了三期扇體,屬于非限制型河流扇。唯一顯著的不同在于霍爾頓河流扇存在倒轉地形。即古河道在水流作用下發生膠結、壓實等成巖作用,從而形成抗侵蝕表面。在環境變得干旱且剝蝕作用強烈時,該侵蝕面使河道沉積不易被侵蝕,形成了脊狀凸起[35-36]。

3.2 火星河流扇的發育參數分析

火星霍爾頓河流扇按照地球河流扇的分類,屬于非限制型河流扇。基于地球河流扇數據集,緯度、坡度、年徑流量、集水區面積、距山前距離被認為直接影響著河流扇的發育面積。上述提到的五個控制因素中,除了年徑流量,其他四個參數都可以通過衛星圖片進行測量得到。

通過將火星霍爾頓扇體發育面積與四個參數的相關性投影到地球非限制型扇體的參數關系圖上,大部分參數點都位于正常區間。這表明火星霍爾頓河流扇各參數與面積的關系與地球河流扇大致相同。其中緯度(圖7a)、坡度(圖7c)、山前距離(圖7d)與面積的關系與地球河流扇幾乎一致,反映了其控制機理與地球河流扇完全相同。二級集水區面積與扇體面積的關系與地球河流扇也幾乎相同,但一級集水區面積與扇體面積關系顯然超出了地球河流扇的范圍(圖7b)。因此,可以推測火星霍爾頓河流扇主要受二級集水區控制。也就是說,在一級集水區內部存在沉積物卸載區,這與衛星圖的顯示相吻合。值得注意的是,與地球相比,雖然緯度對面積控制未超出地球河流扇的范圍,但相同緯度上火星霍爾頓河流扇的面積更小。

圖7 火星霍爾頓河流扇參數關系與地球非限制型河流扇參數關系對比(a)非限制型河流扇面積與緯度關系;(b)非限制型河流扇面積與集水區面積;(c)非限制型河流扇面積與坡度;(d)非限制型河流扇面積與山前距離。紅色點為火星,藍色點為地球Fig.7 Interrelationship of parameters for the Holden crater fluvial fan on Mars and for unrestricted fluvial fans(a) relationship between the area of unrestricted fluvial fan and latitude;(b) area of unrestricted fluvial fan and catchment area;(c) area and slope of unrestricted fluvial fan;(d) area of unrestricted fluvial fan and piedmont distance

4 火星河流扇發育期的地表水徑流量預測

4.1 徑流量預測

在西方紀晚期的火星峽谷中可能出現了湖泊等大規模水體,且持續地流向鄰近洼地,與外流通道相融合[37-38]。水的出現可能是西方紀氣候變暖導致的降雨結果,或者是行星傾角改變時期氣候的改變導致冰雪融化的結果[39-40]。火星地表在經歷了后期水流和風等外界環境的不斷侵蝕和改造,逐漸形成了現有的復雜地貌[41-43]。同時,地表徑流會促使沉積物在地表的輸移和沉積,從而形成沖積扇、河流扇、河流、三角洲等沉積體系[44-46]。這些沉積體系都受地表水徑流的影響,均可以用來分析地表徑流量。要實現對徑流量的定量分析,需要滿足兩個條件,一是流域和集水區要相對集中,便于分析和計算。二是要盡量選擇水活動穩定,受突發的事件性沉積影響小的穩定沉積體。沖積扇和河流扇的流域和集水區相對集中,更容易用來預測局部的徑流量。同時,在火星的南部高地確實發現了大量類似沖積扇和河流扇的網狀沉積體系。相比較而言,沖積扇受突發性洪水和重力流沉積影響大,而河流扇相對比較穩定,受地表徑流作用時間長,可以更好地反映地表徑流的影響。

采用上文建立的地球非限制型河流扇面積的神經網絡預測模型,將火星上霍爾頓河流扇的面積及其他要素輸入,即可反向求取火星河流扇的區域年徑流量。霍爾頓河流扇的面積、坡度、緯度、山前距離等參數可以比較準確地獲取。但集水區面積存在一定的不確定性,從圖6a可以看出,霍爾頓河流扇的集水區存在著兩級分水嶺。如果在一級分水嶺控制的集水區內部不存在沉積物卸載區,如湖泊或小型扇體,則一級分水嶺對應的集水區為霍爾頓河流扇的集水區。反之,則應采用二級分水嶺控制的集水區面積。雖然現有的圖像分辨率難以確定一級分水嶺的集水區內部是否存在扇體或湖泊,但通過前文參數分析,傾向于二級分水嶺和二級集水區控制。最終,求得形成霍爾頓河流扇所需的地表水年徑流量為1.882 8×106m3。

如果徑流量完全為降水構成,則該地區的年平均降水量可達2 695.8 mm。而地球只有赤道地區年平均降雨量才能達到2 000 mm。而火星霍爾頓河流扇處于南緯23°,所以,參考地球情況,火星上該地區降水難以達到如此豐度。參考諾亞紀和西方紀的氣候變遷,推測該地表徑流應該由雪山(冰川)融水和降水兩部分構成。

據此,推測在火星的諾亞紀晚期至西方紀早期,由于氣候的短暫變暖,產生了火星大氣水和地表水循環[33]。在南部高地許多山前地帶,降水和雪山(冰川)融水產生的地表徑流攜帶著侵蝕的沉積物由局部高地進入了鄰近洼地。由于地勢的突然平緩,河道由頂點向下游發散,大量沉積物開始卸載,促使了火星河流扇沉積的形成(圖8)。

圖8 火星南部中緯度地區河流扇沉積示意圖[35]Fig.8 Schematic diagram of fluvial fan deposition in the mid-latitude area of southern Mars[35]

4.2 研究的局限性和展望

火星表面的扇形沉積體一般被認為是火星表面液態水對地表進行侵蝕與改造而形成,但也有人認為是冰川作用等其他作用形成,對于液態水在火星表面是長期還是短期穩定存在的問題也具有爭議[21,26-27]。通過與地球上的河流扇作對比,推斷火星上的霍爾頓河流扇是地表液態水徑流持續作用的產物。由于火星與地球在壓力、溫度、重力等方面的差異,火星上液態水存在的狀態可能與地球存在差異,必然使對比研究存在不可避免的誤差。這種誤差只有對火星液態水存在的性質和狀態進行深入研究后才能消除。同時,推測了形成霍爾頓河流扇的地表水徑流量,并指出其地表徑流不可能由單一降水構成。其液態水可能為降雨、降雪、雪山(冰川)融水、地下水等多種來源。但是,也不能排除火星出現極端氣候條件的情況,這些只能在采集火星土壤和巖石樣品并分析研究后才能得到準確答案。

5 結論

基于地球河流扇沉積的數據集,對火星上的霍爾頓河流扇進行了識別和分析,并對其形成區域的地表水徑流量及其來源進行了預測,取得了如下認識。

(1)地球河流扇在世界范圍內不均勻分布,呈現地域分區,緯度分帶特征。亞洲中部及中國西部是河流扇的主要發育區。30° N~50° N 是河流扇的主要發育帶。

(2)地球河流扇的外部形態主要受控于氣候以及局部地形地貌。影響河流扇發育面積的直接控制因素有緯度、坡度、集水區面積、區域年徑流量、山前距離等。

(3)通過對比分析,火星上的霍爾頓河流扇屬于非限制河流扇,預測扇體形成所需的地表水徑流量為1.882 8×106m3。

(4)火星霍爾頓地區的大氣降水難以提供該河流扇形成所需的足夠徑流量,由此推測該區域地表水為多種來源,如雪山(冰川)融水和大氣降水混合來源。

(5)火星上河流扇與地球河流扇對比分析將推動行星際沉積學研究,并為火星水活動歷史研究提供新的方法和思路。

致謝 該研究是在中國地質大學(北京)科學研究院開展完成的。研究中所用的數據可從通信作者處獲得。由于保密方面的限制,部分底層數據無法公開。

主站蜘蛛池模板: 国产美女在线观看| 十八禁美女裸体网站| 国产99热| 国内精品自在自线视频香蕉| 99精品高清在线播放| 久久久久人妻精品一区三寸蜜桃| 国产欧美日韩专区发布| 欧美第九页| 欧美人与性动交a欧美精品| 91青青草视频在线观看的| 国产麻豆aⅴ精品无码| 久久精品人人做人人综合试看| 伊人久久大香线蕉aⅴ色| 成年A级毛片| 精品综合久久久久久97超人| 亚洲无码高清一区二区| 狠狠ⅴ日韩v欧美v天堂| 国产青榴视频| 欧美成人怡春院在线激情| 久久久噜噜噜久久中文字幕色伊伊| 美女无遮挡免费视频网站| 区国产精品搜索视频| 亚洲视频四区| 国产精品亚洲一区二区三区z | 国产精品浪潮Av| 欧美在线伊人| 97se亚洲综合| 在线观看的黄网| 欧美综合成人| 成人无码区免费视频网站蜜臀| 国产情精品嫩草影院88av| 狠狠干综合| 综合色88| 99视频在线精品免费观看6| 午夜国产不卡在线观看视频| 97青青青国产在线播放| 亚洲美女AV免费一区| 中文字幕在线观看日本| 亚洲自拍另类| 亚洲精品动漫在线观看| 国产高清色视频免费看的网址| 第九色区aⅴ天堂久久香| 国产精品蜜臀| 奇米影视狠狠精品7777| h视频在线播放| 国产男女免费视频| 日本国产一区在线观看| 亚洲欧洲自拍拍偷午夜色| 国产自在线播放| 国产丝袜第一页| 国产特一级毛片| 免费AV在线播放观看18禁强制| 美女裸体18禁网站| 日韩成人午夜| 婷婷综合色| 亚洲免费福利视频| 国内精品一区二区在线观看| 精品无码一区二区在线观看| 日韩福利视频导航| 国内精品自在欧美一区| 国产白浆在线| 日本少妇又色又爽又高潮| 国产精品 欧美激情 在线播放| 精品撒尿视频一区二区三区| 久久久精品久久久久三级| 成人免费视频一区二区三区 | 91啪在线| 99伊人精品| 91精品啪在线观看国产| 特级aaaaaaaaa毛片免费视频| 欧美色香蕉| 国产精品成人第一区| 亚洲人成影院在线观看| 亚洲资源站av无码网址| 精品黑人一区二区三区| 久久精品无码专区免费| 中字无码精油按摩中出视频| 国产无遮挡猛进猛出免费软件| 国产福利免费在线观看| 在线国产91| 久久精品国产一区二区小说| 亚洲不卡影院|