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

基于船位監控系統的拖網捕撈努力量提取方法研究

2016-07-11 08:54:47張勝茂唐峰華黃華文
海洋科學 2016年3期

張勝茂,張 衡,唐峰華,樊 偉,黃華文

?

基于船位監控系統的拖網捕撈努力量提取方法研究

張勝茂1,張衡1,唐峰華1,樊偉1,黃華文2

(1.中國水產科學研究院東海水產研究所,農業部東海與遠洋漁業資源開發利用重點實驗室,上海 200090;2.上海普適導航技術有限公司,上海 201112)

摘要:為了基于船位監控系統提取拖網捕撈努力量,通過統計航速獲得 3個峰值,拖網作業在第 2個峰值,即1~2.1 m/s,拖網作業航向差一般在-50°~50°。利用航速、航向差閾值設定,把拖網船狀態劃分為慢速、作業、航行,然后提取出捕撈作業狀態點,1 423艘拖網船共提取到處于捕撈狀態的點318 433個,合計拖網捕撈時間15 921 h,利用反距離加權插值法生成捕撈強度分布變化趨勢圖。捕撈努力量在漁業資源研究中是重要的參考值之一,與傳統的捕撈努力量計算方法相比,該方法具有實時、大范圍、快速、分辨率高的特點,能夠用于輔助漁業資源保護。

關鍵詞:北斗衛星; 船位監控系統; 捕撈努力量; 航向; 航速

[Foundation: Fundamental Research Funds in Central Public-interest Scientific Institution,No.ECSFRI2014T13; Twelfth Five-year Scientific and Technological Support Project,No.2013BAD13B01; Shanghai Municipal Science and Technology Commission,No.12511501200]

捕撈努力量在漁業資源變動研究中是一個重要的參數,即在捕撈中所做的功,傳統計算方法是由投入生產的漁船數量、總噸位、功率,以及作業人數、作業時間、技術與工藝狀況、投網次數等折算獲得[1-4]。隨著近年來漁船船位實時監控系統逐步得到應用推廣,可獲取高時空精度的漁船船位數據,通過船位數據統計分析能從一個新的角度獲得捕撈努力量。我國南海、東海等海區,相繼應用北斗導航衛星構建了漁船監控系統。通過船位監控系統的船位數據能夠計算出累計捕撈時間,把它作為拖網捕撈努力量,這種方法具有近實時、范圍廣、自動、快速等特點,可以獲得高時空分辨率的累計捕撈時間。聯合國糧農組織(FAO)采用每年總的發動機功率和捕撈作業時間(kW·d)表達全球捕撈努力量,在漁船捕撈方式、捕撈漁區、捕撈魚種確定,并且在時段一定的情況下,累計捕撈時間與漁獲量成正相關,Lee等[5]的研究對此也做了一些驗證。國外[6-7]已經用 VMS(Vessel Monitoring System)信息計算捕撈努力量,用于漁業資源評估。北斗衛星導航系統在漁業中的應用研究經過多年的發展[8-9],已經初具規模,據統計到目前為止安裝北斗衛星導航系統終端的漁船已經有 4萬余艘,通過獲取的船位數據可以做到對捕撈努力量的計算[10-11]。本研究利用拖網船的船位數據,通過航速、航向閾值區分出捕撈狀態,把捕撈作業狀態的累計捕撈時間作為捕撈努力量,將某個區域內漁船累計捕撈小時數作為該區域的捕撈努力量,建立中國近

海漁船實時快速的捕撈努力量估算方法。

1 材料與方法

1.1數據來源

北斗漁船船位數據來源于上海的北斗民用分理服務商,數據主要包括漁船的北斗卡號、經緯度、航速、航向、發報時間。船位數據時間分辨率為3 min,空間分辨率約為10 m。數據庫選用SQL Server 2008,它具有管理分析空間數據的功能,可以確定幾何圖形實例之間的空間關系。

