一种基于计数比例加权平均的波浪方向浮标主波向计算方法与流程

文档序号:17442826发布日期:2019-04-17 05:01阅读:377来源:国知局
一种基于计数比例加权平均的波浪方向浮标主波向计算方法与流程

本发明属于海洋监测技术领域,具体涉及一种基于计数比例加权平均的波浪方向浮标主波向计算方法。



背景技术:

当今世界,科学与经济技术迅猛发展,人口急剧增加,地球上陆地自然资源日趋紧张,能源和矿产资源匮乏成为如今各国可持续发展所需解决的重要问题之一。其中蕴藏量丰富的海洋已经成为目前研究最多的科学领域之一,并且已经获得了一些可观的成果。海洋能源包括潮汐能、波浪能、海水温差能、海流能及盐度差能[1]。因此,发展与开发海洋能源,成为了世界上众多沿海国家的战略性选择。海浪监测是人类认识、研究、开发利用海洋的重要任务之一,现代航海、造船、军事、海洋工程、灾害监测预警等领域,都需要获取高精度的实时海浪资料,用来描述海浪特征的量,如波高、波向、周期、相速、波长等,统称为海浪要素[2]。其中海浪的波向参数特性的研究对于广泛的海洋工程和海洋科学研究来说是十分重要的,例如对海洋工程的开发、泥沙运输计算和海浪预测模型的验证都具有重大意义[3]。由于海洋开发己经向较深的水域挺进,因此,海洋工程设计中采用过份冗余的安全系数设计所付出的代价越来越难以承受,在设计要求中所使用的环境数据必须体现海洋要素与波浪传播方向的关系,而不能笼统地假设在各个方向上都会遇到最恶劣的情况,例如利用波浪能量进行发电,在可行性研究中需要考虑海浪波向资料;对海岸开发来说,在设计诸如防波堤之类的建筑物时,波向资料亦十分重要[4]。见参考文献[1-4](王传崑,卢苇.海洋能资源分析方法及储量评估[m].海洋出版社,2009.中国海洋大学,2007.齐勇,闫星魁,郑姗姗,等.海浪监测技术与设备概述[j].气象水文海洋仪器,2015,33(3):113-117.kuikaj,vanvleddergp,holthuijsenlh.amethodfortheroutineanalysisofpitch-and-rollbuoywavedata[j].journalofphysicaloceanography,1988,18(7):1020-1034.孙仲汉.marex公司的波向浮标研究[j].海洋技术学报,1987(2):63+94-98.)

目前国内外进行海浪的波向数据观测主要是采用波浪方向浮标的方式。波浪方向浮标是一种无人值守的自动、定点、定时的对海面波浪的各要素进行遥测的小型浮标测量系统。在测量波向方面,因为浮标在随波浪升沉的同时也在进行水平方向的运动,所以能够反映出海浪方向的变化,且不易受其它因素的影响[5]。见参考文献[5](朱光文.海洋监测技术的国内外现状及发展趋势[j].气象水文海洋仪器,1997(2):1-14.)

国内外公开文献上能够查到的利用波浪浮标进行波向测量的仪器原理上主要分成如下两类[6][7]:第一类是利用浮标体内的三轴加速度传感器和电子罗盘,测得浮标体的方位角和水平两个垂直方向的瞬时位移,进而计算出主波向,国产的sbf型波浪浮标即属于该类型;第二类是利用浮标内的波高倾斜一体化传感器和方位传感器,分别测得浮标体的瞬时方位角、横倾角和纵倾角信息,进而计算出主波向,国产的szf型波浪浮标即属于该类型。见参考文献[6-7](唐原广,康倩.波浪浮标测波方法比较[j].现代电子技术,2014,v.37;no.422(15):121-122.毛祖松.我国近海波浪浮标的历史、现状与发展[j].海洋技术学报,2007,26(2):23-27.)

