?

GNSS 漂流浮標實時在線波浪測量技術及其軟件實現

2023-12-27 15:12王亞彬劉焱雄喬方利蔣暑民王巖峰
海岸工程 2023年4期
關鍵詞:波高浮標波浪

王亞彬,劉 楊,2,劉焱雄,2*,喬方利,蔣暑民,王巖峰

(1. 自然資源部 第一海洋研究所, 山東 青島 266061;2. 自然資源部 海洋測繪重點實驗室, 山東 青島 266061)

波浪是一種重要且復雜的海洋水文要素。有效、穩定、長期地監測波浪,對海洋水文和氣象研究、海洋資源開發利用,以及沿海農業、漁業生產等具有重要意義(戴洪磊等, 2014)。目前,波浪浮標是實現海洋波浪定點長期觀測的主要設備,波浪觀測技術的進步能夠提高觀測效率和觀測精度。所以,研究和發展波浪浮標及波浪觀測技術是提高波浪觀測質量、組建波浪觀測網、豐富波浪觀測資料、促進波浪觀測行業發展的必然需求(毛祖松, 2007)。

海洋觀測浮標作為一種現代化海洋監測設備,因其具有定點、全天候長期觀測的特點,在近些年海洋探索的進程中蓬勃發展。發達國家在這一方面的研發起步較早,起始于20 世紀30 年代末,至今一直處于領先位置。目前,在波浪測量領域,美國、英國、加拿大、荷蘭、挪威等國家都在研發和使用波浪浮標,并布放于各自臨海區域組建成網(王軍成, 1998)。國外有代表性的波浪浮標主要有荷蘭Datawell 公司的Directional Waverider 系列和法國的TRIAXYS 系列浮標。然而,我國海洋觀測浮標的研制于20 世紀60 年代才起步,后續雖然已研制出各種類型的觀測浮標及裝備,并且覆蓋了相當一部分海域,但搭載的儀器性能、測量精度和工作穩定性等方面,跟國外相比還有一定差距(毛祖松, 2007),這也在一定程度上阻礙了海洋觀測工作的有效開展,進而影響到海洋資料的獲取和積累。從國內海洋觀測浮標研制開始,國內多家企業、高校和科研院所等對波浪測量領域進行了研究探索,并取得寶貴的經驗和成果(王亞洲等, 2006; 唐原廣等, 2008; 劉國棟, 2011; 齊占輝, 2019; 王斌等, 2021; 張新文等, 2023)。

關于國內波浪浮標,比較有代表性的主要有3 種(毛潤雨, 2021):SZF 系列波浪浮標、SBF3 系列測波浮標和OSB 系列浮標。其中,SBF3 系列測波浮標還可利用橡膠彈性材料的錨系繩纜來進一步接近波浪的真實運動軌跡(孫金偉等, 2012)進而提高觀測質量。截至2016 年底,測波浮標仍為國內主要的業務化運行現場自動測波儀器,我國海洋站中約87.5%的站點使用小型測波浮標進行波浪測量(常怡婷等, 2021)。

基于全球導航衛星系統(Global Navigation Satellite System, GNSS)的測波浮標取消了加速度傳感器等元件,浮標位置與波浪要素依靠同一定位模塊獲得,以使浮標成本更低且更為小巧輕便易于布放。由于GNSS 模塊逐歷元定位,不存在系統誤差的累積,不需要定期調校,所以,GNSS 測波浮標特別適合無人值守的長期監測。de Vries 等(2003)證明基于全球定位系統(Global Positioning System, GPS)多普勒測速原理的Datawell DWR-G 型號波浪浮標的測量精度與傳統加速度計式波浪浮標相當,且能更好地用于研究遠洋涌浪、港口或航道共振等形成的長周期波。同樣,基于GPS 技術的漂流式波浪浮標Spotter,以其低成本、易操作、實時傳輸和緊湊的外形得以組成“實時海洋傳感器陣列”,實現了在全球范圍內布放(Raghukumar et al, 2019; Pierik et al, 2020)。但傳統的基于動態后處理(Post Processed Kinematic, PPK)技術、實時動態定位(Real Time Kinematic, RTK)技術、精密單點定位(Precise Point Positioning, PPP)技術及多普勒測速技術的GNSS 波浪測量技術和波浪浮標,總會受到基站距離、通信成本、時效性及結果精度等方面的限制。本文研究了基于GNSS 載波相位歷元差分算法的GNSS 漂流浮標實時在線波浪測量技術及其實現方法,給出了該技術硬件環境的選型、組成,以及實時在線軟件的設計思路和工作原理,并進行軟件實現。研究的GNSS 載波相位歷元差分算法的漂流浮標實時在線波浪測量技術及其軟件,能夠在基于廣播星歷、無需基站及高額改正服務的基礎上進行高精度實時波浪測量。并且,結合我國北斗衛星導航系統的發展,研究與北斗技術相結合的漂流浮標實時在線波浪測量技術,對實現波浪浮標國產化、打破國外設備壟斷局面具有重要意義。