捕撈類型、船名等信息主要來源于“中央財政國內海洋漁船油價補助公示”資料,經過與北斗數據的匹配,確定了3 333艘漁船的類型,其中數量較多的是拖網和刺網。拖網船有 2 212艘,占總量的 66%,本研究以拖網為例進行研究。

1.2方法

漁船是否處于捕撈的狀態通過式(1)判斷,當航速和航向處于某個閾值范圍之內時,處于捕撈狀態。其中,V1和V2是捕撈狀態的航速閾值,D1和D2是捕撈狀態的航向差閾值,兩個閾值通過拖網漁船狀態分析與航速統計來確定。

一個漁區格網內可能有多艘漁船,一艘拖網漁船捕撈分為多個網次,一般每個網次持續幾個小時,一個網次結束后間隔一段時間,然后是下一網次。每個網次又有離散的多個船位點組成,因此某個格網內的累計捕撈時間用式(2)計算。其中,Zi是第i個漁區格網的累計時間,Pi,j,k、Pi,j,k-1是某漁船相鄰的兩個船位點的時間,兩者的差是時間長度。第一次求和是一個網次內累計捕撈時間,第二次求和是一艘拖網漁船一段時間多個網次的累計捕撈時間,第三次求和是第i個漁區格網內所有拖網漁船的累計捕撈時間。

捕撈對象是游動的魚類,當在某個漁區捕撈時,會對周邊的漁業資源量產生影響,某漁區的累計捕撈時間越大,周邊漁業資源下降越快,幾個相鄰點的累計捕撈時間在周邊范圍產生一個資源量影響的趨勢面。通過插值生成趨勢面專題圖,能夠輔助資源量全局變化趨勢分析,本文采用反距離加權插值法[12]生成變化趨勢專題圖。

假定每個待插值點都受到局部影響,而這種影響會隨著距離的增大而減小。插值是以權重進行加權平均,所以該平均值不會大于最大輸入值或小于最小輸入值[13]。如果采樣對于正在嘗試模擬的已有變量來說足夠密集,則基于反距離權重會獲得最佳結果[14-15]。漁區格網中心離散點分布較均勻,在制圖分辨率要求較低的情況下,其密集程度可以滿足局部表面的變化分析[13],公式如下:其中,Zo為O 點的估計值; n 為在估算中用到的控制點數目; λi為預測計算過程中使用的已知點的權重,該值隨著樣點與預測點之間距離的增加而減??; Zi為控制點i 的Z 值。確定權重的計算公式為:

其中,r 為指定的指數,本文反距離加權插值法中的r默認值設為2; di為控制點i 與點O 間的距離。

采樣點隨著與預測值點之間距離的增大,其對預測點影響的權重按指數規律減小。在預測過程中,各采樣點對預測點作用的權重大小是成比例的,這些權重的總和為1。式(4)代入式(3),經過轉換可以得到式(6),文中插值所用公式為式(6)。

2 研究結果

2.1拖網漁船狀態分析

拖網屬于過濾性的運動漁具,依靠漁船的運動拖曳網具,在海底或海水中前進,迫使漁具經過水域中的魚蝦蟹等捕撈對象進入網囊,達到捕撈的目的,捕撈作業過程有放網、拖曳網具、起網3個階段。2012年10月10日浙江奉化北斗卡號為300585的拖網船(以下稱“拖網船300585”),在00: 00~24: 00的軌跡如圖1 a~h點,通過調研以及數據分析可知: 該漁船在00: 01從a點開始是處于航行狀態,02: 30航行到b點(30°33′N ,123°41′E),并在b點附近拋錨休息,05: 07開始拖網捕撈作業,經過b~h點共拖網作業6次,在21: 09到h點(30°21′N ,123°43′E),并在h點附近拋錨休息。

圖1 2012年10月10日拖網船300585軌跡Fig.1 Tracks of trawlers 300585 on October 10,2012