针对第一类工作原理的浮标,目前能够从公开文献上查到的测量波向方法如下:以国产的sbf型波浪浮标为例,该浮标是由山东省科学院海洋仪器仪表研究所研制完成的,测量波向方法是采用浮标体内的三轴加速度传感器和电子罗盘相结合的方式,在一个波浪周期内的跨零点处,计算浮标的水平矢量合成加速度方向作为单个波浪周期的参考方向。以加速度x轴正轴和y轴正轴建立相对正交坐标系,计算出水平矢量合成加速度与x轴的夹角即为相对波向,再结合电子罗盘计算得到的方位角就可得到波浪的地理波向[8][9]。见参考文献[8-9](赵杰,惠力,初士博,等.基于三轴加速度的波浪测量技术研究[j].海洋技术学报,2015,34(5):66-70.刘国栋.波浪浮标数据处理方法研究[d].天津大学,2012:49-50.)

国内针对第二类工作原理的浮标,目前能够从公开文献上查到的测量波向方法如下:以国产的szf型波浪浮标为例,该浮标是由中国海洋大学研制生产的,测量波向方法与国外文献公布的方法类似,首先根据上跨零点来确定单个波浪周期,再找到每一个单个波浪周期的下跨零点,利用浮标测得的这一时刻瞬时方位角、横倾角和纵倾角信息,计算得到该单个波浪周期的参考波向;然后把各个波浪周期测得的单个波浪波向进行归并统计,将360°分成16个角度区间,统计16个角度区间中每个区间内单个波浪波向出现的个数,进而得到波向出现率,最后,选择波向出现率最大的角度区间的中值作为主波向。[10][11][12][13][14]。见参考文献[10-14](唐原广,王金平.szf型波浪浮标系统[j].海洋技术学报,2008,27(2):31-33.王娟.波浪浮标传感器综合测试系统[d].中国海洋大学,2010:28-30.王金平.szf型波浪浮标数据采集系统及数据处理软件设计[d].中国海洋大学,2008:14-15.张莹.szf型波浪浮标数据采集、处理与无线传输系统[d].中国海洋大学,2010:8-9.贺成柱.海浪谱估计方法研究及szf型波浪浮标数据分析软件设计[d].中国海洋大学,2007:42-43.)

国外针对第二类工作原理的浮标,目前能够从公开文献上查到的测量波向方法与上述国内相关文献的方法基本一致,最早查到的文献是由美国国家数据浮标中心公布了一种采用横倾、纵倾传感器进行波浪波向测量的研究工作,但未公布具体的算法[15]。测量波向的方法由kennethsteele,dwwang等人首次公布,该方法首先测得浮标的方位角、横倾和纵倾,然后根据这些参数计算浮标体的北向和东向斜度,最后得到波向[16][17][18]。通过比对,国内针对第二类原理浮标的波向测量方法与kennethsteele方法类似。见参考文献[15-18](steelek,lauj,hsuyh.theoryandapplicationofcalibrationtechniquesforanndbcdirectionalwavemeasurementsbuoy[j].ieeejournalofoceanicengineering,1985,10(4):382-396.steeleke.pitch-rollbuoymeanwavedirectionsfromheaveacceleration,bowmagnetism,andstarboardmagnetism[j].oceanengineering,2003,30(17):2179-2199.steeleke,wangdw,m.d.earle,etal.buoypitchandrollcomputedusingthreeangularratesensors[j].coastalengineering,1998,35(1–2):123-139.steeleke,earlemd.directionaloceanwavespectrausingbuoyazimuth,pitch,androllderivedfrommagneticfieldcomponents[j].ieeejournalofoceanicengineering,1991,16(4):427-433.)

本发明涉及的方法适用于上述第二类工作原理的波浪方向浮标,通过传感器测得的载体原始方位、横倾和纵倾数据,对kennethsteele和唐原广团队等提出的方法进行了改进,实现了一种基于浮标测得的方位角、横倾角和纵倾角信息进行波浪主波向计算的新方法。传统算法简单选择波向出现率最大的角度区间中值作为主波向,分辨率较差,且没有充分考虑到次强出现率波向测量结果的贡献度,而改进方法充分考虑了单个波浪波向测量结果的离散性特征,通过计数比例加权矢量平均的方法计算主波向,不仅有效提高了主波向的测量精度,还提高了波向计算的分辨率。采用国产的szf型波浪浮标在海上进行了现场比测试验,结果表明本发明计算的结果与现场实测波向结果的吻合度更高,改进方法优于浮标自身计算输出结果,即唐原广团队提出的波向测量方法,提高了主波向测量的可靠性和准确性。



技术实现要素:

本发明的目的在于提供了一种基于计数比例加权平均的波浪方向浮标主波向计算方法,目的在于提高主波向测量的可靠性和准确性,更加精确地进行海洋监测。

本发明的目的是这样实现的:

一种基于计数比例加权平均的波浪方向浮标主波向计算方法,具体的实现步骤如下:

步骤1.数据预处理,读取浮标测量记录的一组连续原始数据,包括瞬时波高、方位角、纵倾角、横倾角,根据浮标测得的波高信息,利用上跨零点将该组数据切分成若干个单个波浪周期,并确定单个波浪周期的总个数n;

步骤2.计算单个波浪波向,针对每一个切分好的单个波浪周期,找到下跨零点,再利用浮标测得的该点这一时刻瞬时方位角、横倾角和纵倾角信息,计算得到该周期的单个波浪波向,重复该步骤,直到得到该组数据中所有的单个波浪波向;

步骤3.统计单个波浪波向测量结果在16个角度区间的计数比例,将360°划分成16个角度区间,统计16个角度区间中每个区间内单个波浪波向出现的个数,得到计数比例;

步骤4.计算主波向,从全部数据中选出参与计算主波向的角度区间及其计数比例,再利用计数比例作为权重进行加权矢量平均处理选出的数据,得到正北分量和正东分量,计算出主波向。

所述的步骤1具体包括:

步骤1.1.读取波浪浮标在海上进行海浪监测时测得的一组连续实验数据包括瞬时波高、方位角、纵倾角、横倾角,根据浮标测得的波高信息,采用上跨零点法,把波面上升与零点的交点,即瞬时波高数据由小于0变至大于0的点作为起点,当波浪运行到零点以下,即瞬时波高出现了小于0的情况后,将其后再次波浪上升与零点的相交点作为终点,则起点与终点之间为一个单个波浪周期;

步骤1.2,重复步骤1.1的方法,从第一个上跨零点开始,将该组实验数据切分成若干个单个波浪周期,并统计出该组数据中单个波浪周期的个数为n。

所述的步骤2具体包括:

步骤2.1.针对步骤1切分好的单个波浪周期,在一个波浪周期内,找到该周期内的下跨零点,即波面下降与零点的交点;

步骤2.2.根据得到的下跨零点,利用浮标测得的这一时刻瞬时方位角、横倾角和纵倾角信息,计算出该点瞬时东向斜度和北向斜度,再计算出这一下跨零点的瞬时波向,作为单个波浪波向;浮标体采样点瞬时东向斜度的计算公式为

其中a为方位角,p为横倾角,r为纵倾角;浮标体采样点瞬时北向斜度的计算公式为

单个波浪波向θw的计算公式为

步骤2.3.重复步骤2.1和步骤2.2,直到得到该组数据中所有的单个波浪波向θw1…θwn,其中n为切分的单个波浪周期个数。

所述的步骤3具体包括:

步骤3.1.以22.5°为间隔将360°均匀分成16个角度区间,以地理正北方位为零度,按照顺时针方向来划分;

步骤3.2.统计步骤2得到的所有单个波浪波向θw1…θwn在16个角度区间内分别出现的个数n1…n16;

步骤3.3.计算波浪波向在16个角度区间的计数比例,计数比例pi的计算公式为:

其中ni为单个波浪波向在第i个角度区间出现的个数,n为波浪周期的总个数。

所述的步骤4具体包括:

步骤4.1.根据步骤3得到的计数比例,选出参与计算主波向的数据,选取方法为,首先确定计数比例中的比例最大值k,其次选择所有计数比例不小于0.7k的数据,若选择的数据个数小于3个,则由大到小继续选择相应数据直至补足3个;将计算主波向数据的角度区间中值集合记为m,角度区间所对应的计数比例集合记为p,集合中的元素个数a不小于3;

步骤4.2.根据步骤4.1选出的角度区间中值集合m={m1,m2,…,ma}和计数比例集合p={p1,p2,…,pa},将计数比例值作为比例权重,使用加权矢量平均的方法计算得到主波向的正北分量和正东分量,正北分量δdn的计算公式为

其中pi为计数比例集合p中的第i个计数比例,mi为角度区间集合m中的第i个角度区间中值;正东分量δde的计算公式为

步骤4.3.由主波向正北分量δdn和正东分量δde,计算主波向,主波向θ计算公式为