1 GNSS 實時波浪測量方法

1.1 GNSS 測速方法

基于載波相位歷元差分原理的GNSS 實時波浪要素測量技術,僅使用廣播星歷即可得到高精度速度(單瑞, 2010; 劉會等, 2020; 劉楊等, 2020; Liu et al, 2022),差分后的方程為:

式中: λ為載波波長;ΔΦsr和esr分別為衛星s到接收機r 的相鄰歷元載波相位觀測值差和單位矢量;Δξr為接收機相鄰歷元間的位置變化量; c為光速;Δδtr、Δδts分別為相鄰歷元間的接收機鐘差、衛星鐘差;Δpsr為相鄰歷元間包括衛星軌道、電離層延遲、對流層延遲、相位中心變化、相位纏繞、地球自轉效應和相對論效應等在內的綜合誤差,其中由廣播星歷計算衛星軌道改正,由雙頻無電離層組合消除一階項進行電離層改正,由氣象數據及投影函數進行對流層改正,由天線模型計算相位中心變化,相位纏繞、相對論效應、地球自轉效應等誤差改正項也可以由相應模型算出;Δεsr為其他殘余誤差項與噪聲;v為三維速度; Δt為相鄰歷元之間的時間間隔(劉楊等, 2020; 田力等, 2021;Liu et al, 2022)。

實際上,當求得解算速度并剔除粗差后,受海況、浮標位置、信號誤差和其他噪聲的影響,速度積分的位移結果存在趨勢項(文圣常, 1984; 俞聿修, 2011; Herbers et al, 2012),因此,在剔除速度粗差并得到速度一次積分的位移結果后,還要去除線性趨勢項,進而得到波浪位移信息,在此過程中考慮到波浪周期的范圍(主要位于1~30 s),也可利用高通濾波來消除0.03 Hz 以下的低頻噪聲(Liu et al, 2022)。

1.2 波浪要素反演

在波浪要素反演過程中,須計算3 種主要的波浪要素,分別為波高、波周期和波向。

1.2.1 波高和波周期

海浪可視為各態歷經的隨機過程(俞聿修, 2011),因此可利用統計手段取樣本中足夠長的一段數據來分析總體數據的特性。通常利用上跨零點法對樣本中的不規則波進行定義,故在上跨零點法定義不規則波的基礎上,利用統計手段求取波高、波周期的步驟如下。

首先,將波浪高度序列Hi按照波高值由大到小依次排列:H1,H2,···,HN-1,HN;同理,波周期序列Ti排列為:T1,T2,···,TN-1,TN。然后,以前1/3 大波波高作為有效波高Hs,波列中所有波浪的波周期平均值作為平均周期Tm,則Hs和Tm定義(文圣常, 1984; 俞聿修, 2011)為:

1.2.2 波向

波向可利用方向譜提取。利用傅里葉變換法求解方向譜,計算速度快、不易發散且算法穩定,求解的主要步驟如下(Benoit, 1992; 毛潤雨, 2021)。

首先,可由波浪特性觀測值的傅里葉變換內積求得時間T內的交叉譜Smn(f):

式中:f為頻率;Fm(f)為波浪特性m的傅里葉變換的共軛;Fn(f)為波浪特性n的傅里葉變換;Cmn(f)和Qmn(f) 分別為同向譜和正交譜。

然后,海浪的二維方向譜可以表示為有限階數的傅里葉級數展開:

式中: θ為傳播方向;A0、A1、A2、B1和B2為海浪二維方向譜的前五系數,可分別由6 個交叉譜求得,如:

最后,二維海浪譜S(f,θ)又為一維海浪譜S(f)與方向分布函數D(θ,f)的乘積,即:

一維海浪譜S(f)可由二維方向譜S(f,θ)在角度 θ上積分得到,因此得知二維方向譜和一維海浪譜后,便可以得出方向分布函數D(θ,f),進而得出成分波方向。

2 GNSS 漂流浮標實時在線波浪測量技術實現

GNSS 漂流浮標實時在線波浪測量技術的實現主要依靠實時在線解算軟件及其硬件運行環境。實際測量時,海上測波浮標搭載的傳感器等采集模塊將測量得到的信號轉換成觀測值后,發送至浮標系統的主控模塊進行所需參數的解算和緩存,解算完成后組成特定格式的交換報文由通信模塊發送至岸基設備,完成一次觀測任務。

實時在線解算軟件作為實時在線波浪測量技術的關鍵,負責任務調度及協調各模塊工作,能夠基于GNSS 載波相位歷元差分算法獲得的高精度三維載體速度和位置進行波浪要素的反演(劉楊等,2020; Liu et al, 2022)。此外,基于GNSS 實時精密單點定位算法還可以在位置估算的同時得到天頂對流層延遲,進一步結合浮標實測氣溫、氣壓或者數值天氣預報模型來反演得到大氣水汽含量(劉楊等, 2022)。因此,浮標所采集物理量主要包括浮標三維位置、波高、波周期和波向,根據GNSS 算法及傳感器的不同還可獲得天頂對流層延遲,以及海水溫鹽和大氣溫壓等其他海氣環境參數。

硬件環境作為軟件運行的載體,負責保證高效的數據采集和數據傳輸,提供充足的計算資源和緩存空間,是GNSS 漂流浮標實時在線波浪測量技術的外在支撐。在滿足技術要求的前提下,硬件環境的合理選配有助于降低設備功耗,提高運行壽命。浮標部分硬件環境的結構示意圖如圖1所示。岸基部分主要負責接收、分析、存儲浮標回傳數據以及向浮標發送指令,由普通主機設備搭配信息接收模塊和數據管理軟件即可完成消息的接收與管理。

圖1 GNSS 漂流浮標實時在線波浪測量技術硬件環境結構Fig. 1 Hardware environment structure of the real-time online wave measurement technology of GNSS drifting buoy

2.1 硬件設計

2.1.1 主控板卡

GNSS 漂流浮標設計采用商業級低功耗主控板,板卡內核基于ARM Cortex-A7 架構,采用2 組優質80 Pin 板對板連接器加強對外擴展,支持8 路UART、2 路以太網以及2 路CAN 等工業級總線接口,支持外接TF(Trans Flash)卡。板載處理器能夠滿足實時采集和處理的任務需求,多樣的多路通信接口符合與多傳感器的外部通信及后續的擴展需要,可擴展存儲空間則保證了系統、軟件以及長時間數據的緩存及存儲。

2.1.2 GNSS 板卡

GNSS 板卡在GNSS 觀測中起到關鍵作用。綜合考慮觀測質量、成本、國產化以及北斗PPPB2b 精密單點定位改正信號(涂滿紅等, 2022),可設計采用國產全系統全頻GNSS 板卡。該板卡支持BDS-3 全球信號,以及BDS-2、GPS、GLONASS、Galileo、IRNSS、QZSS 和SBAS 等衛星信號,支持L-Band 星際增強,且支持抗干擾;標稱測量準確度可達偽距≤10 cm,載波相位精度≤0.005c(c為載波波長,單位為m);采樣頻率最高支持20 Hz,支持RTCM3、RTCM2 以及NMEA-0183 等標準數據格式的輸出。

板卡設置輸出通用二進制RTCM3 格式信息,通過串口向核心板輸入GNSS 觀測數據與星歷,數據流內部配置的輸出報文如表1 所示。

表1 GNSS 板卡所配置的RTCM3 消息類型和內容Table 1 RTCM3 message type and contents configured for GNSS board

2.1.3 通信模塊