表1是對圖1中2012年10月10日拖網船300585 從00: 00~24: 00進行的船位統計,各節點是漁船狀態變化的點,各分段是漁船所處的狀態,平均速度是各段內的漁船速平均值,累計行程是各段中,點連線的曲線長度,相對距離是兩個起止節點的最短距離,長度比是累計行程與相對距離的比值。A段平均速度較大,長度比最?。?B段和I段船員休息處于拋錨狀態,航速因海流與風的影響而產生,B段2.5 h移動了150 m,I段約3 h移動了72 m; C~H段平均速度在1.3~1.8 m/s,長度比在1.3左右,F段由于漁船彎曲航行捕撈,因此比值較大。

表1 2012年10月10日拖網船300585狀態Tab.1 Status of trawlers 300585 on October 10,2012

圖2是2012年10月10日00: 00~24: 00拖網船300585航向變化,圖中的 A~I段按照狀態劃分,其分段與圖1和表1中相同,漁船的航向(方位角)是指在水平面上以漁船位置為中心,從該點的指北方向線起,依順時針方向到目標方向線之間的水平夾角,航向值在0°~360°之間。

圖2 2012年10月10日00: 00到24: 00 拖網船300585航向變化Fig.2 Heading changes of trawlers 300585 from 00: 00 to 24: 00 on October 10,2012

當漁船航向在 0°或 360°附近變化時,航向值會出現較大的變化幅度,圖 2中航向變化較大的有 B段、G段、I段,為了進一步分析航向的實際變化狀況,對航向角進行了差值計算。航向差是兩個相鄰時間,后一時間航向與前一時間航向的差值,差值結果中的正值表示航向順時針轉動,負值表示漁船航向逆時針轉動。圖2中航向差主要在0°左右變化,圖中 G段為捕撈狀態,但航向變化大,經過相鄰航向相減處理成航向差后與其他捕撈狀態時的航向差相似。為了避免漁網纏繞,拖網船作業時在3 min內一般不會出現較大的航向變動。

圖3是2012年10月10日00: 00~24: 00拖網船300585航速的變化,圖中的 A~I段按照狀態劃分,其分段與圖1和表1中相同,基本可以分為,拋錨的B段和I段,捕撈的C~H段、航行的A段,3種狀態。在C~H段的6次拖網作業中,每段都有航速的突變,段開始有航速峰值,段結束有航速的谷值,兩者之間是一個持續的較穩定航速,這是因為拖網作業特點有放網、拖曳網具、起網 3個階段。放網時漁船從尾滑道放出網具,然后漁船快速前進逐步放出曳綱,出現航速峰值的突變; 當曳綱放出預定長度后,漁船按預定的航向和航速拖網前進,航速較穩定;起網時漁船慢速前進收絞曳綱、網具,將其自尾滑道拖到甲板上,取出漁獲物,出現航速谷值的突變,航向差變化較大。

累計變化角度是航向差值的累加和,按照航向差順時針加,逆時針減的方式計算。累計變化角度可以反映出漁船方向在一段時間角度總體的變化值。圖3中拋錨B段和I段,累計角度變化較大,捕撈C~H段、航行 A段累計角度變化小,拖網船 300585 在2012年10月10日從00: 00~24: 00,累計角度是-1 066°,即船在這24 h內相對初始航向逆時針轉了近3圈。

圖3 2012年10月10日00: 00~24: 00拖網船300585船速與累計變化角度分布圖Fig.3 Speed distribution and cumulative change angle of trawlers 300585 from 00: 00 to 24: 00 on October 10,2012

由以上分析拖網船可以分為拋錨或慢速、作業、航行 3種狀態,結合航速和航向差可以判斷出這些狀態(表2),2012年10月10日00: 00到24: 00,拖網船300585作業時間15.9 h,航行2.5 h,拋錨5.6 h。

表2 拖網船300585作業狀態Tab.2 300585Boat Status

2.2拖網漁船航速統計