本发明的有益效果在于:本发明充分考虑了单个波浪波向测量结果的离散性特征,通过计数比例加权矢量平均的方法计算主波向,不仅有效提高了主波向的测量精度和可靠性,还提高了波向计算的分辨率,更加精确地进行海洋监测。

附图说明

图1为上跨零点划分波浪周期方法。

图2为某一个波浪周期中的下跨零点。

图3为计数比例柱状图。

图4为本发明实施方式流程图。

具体实施方式

下面将结合附图对本发明提出的一种基于计数比例加权平均的波浪方向浮标主波向计算方法作进一步的详细说明。

本发明公开了种一种基于计数比例加权平均的波浪方向浮标主波向计算方法,目的在于提高主波向测量的可靠性和准确性,更加精确地进行海洋监测。其实施步骤如下:

步骤1,数据预处理:

读取浮标测量记录的一组连续原始数据(包含瞬时波高、方位角、纵倾角、横倾角等),根据浮标测得的波高信息,利用上跨零点将该组数据切分成若干个单个波浪周期,并确定单个波浪周期的总个数。

步骤2,计算单个波浪波向:

针对每一个切分好的单个波浪周期,先找到下跨零点,再利用浮标测得的该点这一时刻瞬时方位角、横倾角和纵倾角信息,计算得到该周期的单个波浪波向。重复该步骤,直到得到该组数据中所有的单个波浪波向。

步骤3,统计上述单个波浪波向测量结果在16个角度区间的计数比例:

先将360°划分成16个角度区间,统计16个角度区间中每个区间内单个波浪波向出现的个数,进而得到计数比例。

步骤4,计算主波向:

从步骤3获得的全部数据中选出参与计算主波向的角度区间及其计数比例,再利用计数比例作为权重进行加权矢量平均处理这些选出的数据,得到正北分量和正东分量,进而计算出主波向。

1.所述步骤1包括以下步骤:

步骤1.1,读取波浪浮标在海上进行海浪监测时测得的一组连续实验数据(包含瞬时波高、方位角、纵倾角、横倾角等),根据浮标测得的波高信息,采用上跨零点法,把波面上升与零点的交点,即瞬时波高数据由小于0变至大于0的点作为起点,当波浪运行到零点以下,即瞬时波高出现了小于0的情况后,将其后再次波浪上升与零点的相交点作为终点,则起点与终点之间为一个单个波浪周期。

步骤1.2,重复步骤1.1的方法,从第一个上跨零点开始,将该组实验数据切分成若干个单个波浪周期,并统计出该组数据中单个波浪周期的个数为n。

2.所述步骤2包括以下步骤:

步骤2.1,针对步骤1切分好的单个波浪周期,在一个波浪周期内,找到该周期内的下跨零点,即波面下降与零点的交点。

步骤2.2,对步骤2.1得到的下跨零点,利用浮标测得的该点这一时刻瞬时方位角、横倾角和纵倾角信息,先计算出该点瞬时东向斜度和北向斜度,再计算出这一下跨零点的瞬时波向,作为单个波浪波向。

浮标体采样点瞬时东向斜度的计算公式为:

式中:a为方位角,p为横倾角,r为纵倾角;

浮标体采样点瞬时北向斜度的计算公式为:

式中:a为方位角,p为横倾角,r为纵倾角;

单个波浪波向θw的计算公式为:

式中:为浮标体采样点瞬时东向斜度,为浮标体采样点瞬时北向斜度;

步骤2.3,重复步骤2.1和步骤2.2,直到得到该组数据中所有的单个波浪波向θw1…θwn。式中n为切分的单个波浪周期个数。

3.所述步骤3包括以下步骤:

步骤3.1,以22.5°为间隔将360°均匀分成16个角度区间,以地理正北方位为零度,按照顺时针方向来划分。

步骤3.2,统计步骤2得到的所有单个波浪波向θw1…θwn在16个角度区间内分别出现的个数n1…n16。

步骤3.3,计算波浪波向在16个角度区间的计数比例。

计数比例pi的计算公式为:

式中:ni为单个波浪波向在第i个角度区间出现的个数,n为波浪周期的总个数;

4.所述步骤4包括以下步骤:

步骤4.1,根据步骤3得到的计数比例,选出参与计算主波向的数据,选取方法如下:

(1)确定计数比例中的比例最大值k;