通信部分主要包括浮標系統內部通信及浮標與岸基的通信。板卡級別的內部通信設計采用標準通信串口進行異步通信,將主控板分別與各類數據的采集模塊及對外通信模塊鏈接。對外通信模塊主要負責將主控板處理和解算完成后組建的消息報文回傳至岸基接收設備,或者將岸基指揮處的運行指令下達至浮標??紤]到遠海使用場景,排除目前近岸短距離波浪浮標的2 種通信方式:海上無線短波通信和岸基移動通信(黨超群等, 2016),可采用以北斗短報文通信為主,銥星、天通等衛星通信為輔的設計,通過衛星通信方式來保證全球范圍內的消息回傳。

2020 年,我國第三代北斗衛星導航系統——北斗三號衛星導航系統組網完成。北斗系統別具特色的短報文通信服務具備低成本、廣覆蓋、高可靠和隨遇接入等特點,逐漸成為現今我國海上實時數據傳輸的重要方式,更成為許多浮標實時通信的選擇(楊軍平等, 2019; 姬生月等, 2021),可作為浮標的主要通信手段。天通一號衛星移動通信系統同樣是我國自主研發的系統,與北斗系統共同為我國打破國外技術封鎖,實現自主可控、國產化、低成本海洋設備提供了技術基礎(李聽聽等,2021),保證了通信的國產化及信息安全性,可作為北斗短報文通信的補充。發展較早的國外銥星衛星通信系統的優勢在于:其由66 顆環繞地球的低軌衛星組網,采用極地軌道和星際鏈路技術提高了通信的信號強度和可靠性(鐘健瑜, 2008),能夠在特殊地區或特殊環境下開展有效通信,但是通信費用較高,可考慮作為個別補充手段。

不同通信系統具有不同等級的報文傳輸限制。最新的北斗三號系統區域短報文和全球短報文服務的單次報文最大長度分別為14 000 bit 及和560 bit,分別約1 000 漢字和40 漢字。若采用北斗三號系統短報文通信為主要通信手段,可依據最大報文長度僅為560 bit 的全球短報文服務設計報文,報文回傳GPS 周、GPS 周內秒、經度、緯度、高度、波浪高度、波浪周期、波浪方向以及其他環境參數等,對于更多的信息回傳需求,還可通過壓縮現有參數字節長度或增加消息包內報文數目來保證信息傳輸完整。

2.2 軟件設計

2.2.1 軟件架構設計

波浪實時在線測量軟件是利用C 語言開發的。軟件占用資源少,運行速度快,適于低級硬件平臺的移植開發工作。軟件采用命令行用戶交互(Command User Interface, CUI)方式,運行在無圖形桌面的主控板系統中,實時解算的配置選項可通過讀入配置文件或命令行手動設置,整個程序采用多線程機制進行實時處理,整個軟件的運行流程如圖2 所示。

圖2 軟件流程示意圖Fig. 2 A schematic diagram of the software flow

程序開始運行后,首先進入主線程做參數配置,并開啟各子線程進行任務調度。當主控板得到GNSS 板卡的輸入數據后,對數據解碼和預處理,待一次采集完成后,由廣播星歷和觀測值求取當前時刻載體速度,并將速度做一次積分,進而得到位移。之后,以設計的滑動窗口(例如窗口大小為30 min)進行滑動,對位移時間序列做去趨勢或高通濾波處理,進而求取波浪要素的結果。在求取波浪要素的同時,其他子線程負責協同完成信息回傳、時間同步和數據緩存等工作。

2.2.2 實時性能分析

實時處理軟件采用多線程同步運行的方式,多線程包括數據處理線程和報文組建發送線程等,實時處理的耗時由其中耗時最長的線程決定。在單核CPU 主頻1.8 GHz,內存512 MB 的虛擬機環境下進行處理時間測試,結果表明,數據處理線程運行3 000 次的單次平均運行時間為0.08 s,小于0.2 s 的GNSS 數據采樣間隔,報文組建發送線程運行10 次的單次平均運行時間為0.05 s,軟件滿足實時運行需求。實際海上試驗表明,在浮標硬件配置的實際運行環境中,在線軟件運行穩定,未出現數據處理速度與獲取速度不匹配和延遲的情況,軟件滿足實時運行需求。

3 試驗數據

為驗證GNSS 漂流浮標實時在線波浪測量技術及其軟件的可行性和精確性,分別在近海和遠海開展1 次海試試驗,并對測試浮標實時回傳結果進行比對驗證。