圖4和圖5是拖網船300585在2012年1~5月、9~12月及全年船位點記錄數量隨速度的變化,它是對2012年該船全年位置點的統計。由于漁船進漁港停泊時也發送船位數據,因此產生了很多值為0 m/s的數據,這些數據對分析沒有意義,統計中去掉了這些數據,6~8月是拖網的休漁期,統計中沒有統計這幾個月。

圖4 2012年1~5月點記錄數量隨速度的變化Fig.4 Number of point changes according to speed in 2012 from January to May

圖5 2012年9~12月及全年點記錄數量隨速度的變化Fig.5 Number of point changes according to Speed in 2012 from September to December and 2012 year

圖4和圖5曲線類似,主要存在3個峰,第一個在0~1 m/s,漁船處于拋錨或漂流; 第二個在1.3~1.9 m/s,漁船處于捕撈狀態; 第三個在3.3~4.3 m/s,漁船處于航行。

每個拖網漁船全年統計曲線與拖網船300585在2012年匯總曲線基本相同,利用拖網船 300585的全年統計數據與數據庫中其他漁船的統計數據計算相關系數R,可以判斷出拖網船,結果如表3,相關系數在 0.9時,可以正確判斷出拖網漁船 781艘,判斷錯誤的漁船59艘,正確率在93%,因此拖網船300585的統計數據具有代表性。

對2012年拖網船300585的100390個船位點,進行統計繪圖(圖4、圖5),其他拖網船全年點記錄數量隨速度的變化曲線相近,根據圖中第二個峰值位置,設定公式(1)中V1為1 m/s,V2為2.1 m/s。各船繪制航向差分布圖,與拖網船 300585相似,設定漁船拖網作業時的航向差 D1為-50°,D2為50°。

表3 拖網漁船判斷Tab.3 Trawler judgment

2.3累計捕撈時間專題圖

項目選用ArcGIS制作累計捕撈時間專題圖,以2012年10月10日的插值圖為例,把海域按0.1°×0.1°劃分格網,按照公式(1)計算格網中所有漁船在 2012 年10月10日0~24 h捕撈狀態的累計捕撈時間(單位h),作為捕撈努力量。根據拖網船航速、航向閾值提取出處于作業狀態的點,共提取到1 423艘拖網船處于捕撈狀態的點318 433個,合計拖網捕撈時間15 921 h。把點分布到0.1°×0.1°的格網中,計算每個格網中的累計捕撈時間,圖 6中的點是 0.1°×0.1°格網中心點(格),顏色變化代表累計捕撈時間的長短,圖中對格點值采用公式(6)進行了插值,生成2012年10月10日插值圖,反應拖網捕撈投入的分布趨勢面。

圖6 2012年10月10日累計時間捕撈強度分布Fig.6 Distribution of cumulative time of fishing intensity on October 10,2012

在傳統的漁場格網中,可以看出累計捕撈時間最高的漁場在舟外漁場和江外漁場南部,大沙漁場、長江口漁場、舟山漁場與魚山漁場交界的區域也有較高的累計捕撈時間分布。

3 討論與分析

國外學者基本都通過漁船的速度和角度區分漁船作業狀態,Skaar,J?rgensen等[16]根據船速分析巴倫支海的鱈魚拖網捕魚活動,Mills,Townsend等[17]在英國北海利用拖網漁船的速度和方向識別作業和航行狀態。另外有學者通過船位數據計算捕撈努力量,2005年Deng和 Dichmont等[18]用澳大利亞北部對蝦作業的船位數據估算拖網捕撈強度,還有學者利用船位數據研究捕撈行為,Walker,Beza等[19-20]利用貝葉斯模型研究捕撈過程。

拖網船有拋錨或慢速、作業、航行3種狀態,拖網作業的每個網次有放網、拖曳網具、起網 3個階段,通過分析拖網漁船的航速、航向差的特征可以判斷漁船的狀態,掌握漁船的作業規律,根據長時間的漁船航速特點可以獲取到漁船各狀態的航速閾值,再通過該值結合航向差能夠提取出漁船處于捕撈狀態的船位點,根據點的數量,處于捕撈狀態點的持續時間,計算出在某個區域累計捕撈時間,制作出累計捕撈時間格點圖,插值出累計捕撈時間變化趨勢圖。