(2)首先选择所有计数比例不小于0.7k的数据;

(3)若上一步选择的数据个数小于3个,则由大到小继续选择相应数据直至补足3个。

根据上述方法得到的参与计算主波向数据为:角度区间中值集合记为m,角度区间所对应的计数比例集合记为p,集合中的元素个数a不小于3。

步骤4.2,根据步骤4.1选出的角度区间中值集合m={m1,m2,…,ma}和计数比例集合p={p1,p2,…,pa},将计数比例值作为比例权重,使用加权矢量平均的方法计算得到主波向的正北分量和正东分量。

正北分量δdn的计算公式为:

式中:pi为计数比例集合p中的第i个计数比例,mi为角度区间集合m中的第i个角度区间中值;

正东分量δde的计算公式为:

式中:pi为计数比例集合p中的第i个计数比例,mi为角度区间集合m中的第i个角度区间中值;

步骤4.3,由得到的主波向正北分量δdn和正东分量δde,计算主波向。

主波向θ计算公式为:

式中:δdn为正北分量,δde为正东分量。

本发明实施方式流程图见图4,具体可以分为以下几步,第一步为数据预处理,第二步为计算单个波浪波向,第三步为统计单个波浪波向测量结果在16个角度区间的计数比例,第四步为计算主波向。

本发明实施例所用的波浪浮标为国产的szf型波浪浮标,部署在平潭海域进行海上测波实验。工作方位为1小时定时测量方式,浮标在每天的整时进行波浪测量工作,采样间隔为0.5s,采样点为2048,浮标每1小时会连续采集1024秒的波浪数据,即一组实验数据。

上述szf型波浪浮标的主要技术参数如表一所示:

表一szf型波浪浮标的技术参数

结合附图1~4,本发明具体实施步骤为:

第一步为数据预处理。包括以下步骤:

步骤1.1,读取波浪浮标在海上进行海洋监测测得的实验数据,本例选取的一组连续实验数据观测,时间为2010年12月24日23:59,根据浮标测得的波高信息,采用上跨零点法,把波面上升与零点的交点,即瞬时波高数据由小于0变至大于0的点作为起点,当波浪运行到零点以下,即瞬时波高出现了小于0的情况后,将其后再次波浪上升与零点的相交点作为终点,则起点与终点之间为一个单个波浪周期。附图1为上跨零点划分波浪周期方法,图示中标注的圆圈所示的点为上跨零点,双箭头所示的区间为上跨零点法确定的单个波浪周期。

步骤1.2,利用上跨零点法,从第一个上跨零点开始,将该组实验数据切分成155个单个波浪周期,即统计该组数据中单个波浪周期的个数为155。

第二步为计算单个波浪波向。包括以下步骤:

步骤2.1,根据步骤1切分好的单个波浪周期,在一个波浪周期内,找到该周期内的下跨零点,即波面下降与零点的交点。附图2为某一个单个波浪周期中的下跨零点,图示中标注的圆圈所示的点为该单个波浪周期的下跨零点。

步骤2.2,对步骤2.1得到的下跨零点,利用浮标测得的该点这一时刻瞬时方位角、横倾角和纵倾角信息,先计算出该点瞬时东向斜度和北向斜度,再计算出这一下跨零点的瞬时波向,作为单个波浪波向,本例此时单个波浪波向计算结果θw=91°。

浮标体采样点瞬时东向斜度的计算公式为:

式中:a为方位角,p为横倾角,r为纵倾角;

浮标体采样点瞬时北向斜度的计算公式为:

式中:a为方位角,p为横倾角,r为纵倾角;

单个波浪波向θw的计算公式为:

式中:为浮标体采样点瞬时东向斜度,为浮标体采样点瞬时北向斜度;

步骤2.3,重复步骤2.1和步骤2.2,直到得到该组实验数据中所有的单个波浪波向。

第三步为统计波浪波向在16个角度区间的计数比例。包括以下步骤:

步骤3.1,以22.5°为间隔将360°均匀分成16个角度区间,以地理正北方位为零度,按照顺时针方向来划分。角度区间划分标准如表二所示。

表二角度区间划分标准

步骤3.2,统计步骤2得到的所有单个波浪波向θw1…θwn在16个角度区间内分别出现的个数为8,9,8,23,18,10,19,11,5,7,7,6,4,7,8,5。