2022 年11 月7 日至10 日,在青島近岸海域(120°44′E, 36°14′N)進行近海海試試驗,稱為試驗一。試驗設備包括2 套內置實時在線軟件的測試浮標(編號分別為1 號和2 號)以及1 套Datawell 公司的Directional Waverider 波浪浮標(型號為DWR-G4)。3 套浮標直徑均為40~50 cm,布放現場如圖3 所示。海試比測波浪參數包括有效波高Hs、平均周期Tm以及主波向D。

圖3 2022 年青島近海2 套測試浮標(1 號、2 號)與1 套波浪浮標(DWR-G4)布放現場Fig. 3 Site for the testing of 2 GNSS buoys (No.1 and No.2) and 1 wave buoy (DWR-G4)in the offshore area of Qingdao in 2022

2022 年2 月,選擇南海及印度洋相關海域進行遠海試驗,稱為試驗二。在自然資源部第一海洋研究所遠海航次中,于南海及印度洋相關海域布放4 套測試浮標,編號分別為3 號、4 號、5 號和6 號,測試浮標布放海域及軌跡如圖4 所示。其中3 號、4 號和5 號浮標設置為連續采集工作模式,6 號浮標設置為間隔采集工作模式,即每間隔6 h 采集1 次,每次采集時間為1 h。

圖4 2022 年4 套遠海測試浮標(3 號、4 號、5 號和6 號)的軌跡Fig. 4 Tracks of 4 sets of the test buoys (No.3, No.4, No.5 and No.6) in the open sea in 2022

4 結果與分析

4.1 近海試驗結果與分析

對比青島近海海試(試驗一)中1 號、2 號測試浮標實時回傳結果與Datawell DWR-G4 浮標事后導出結果(圖5)可知,總體而言,3 套浮標的有效波高、平均周期及主波向變化趨勢基本一致。2套測試浮標與DWR-G4 波浪要素差值的絕對值最大值、差值的均方根結果(表2)顯示:在3 d 左右的數據對比中,1 號測試浮標與DWR-G4 有效波高差值的絕對值的最大值為0.23 m,差值均方根為0.06 m;2 號測試浮標與DWR-G4 有效波高差值的絕對值最大值為0.18 m,差值均方根為0.04 m;1 號、2 號測試浮標與DWR-G4 平均周期差值的絕對值最大值分別為1.43 s 和1.12 s,均方根分別為0.48 s 和0.49 s。對比分析近海試驗中2 套測試浮標實時回傳結果與DWR-G4 事后導出結果表明,2套測試浮標要素反演性能良好,其中波高測量準確度符合主流波浪浮標技術指標,即0.1 m+5%H,H為波高(常怡婷等, 2021)。

表2 近海海試實時回傳有效波高與平均周期反演結果Table 2 The retrieved results of the significant wave height and mean wave period that returned real-time in the offshore test

圖5 2022 年青島近海試驗測試浮標實時回傳結果與DWR-G4 對比Fig. 5 Comparison between the results returned real-time by the GNSS test buoys and those by the DWR-G4 buoy in the offshore test in Qingdao in 2022

4.2 遠海試驗結果與分析

因現場無其他設備比測,采用實時在線回傳數據與歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecasts, ECMWF)再分析波浪數據進行對比。ECMWF 再分析對比數據是1950 年1 月至今全球氣候的第五代大氣再分析數據集(https://cds.climate.copernicus.eu),本次試驗選取數據集中風浪、涌浪共同作用下的有效波高作為本文有效波高Hs的參考值,平均周期作為本文平均周期Tm的參考值,數據空間分辨率為0.25°,時間分辨率為1 h,區域范圍(25°~124°E,66°S~30°N)基本覆蓋遠海試驗(試驗二)4 套測試浮標在南海及印度洋的布放區域,時間從2 月到7 月。

2022 年遠海試驗(試驗二)4 套測試浮標3 號、4 號、5 號、6 號實時回傳結果與ECMWF 再分析數據的對比結果(圖6)顯示,測試浮標與再分析數據的有效波高Hs(圖6a~圖6d)和平均周期Tm(圖6e~圖6h)的變化趨勢基本一致,說明測試浮標能夠在遠海海域實時回傳波浪要素結果反映當地海浪變化趨勢。

圖6 2022 年遠海試驗4 套測試浮標(3 號、4 號、5 號和6 號)實時回傳結果與ECMWF 產品的比較Fig. 6 Comparison between the results returned real-time by 4 sets of the test GNSS buoys (No.3, No.4, No.5 and No.6)and those by the ECMWF products in the open sea test in 2022

2022 年遠海試驗4 套測試浮標實時回傳數據與ECMWF 再分析數據波浪要素之間差值的絕對值最大值、差值的平均值、差值的均方根結果(表3)顯示,4 套測試浮標與對比數據的有效波高Hs差值的均方根最小為0.21 m,最大為0.47 m,平均周期差值的均方根最小為1.1 s,最大為3.4 s。平均周期Tm差值的平均值均為負值,結合平均周期(圖6e~圖6h)變化趨勢,發現測試浮標反演的平均周期多數小于對比數據的平均周期??紤]到對比數據的精度和時空窗口限制,不能精確評估測試浮標的結果精度。后續需要進一步開展與現場實測數據的比測工作。