拖網捕撈的累計時間越長,投入的功越多,累計捕撈時間在一定程度上能夠反映出捕撈努力量,即一個時間段內多艘拖網船投入捕撈小時數的累加。但是由于拖網船功率不全相同,單位時間投入功的大小存在差別,因此在掌握漁船功率數據的條件下,將進一步研究累計捕撈時間與功率的乘積,獲得累計捕撈(kW·h),把它作為捕撈努力量。拖網漁船在網具放置水層深度不同、捕撈魚種不同、功率不同的情況下作業的航速有些細微差別,因此下一步研究中將根據每艘船的航速閾值分別判斷漁船狀態。

下一步將從時空角度分析累計捕撈時間變化規律,比較累計捕撈時間與傳統捕撈努力量計算方式的優缺點,做進一步改進。通過編程實現累計捕撈時間專題圖的自動制圖,逐漸實現業務化,為漁政管理、漁業資源保護提供指導和借鑒。

致謝: 中國水產科學研究院東海水產研究所嚴利平研究員,張寒野、劉勇、馮春雷副研究員等為本文撰寫提出了積極建議,謹此致謝!

參考文獻:

[1]張衡,張勝茂.東南太平洋智利竹筴魚漁場及單位捕撈努力量的時空分布[J].生態學雜志,2011,30(6):1142-1146.Zhang Heng,Zhang Shengmao.Spatiotemporal distri-bution pattern of Chilean jack mackerel(Trachurus murphyi)fishing grounds and catch yield per unit effort in Southeast Pacific Ocean[J].Chinese Journal of Ecology,2011,30(6): 1142-1146.

[2]田思泉,陳新軍.不同名義CPUE計算法對CPUE標準化的影響[J].上海海洋大學學報,2010,19(2):240-245.Tian Siquan,Chen Xinjun.Impacts of different calculating methods for nominal CPUE on CPUE standardization[J].Journal of Shanghai Ocean University,2010,19(2): 240-245.

[3]方水美.福建沿海張網作業捕撈能力的計算分析[J].中國水產科學,2005,12(3): 321-328.Fang Shuimei.Calculated analysis on fishing capacity of swing net in Fujian coastal sea [J].Journal of Fishery Sciences of China,2005,12(3): 321-328.

[4]鄭國富.底拖網作業捕撈努力量標準化方法研究[J].福建水產,2000,85(2): 28-34.Zheng Guofu.An approach to the standardization of fishing effort of bottom trawl[J].Journal of Fujian Fisheries,2000,85(2): 28-34.

[5]Lee J,South A B,Jennings S.Developing reliable,repeatable,and accessible methods to provide highresolution estimates of fishing-effort distributions from vessel monitoring system(VMS)data [J].ICES Journal of Marine Science: Journal du Conseil,2010,67(6):1260-1271.

[6]Mullowney D R,Dawe E G.Development of performance indices for the Newfoundland and Labrador snow crab(Chionoecetes opilio)fishery using data from a vessel monitoring system[J].Fisheries Research,2009,100(3): 248-254.

[7]Stelzenmüller V,Maynou F,Bernard G,et al.Spatial assessment of fishing effort around European marine reserves: Implications for successful fisheries management[J].Marine Pollution Bulletin,2008,56(12):2018-2026.

[8]郭飚,薛元宏.北斗系統在海洋漁業信息化建設中的關鍵技術與實現途徑[J].現代漁業信息,2004,19(5):13-14.Guo Biao,Xue Yuanhong.Key technique and realized way of Beidou System(Staellite Navigation System)in the construction of marine fisheries information[J].Modern Fisheries Information,2004,19(5): 13-14.