步骤3.3,计算波浪波向在16个角度区间的计数比例,附图3为计数比例柱状图。

计数比例pi的计算公式为:

式中:ni为单个波浪波向在第i个角度区间出现的个数,n为波浪周期的总个数,本例中为155;

第四步为计算主波向。包括以下步骤:

步骤4.1,根据步骤3得到的计数比例,选出参与计算主波向的数据,选取方法如下:

(1)确定计数比例中的比例最大值k;

(2)首先选择所有计数比例不小于0.7k的数据;

(3)若上一步选择的数据个数小于3个,则由大到小继续选择相应数据直至补足3个。

根据上述方法得到的k=14.8%,此区间为第4角度区间,则选择的角度区间的计数比例不小于10.4%,得到两个角度区间为第5角度区间和第7角度区间,计数比例分别为11.6%和12.3%,且选取的角度区间个数为3个,已满足选取条件。则本例角度区间中值集合为m=[67.5,90,135],和角度区间所对应的计数比例集合为p=[14.8%,11.6%,12.3%]。

步骤4.2,根据步骤4.1选出的角度区间中值集合m=[67.5,90,135]和计数比例集合p=[14.8%,11.6%,12.3%],将计数比例值作为比例权重,使用加权矢量平均的方法计算得主波向的正北分量为-0.0303,计算得主波向的正东分量为0.3397。

正北分量δdn的计算公式为:

式中:pi为计数比例集合p中的第i个计数比例,mi为角度区间集合m中的第i个角度区间中值;

正东分量δde的计算公式为:

式中:pi为计数比例集合p中的第i个计数比例,mi为角度区间集合m中的第i个角度区间中值;

步骤4.3,由得到的主波向正北分量δdn=-0.0303和正东分量δde=0.3397,计算得到主波向为95°。

主波向θ计算公式为:

式中:δdn为正北分量,δde为正东分量;

在本海区的长期观测结果表明,利用wavex雷达测波仪测量的波峰峰向与具有多年观测经验的观测员人工观测结果一致性非常好,因此利用wavex实测的波峰峰向数据和本算法计算的浮标自测结果进行了比测,以检验算法性能。

1、在该时刻,wavex实测结果为103°。

2、传统的浮标测量波向方法如下:首先根据上跨零点来确定单个波浪周期,再找到每一个单个波浪周期的下跨零点,利用浮标测得的这一时刻瞬时方位角、横倾角和纵倾角信息,计算得到该单个波浪周期的参考波向;然后把各个波浪周期测得的单个波浪波向进行归并统计,将360°分成16个角度区间,统计16个角度区间中每个区间内单个波浪波向出现的个数,进而得到波向出现率,最后,选择波向出现率最大的角度区间的中值作为主波向。该算法计算的结果为67.5°。

误差分析表明,本发明反演结果与wavex实测结果的绝对误差e1=8°,浮标传统算法测量结果与wavex实测结果的绝对误差e2=35.5°。本发明反演结果与wavex实测结果基本吻合,表明本改进方法有效,且优于传统算法。

本发明反演结果与wavex实测结果的绝对误差e1计算公式为:

e1=|θ-θwavex|

式中:θ为本发明反演结果,θwavex为wavex实测结果;

浮标传统算法测量结果与wavex实测结果的绝对误差e2计算公式为:

e2=|θbuoy-θwavex|

式中:θbuoy为浮标传统算法测量结果,θwavex为wavex实测结果;

为避免偶然性,选取不同时段的多组实验数据,按照上述步骤进行实验验证。本用例选取了35组不同时段的波浪浮标实验数据进行实验验证,本发明反演结果、浮标传统算法测量结果和wavex实测结果对比如表三所示。

表三本发明反演结果,浮标传统算法测量结果和wavex实测结果对比

实验结果表明,与wavex实测结果对比,本发明反演结果的精度为9°,浮标传统算法计算结果的精度为45°。

本发明所提出的一种基于计数比例加权平均的波浪方向浮标主波向计算方法,不仅在实际测量中有效,与人工观测结果吻合度更高,且在精度上更加优于浮标传统算法测量的主波向结果,提高了波浪浮标测量主波向的稳定性和准确率,可在海洋观测设备中进行广泛的推广和应用。

当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1