表3 遠海海試實時回傳結果Table 3 Results returned real-time in the open sea tests

5 結 語

傳統高精度GNSS 方法受時間、距離和成本的制約,本文主要研究了GNSS 漂流浮標實時在線波浪測量技術及其軟件實現方法。結合載波相位歷元差分測速模型及波浪要素反演方法,給出該技術硬件環境的選型、結構設計和實時在線軟件的實現策略,并進行了軟件實現和實時處理速度測試。選擇青島近海(120°44′E, 36°14′N)、南海及印度洋(25°~124°E, 66°S~30°N)分別作為近海試驗(試驗一)和遠海試驗(試驗二)的研究區域,對比分析了1 號、2 號(近海)測試浮標實時回傳結果與DWR-G4 事后導出結果,以及3 號、4 號、5 號、6 號(遠海)測試浮標實時回傳結果與ECMWF再分析產品的開源數據結果,從而驗證了GNSS 漂流浮標實時在線波浪測量技術及軟件的可行性和精確性。本文主要結論如下。

1)硬件環境設計以主控板為核心,通過標準串口鏈接GNSS 板卡、通信模塊及其他傳感器,主控板讀取GNSS 板卡及各傳感器的觀測數據,待數據解算完畢后組成特定報文,通過衛星通信技術發送至岸基接收設備。該硬件環境的組成和結構能夠滿足GNSS 漂流浮標實時在線波浪測量的任務需求。

2)實時在線波浪測量軟件采用多線程機制。主線程負責運行參數設置、任務調度;各子線程同步運行,協調完成高精度位置和速度的解算、波浪要素的解算、報文回傳、數據緩存以及時間同步等工作。軟件實時處理速度滿足實時采集和解算的任務需要。

3)近海1 號和2 號測試浮標實時回傳結果與Datawell DWR-G4 事后導出結果表明,有效波高差值的均方根分別為0.06 m 和0.04 m,波高測量誤差優于主流測波浮標產品技術指標(0.1 m+5%H,H為波高)。周期和方向的量值及變化趨勢與對比的浮標結果基本一致,技術指標與國際主流產品相當。

4)遠海4 套測試浮標(3~6 號)的實時回傳結果與ECMWF 再分析產品比較結果顯示,二者有效波高和平均周期變化趨勢吻合良好,能夠反映當地海浪實時變化趨勢。相較于傳統高精度GNSS 方法,GNSS 漂流浮標實時在線波浪測量技術及其軟件無需基站配合,能夠克服測量范圍限制和避免差分改正服務費用。該技術只需衛星廣播星歷,即可在保證實時性的同時獲取較高的測速精度,進行準確的波浪要素反演。

多次近海及遠海測試結果驗證了GNSS 漂流浮標實時在線波浪測量技術及其軟件實現的可行性和精確性。但是,目前GNSS 漂流浮標測波仍有很多改進空間,如減少軟件資源的占用以降低硬件配置需求,改進算法精度以提高反演結果的可靠性等。

猜你喜歡
波高浮標波浪
受了委屈的浮標君
受了委屈的浮標君
受了委屈的浮標君
受了委屈的浮標君
波浪谷和波浪巖
珊瑚礁地形上破碎波高試驗研究
基于漂流浮標的南大洋衛星高度計有效波高研究
非平整港池的多向不規則波試驗研究
波浪谷隨想
海堤設計波高計算
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合