[9]張勝茂,王曉璇,周為峰.基于漁船監測系統的近海捕撈水產品溯源[J].電腦開發與應用,2014,24(4):16-19.Zhang Shengmao,Wang Xiaoxuan,Zhou Weifeng.Offshore fishing aquatic products traceability based on vessel monitoring system[J].Computer Development and Applications,2014,24(4): 16-19.

[10]胡剛,馬昕,范秋燕.北斗衛星導航系統在海洋漁業上的應用[J].漁業現代化,2010,37(1): 60-62.Hu Gang,Ma Xin,Fan Qiuyan.The apllications of COMPASS Navigation Satellite System to marine fishing industry[J].FisheryModernization,2010,37(1):60-62.

[11]居禮.北斗衛星導航系統在海洋漁業的應用[J].衛星與網絡,2013(3): 16-22.Ju Li.Application of Beidou satellite navigation system in marine fishery[J].Satellite and Network,2013(3):16-22.

[12]張勝茂,樊偉.海洋次表層FIDW溫鹽影像插值算法[J].計算機工程與應用,2012,48(26): 205-209.Zhang Shengmao,Fan Wei.FIDW algorithm for temperature and salinity image interpolation in sub-surface ocean[J].Computer Engineering and Applications,2012,48(26)205-209.

[13]湯國安,楊昕.ArcGIS地理信息系統空間分析實驗教程[M].北京: 科學出版社,2006.Tang Guoan,Yang Xin.ArcGIS geographic information system spatial analysis experiment course[M].Beijing:Science Press,2006.

[14]Watson D F,Philip G M.A refinement of inverse distance weighted interpolation[J].Geoprocessing,1985,2:315-327.

[15]Philip G M,Watson D F.A precise method for determining contoured surfaces[J].Australian Petroleum Exploration Association Journal,1982,22(1): 205-212.

[16]Skaar K L,J?rgensen T,Ulvestad B K H,et al.Accuracy of VMS data from Norwegian demersal stern trawlers for estimating trawled areas in the Barents Sea[J].ICES Journal of Marine Science,2011,68(8):1615-1620.

[17]Mills C M,Townsend S E,Jennings S,et al.Estimating high resolution trawl fishing effort from satellite-based vessel monitoring system data [J].ICES Journal of Marine Science: Journal du Conseil,2007,64(2): 248-255.[18]Deng R,Dichmont C,Milton D,et al.Can vessel monitoring system data also be used to study trawling intensity and population depletion? The example of Australia's northern prawn fishery [J].Canadian Journal of Fisheries and Aquatic Sciences,2005,62(3):611-622.

[19]Walker E,Bez N.A pioneer validation of a state-space model of vessel trajectories(VMS)with observers' data[J].Ecological Modeling,2010,221(17): 2008-2017.

[20]Bez N,Walker E,Gaertner D,et al.Fishing activity of tuna purse seiners estimated from vessel monitoring system(VMS)data[J].Can J Fish Aquat Sci,2011,68(11): 1998-2010.

Received: Feb.17,2014

(本文編輯: 劉珊珊)

Key words:Beidou satellite; vessel monitoring system; fishing effort; heading; speed

Method of extracting trawling effort based on vessel monitoring system

ZHANG Sheng-mao1,ZHANG Heng1,TANG Feng-hua1,FAN Wei1,HUANG Hua-wen2
(1.Key Laboratory of East China Sea & Oceanic Fishery Resources Exploitation and Utilization,Ministry of Agriculture,China,East China Sea Fisheries Research Institute,Chinese Academy of Fishery Sciences,Shanghai 200090,China; 2.Shanghai Ubiquitous Navigation Technologies Ltd.,Shanghai 201112,China)

Abstract:The purpose of this study was to extract trawling effort from the statistics of a vessel monitoring system.Using speed data,the results show that there are three peaks in trawling effort.The second peak is in the fishing state,and generally occurs at speeds between 1 m/s to 2.1 m/s.The heading difference is between -50 and 50 degrees.Trawlers have three states,including low speed,fishing,and sailing,which can be distinguished by their thresholds of speed and heading differences.We extracted a total of 318 433 fishing-state points from 1 423 fishing vessels,and the total trawling time was 15 921 h.We used an inverse distance weighting interpolation method to generate change trends in the fishing intensity distribution.Fishing effort is an important reference value in fishery resource research.Compared with traditional fishing effort calculations,this method is characterized by real-time,large-scale,rapid,and high resolution.The results are applicable in the conservation of fisheries resources.

中圖分類號:S975

文獻標識碼:A

文章編號:1000-3096(2016)03-0146-08

doi:10.11759/hykx20140217002

收稿日期:2014-02-17; 修回日期: 2014-06-28

基金項目:中央級公益性科研院所基本科研業務費專項資金項目(東海水產研究所 2014T13); 十二五國家科技支撐計劃項目(2013BAD13B01);上海市科學技術委員會資助項目(12511501200)

作者簡介:張勝茂(1976-),男,河北吳橋人,副研究員、博士,主研方向: 漁業遙感與地理信息,電話: 021-65682395,E-mail: ryshengmao@126.com; 通信作者,樊偉,男,河南鄭州人,博士,研究員,研究方向為漁業遙感與地理信息,電話: 021-65680117,E-mail:dhyqzh@sh163.net

主站蜘蛛池模板: 综合色婷婷| 成人一级免费视频| 亚洲男人的天堂视频| 国产粉嫩粉嫩的18在线播放91| 国产av无码日韩av无码网站| 久久精品国产国语对白| 99热国产在线精品99| 99热这里都是国产精品| 国产成人午夜福利免费无码r| 天天综合网亚洲网站| 啪啪永久免费av| 日韩国产一区二区三区无码| 国产视频只有无码精品| 欧美激情,国产精品| 国产av色站网站| 性喷潮久久久久久久久| 亚洲国产高清精品线久久| 国产精品流白浆在线观看| 亚洲精品第一页不卡| 亚洲系列中文字幕一区二区| 国产欧美综合在线观看第七页 | 久久久久亚洲AV成人网站软件| 国产高清在线观看| 亚洲资源站av无码网址| 亚洲精品国产综合99| 欧美国产综合色视频| www成人国产在线观看网站| 亚洲人成在线精品| 九色视频线上播放| 亚洲永久色| 成人综合网址| 亚洲精品国产首次亮相| 激情爆乳一区二区| 亚洲人成网站观看在线观看| 亚洲欧美另类中文字幕| 国产日韩欧美视频| 亚洲Av激情网五月天| 亚洲精品国产日韩无码AV永久免费网| 美臀人妻中出中文字幕在线| 国产不卡国语在线| 成色7777精品在线| 亚洲女人在线| 99热在线只有精品| 久久99精品久久久大学生| 久久精品视频亚洲| 国产av剧情无码精品色午夜| 国产成在线观看免费视频| 亚洲国产亚洲综合在线尤物| 99久久国产综合精品2020| 国产精品嫩草影院视频| 国产免费久久精品99re丫丫一 | 日本高清免费不卡视频| 成人中文在线| 久久午夜影院| 国产又色又刺激高潮免费看| www中文字幕在线观看| 国产亚洲精品在天天在线麻豆 | 幺女国产一级毛片| 99国产在线视频| 国外欧美一区另类中文字幕| 欧美日韩一区二区在线免费观看| 亚洲91精品视频| 国产91视频观看| 一级成人a毛片免费播放| 99资源在线| 久久亚洲天堂| 欧美日本中文| 国产综合网站| 2020国产精品视频| 欧美日韩精品一区二区视频| 国产精品欧美激情| 狼友视频国产精品首页| 日本国产精品一区久久久| 国产欧美性爱网| 香蕉久人久人青草青草| 色综合久久久久8天国| 亚洲欧洲自拍拍偷午夜色| 色哟哟国产精品| 欧美区一区二区三| 夜夜操国产| 欧美亚洲综合免费精品高清在线观看| 欧美中文字幕一区二区三区|