一种残膜回收机防缠绕挑膜装置的制 一种秧草收获机用电力驱动行走机构

海面测量系统、海面测量方法以及存储介质

2022-05-27 02:24:03 来源:中国专利 TAG:

海面测量系统、海面测量方法以及存储介质
1.本技术是申请日为2020年01月10日、申请号为202010027415.3、发明名称为“海面测量系统、海面测量方法以及存储介质”的中国发明专利申请的分案申请。
技术领域
2.本发明涉及一种基于通过多个摄像机拍摄到的图像来测量海啸、巨浪、满潮、急潮等海面异常的海面测量系统、海面测量方法以及海面测量程序。


背景技术:

3.现有的海啸的测量(观测)方法能够大致分类为已到达陆地的海啸的观测方法和到达陆地前的海啸的观测方法两种。
4.(1)已到达陆地的海啸的观测方法
5.在已到达陆地的海啸的高度的观测中,除了使用气象局在各地设置的验潮站以外,还使用海啸观测仪、巨大海啸观测仪等海啸观测设备。气象局使用验潮站、海啸观测仪始终观测潮位,将观测到的数据作为潮汐观测资料(速报值)和潮汐观测资料来公开。但是,这些观测系统均被设置在海岸处,因此虽然能够观测已到达陆地的海啸,但是无法观测到达前的远方(几十公里(km)至超过50公里(km))的海啸,从而无法估计到达时刻。
6.(2)到达陆地前的海啸的观测方法
7.在到达陆地前的海啸的观测中,使用海底地震及海啸观测网、gps波浪计以及最近提出的海啸雷达。
8.海底地震及海啸观测网是作为日本海沟海底地震海啸观测网配备业务而由国立研究开发法人防灾科学技术研究所推进并配备的大规模的实时地震及海啸观测网络。在海底地震及海啸观测网中,通过海底光缆连接将地震仪与海啸计形成为一体而得到的观测装置,将该观测装置设置于日本东部的太平洋海底,来实时地且24小时连续地获取观测数据。截至2018年1月,在日本全国沿海的150处设置了观测装置,线缆总长约为5,700km。期待海底地震及海啸观测网直接探测海沟型地震、紧接其后的海啸,并通过迅速且高精度的信息传达来对灾害的减轻、避难行动等防灾对策做出贡献。
9.但是,海底系统、陆地基站的设置需要庞大的工程,且需要自治机构、渔业相关者同意。另外,设备设施的维持、管理也需要巨大的费用。因此,存在难以在全国的海域普及的一面。
10.gps波浪计系统是使用gps卫星测量漂浮于近海域的浮标(gps波浪计)的上下变动来直接观测波浪、潮汐等海面变动的海洋观测设备。截至2017年1月,在日本近海域设置有18台gps波浪计。
11.gps波浪计设置在港湾设备中以获得所需要的近海域的波浪信息,但是由于还能够观测在地震发生时由海啸引起的海面的上下波动,因此非常期待还被有效利用于海啸预防措施。gps波浪计通常被设置于近海域20km前后的地方,因此能够观测几十公里处的海啸并估计到达时刻。
12.但是,gps波浪计是昂贵的设备,在2005年设置之初,一台超过3亿日元,现在也是以亿为单位。另外,gps波浪计只能设置在特定的地点,难以进行广范围的观测。
13.海啸雷达是使用雷达测量技术测量远方的海面来监视海啸的发生的系统。在这些方法中,从通过雷达观测到的海面的流速中提取海啸成分,来估计海啸的波浪高度和到达时刻。
14.但是,在这些方法中,难以从海面的流速中提取海啸成分,特别是难以提前测量长周期的海啸。
15.现有技术文献
16.专利文献
17.专利文献1:日本特开2013-40898号公报
18.专利文献2:日本特开2018-4529号公报
19.专利文献3:日本特开2014-160044号公报
20.专利文献4:日本特开2004-191268号公报
21.专利文献5:日本特开2017-150917号公报
22.专利文献6:日本特开2016-85206号公报


技术实现要素:

23.发明要解决的问题
24.另外,从30年前开始研究图像测量技术,随着数字摄像机技术及计算机制造技术等硬件技术的进步、算法及人工智能等软件技术的发展,近年来图像测量技术的实用化正在发展。当前,图像测量技术被实际使用于印刷品质的管理、产品的品质管理、珠宝及艺术品的鉴定等品质管理领域。另外,还被广泛地实际使用于脸部识别、物体识别的领域。
25.本发明的目的在于提供一种应用了图像测量技术的对海啸等海面异常进行测量的海面测量系统、海面测量方法以及海面测量程序。
26.用于解决问题的方案
27.本发明的海面测量系统包括:摄影单元,其通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;波浪提取单元,其从通过多个摄像机同时拍摄到的各个图像中提取波浪;波浪对应单元,其针对由波浪提取单元提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;三维信息计算单元,其通过对由波浪对应单元对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及海面异常有无判定单元,其基于由三维信息计算单元计算出的海面的高度和波浪的高度,来判定是否发生海面异常,其中,所述波浪提取单元将由所述多个摄像机拍摄到的图像分割为几个邻接的长方形或正方形的小区域,针对包括小区域的左上角、右上角、左下角、右下角这四个角的局部区域,求出用于将图像进行二值化的阈值,基于关注点所存在的小区域的左上角、右上角、左下角、右下角这四个角的小区域阈值,乘以取决于到四个角的绝对距离的权重系数,来求出用于将关注点进行二值化的阈值。
28.另外,本发明的海面测量系统包括:摄影单元,其通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;波浪提取单元,其从通过多个摄像机同时拍摄到的各个图像中提取波浪;波浪对应单元,其针对由波浪提取单元提取出的各波浪,从通过多
个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;三维信息计算单元,其通过对由波浪对应单元对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及海面异常有无判定单元,其基于由三维信息计算单元计算出的海面的高度和波浪的高度,来判定是否发生海面异常,其中,所述波浪对应单元进行以下处理:通过极线几何和二维图像的横纵坐标关系的分析,根据由所述多个摄像机拍摄到的多张图像来求出对应候选区域;生成以在各图像的所述对应候选区域中提取出的各波浪的图像的重心坐标、像素的数量、周长、圆形度以及在各方向上的长度为特征的二维特征向量;针对各图像的同一区域内的所有波浪,基于极线约束找出候选波浪,并求出该候选波浪的二维特征向量的相似度,按照相似度从高到低的顺序制作各波浪的对应波浪的候选波浪列表;以及基于关注波浪和所述候选波浪列表中的相似度高的波浪的重心的图像坐标,求出关注波浪的三维世界坐标(x
关注
,y
关注
,z
关注
)和候选波浪的三维世界坐标(x
候选
,y
候选
,z
候选
),通过判定所述关注波浪的三维世界坐标与所述候选波浪的三维世界坐标之差是否小于事先决定的阈值,来决定对应波浪的候选波浪是否为对应波浪。
29.另外,本发明的海面测量方法包括以下步骤:通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;从通过多个摄像机同时拍摄到的各个图像中提取波浪;针对提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;通过对被对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及基于计算出的海面的高度和波浪的高度,来判定是否发生海面异常,其中,在提取波浪的所述步骤中,将由所述多个摄像机拍摄到的图像分割为几个邻接的长方形或正方形的小区域,针对包括小区域的左上角、右上角、左下角、右下角这四个角的局部区域,求出用于将图像进行二值化的阈值,基于关注点所存在的小区域的左上角、右上角、左下角、右下角这四个角的小区域阈值,乘以取决于到四个角的绝对距离的权重系数,来求出用于将关注点进行二值化的阈值。
30.另外,本发明的海面测量方法包括以下步骤:通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;从通过多个摄像机同时拍摄到的各个图像中提取波浪;针对提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;通过对被对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及基于计算出的海面的高度和波浪的高度,来判定是否发生海面异常,其中,在进行对应的所述步骤中进行以下处理:通过极线几何和二维图像的横纵坐标关系的分析,根据由所述多个摄像机拍摄到的多张图像来求出对应候选区域;生成以在各图像的所述对应候选区域中提取出的各波浪的图像的重心坐标、像素的数量、周长、圆形度以及在各方向上的长度为特征的二维特征向量;针对各图像的同一区域内的所有波浪,基于极线约束找出候选波浪,并求出该候选波浪的二维特征向量的相似度,按照相似度从高到低的顺序制作各波浪的对应波浪的候选波浪列表;以及基于关注波浪和所述候选波浪列表中的相似度高的波浪的重心的图像坐标,求出关注波浪的三维世界坐标(x
关注
,y
关注
,z
关注
)和候选波浪的三维世界坐标(x
候选
,y
候选
,z
候选
),通过判定所述关注波浪的三维世界坐标与所述候选波浪的三维世界坐标之差是否小于事先决定的阈值,来决定对应波浪的候选波浪是否为对应波浪。
31.另外,本发明的存储介质存储有海面测量程序,该海面测量程序使计算机执行上
述的海面测量方法。
32.根据这些发明,通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像,从通过多个摄像机同时拍摄到的各个图像中提取波浪,针对提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应,通过对被对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度,基于计算出的海面的高度和波浪的高度来判定是否发生海面异常。
33.另外,本发明的海面测量程序使计算机作为以下单元发挥功能:摄影单元,其通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;波浪提取单元,其从通过多个摄像机同时拍摄到的各个图像中提取波浪;波浪对应单元,其针对由波浪提取单元提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;三维信息计算单元,其通过对由波浪对应单元对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及海面异常有无判定单元,其基于由三维信息计算单元计算出的海面的高度和波浪的高度,来判定是否发生海面异常。执行本发明的海面测量程序的计算机与上述本发明的海面测量装置同样地发挥功能。
34.另外,本发明的存储介质存储有海面测量程序,该海面测量程序使计算机作为以下单元发挥功能:摄影单元,其通过设置于不同地点处的多个摄像机同时拍摄包括同一海面区域的图像;波浪提取单元,其从通过多个摄像机同时拍摄到的各个图像中提取波浪;波浪对应单元,其针对由波浪提取单元提取出的各波浪,从通过多个摄像机同时拍摄到的各个图像中找出相同的波浪并进行对应;三维信息计算单元,其通过对由波浪对应单元对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及海面异常有无判定单元,其基于由三维信息计算单元计算出的海面的高度和波浪的高度,来判定是否发生海面异常。
35.发明的效果
36.根据本发明,能够基于通过多个摄像机拍摄到的图像,来实时地判定是否发生海啸、巨浪、满潮、急潮等海面异常,在作为海面异常而发生了海啸的情况下,能够实时地计算出海啸的规模和到达时刻、或者自动地向必要的地方和必要的人发送关于海啸等海面异常的警报,能够有助于针对海面异常的避难、灾害的减轻。
附图说明
37.图1是本发明的实施方式中的海面测量系统的硬件结构图。
38.图2是示出图1的海面测量系统的摄像机组的设置状态的说明图。
39.图3是示出基于立体观察的三维图像测量的坐标关系的说明图。
40.图4是示出图1的海面测量系统的结构的框图。
41.图5是示出由图1的海面测量系统进行的海面测量的流程的流程图。
42.图6是大气模型的影像图。
43.图7是由雨影响减轻单元进行的处理的流程图。
44.图8是波浪提取的影像图。
45.图9是示出使用两张图像进行的波浪的对应处理的流程的流程图。
46.图10是波浪的特征向量的影像图。
47.图11是使用左右两台摄像机的摄影图像计算海面的高度的影像图。
48.图12是示出使用左右两台摄像机的摄影图像计算海面的高度的流程的流程图。
49.图13是示出圆柱状目标物的一例的图。
50.图14是示出由远距离摄像机系统校准单元进行的校准的流程的流程图。
51.图15是示出是否发生海啸的判定的流程的流程图。
52.图16是各时刻的海面的高度的三维图像测量的影像图。
53.图17是摄像机的视线的机械控制的影像图。
54.图18是通过两台摄像机的图像测量控制进行的精密控制的影像图。
55.图19是示出具有两台摄像机的情况下的通过图像测量控制进行的精密控制的流程的流程图。
56.附图标记说明
57.1:服务器;2a、2b:客户端计算机;3a、3b:摄像机组;4a、4b:摄像机;5a、5b:摄像机固定调节部;10:摄影单元;11:波浪提取单元;12:波浪对应单元;13:三维信息计算单元;14:海面异常有无判定单元;15:远距离摄像机系统校准单元;16:海啸行进速度计算单元;17:海啸规模推测单元;18:海啸到达时刻推测单元;19:海啸信息发送单元;20:存储单元;21:雾影响减轻单元;22:雨影响减轻单元;23:两级式控制单元。
具体实施方式
58.下面,使用附图说明本发明的实施方式中的海面测量系统。图1是本发明的实施方式中的海面测量系统的硬件结构图,图2是示出图1的海面测量系统的摄像机组的设置状态的说明图,图3是示出基于立体观察的三维图像测量的坐标关系的说明图,图4是示出图1的海面测量系统的结构的框图。
59.如图1所示,本实施方式中的海面测量系统包括作为计算机的服务器1、分别与服务器1连接的客户端计算机2a、2b以及分别与客户端计算机2a、2b连接的第一摄像机组3a及第二摄像机组3b。第一摄像机组3a和第二摄像机组3b各自包括摄像机4a、4b和旋转台等摄像机固定调节部5a、5b,该摄像机4a、4b包括摄像机主体和镜头,该摄像机固定调节部5a、5b将摄像机4a、4b以能够调节姿势的方式固定。
60.如图2所示,第一摄像机组3a和第二摄像机组3b设置在接近海洋的海岸的陆地上的基线长度b且设置高度h的位置中的不同的地点c1、c2处。此外,严格地说,如图3所示,将分别设置在地点c1、c2处的第一摄像机组3a和第二摄像机组3b各自的摄像机镜头中心间的距离设为基线长度b。另外,在下面的说明中,将基线长度b延伸的方向设为x方向,将高度方向设为y方向,将与x方向及y方向正交的方向设为z方向。
61.如图4所示,服务器1包括:摄影单元10,其通过设置于不同地点处的多个摄像机4a、4b同时拍摄包括同一海面区域的图像;波浪提取单元11,其从通过多个摄像机4a、4b同时拍摄到的各个图像中提取波浪;波浪对应单元12,其针对由波浪提取单元11提取出的各波浪,从通过多个摄像机4a、4b同时拍摄到的各个图像中找出相同的波浪并进行对应;三维信息计算单元13,其通过对由波浪对应单元12对应起来的图像中的各波浪进行三维分析,来计算海面的高度和波浪的高度;以及海面异常有无判定单元14,其基于由三维信息计算单元13计算出的海面的高度和波浪的高度,来判定是否发生海啸等海面异常。
62.另外,本实施方式中的海面测量系统还包括:远距离摄像机系统校准单元15,其用于提高海面的三维信息计算精度;海啸行进速度计算单元16,其在作为海面异常而测量出发生海啸的情况下,计算海啸的行进速度;海啸规模推测单元17,其推测海啸的规模;海啸到达时刻推测单元18,其推测海啸的到达时刻;海啸信息发送单元19;以及存储单元20。并且,在后面记述详细内容,但是本实施方式的海面测量系统还包括雾影响减轻单元21、雨影响减轻单元22以及两级式控制单元23。服务器1通过执行海面测量程序,来作为上述各单元10~23发挥功能。
63.图5是示出由图1的海面测量系统进行的海啸测量的流程的流程图。
64.为了实现作为海面异常的海啸的图像测量,首先将两台以上的摄像机4a、4b设置于不同的地点处(s101)。
65.接着,通过多个摄像机4a、4b同时拍摄包括同一海面区域的图像(s102)。另外,在夜晚、雨、雾等恶劣环境的摄影中使用红外线摄像机、雾影响减轻单元21以及雨影响减轻单元22,在后面记述详细内容(s201)。
66.之后,使用波浪提取单元11从各摄像机4a、4b拍摄到的图像中分别提取海啸测量所需要的波浪(s103)。
67.在针对各图像提取出波浪之后,使用波浪对应单元12进行波浪的对应,来进行从各摄像机4a、4b的图像中提取出的波浪的对应,并求出相同波浪的组(s104)。
68.针对进行对应得到的各波浪的组,使用海面的三维信息计算单元13计算海面的三维信息。首先,计算所检测出的各波浪的三维世界坐标(x,y,z),基于该三维世界坐标(x,y,z)计算海面的三维信息、即海面的高度和波浪的高度(s105)。
69.此外,为了确保三维信息计算精度,使用远距离摄像机系统校准单元15进行摄像机系统的校准(s202)。
70.基于所计算出的海面的高度和波浪的高度,使用是否发生海啸的判定单元判定是否发生海啸(s106、s107)。
71.在判定为发生海啸的情况下,使用海啸行进速度计算单元16、海啸规模推测单元17以及海啸到达时刻推测单元18估计海啸的规模(s108),并估计到达时刻(s109)。
72.最后,通过海啸信息发送单元19向必要的地方和必要的人发送所测量出的海啸信息(s110)。
73.接着,对各单元进行详细说明。
74.[摄影单元10]
[0075]
摄影单元10通过设置于不同地点处的多个摄像机4a、4b同时拍摄包括同一海面区域的图像。如上所述,包括摄像机4a的第一摄像机组3a和包括摄像机4b的第二摄像机组3b如图2所示那样设置在接近海洋的海岸的陆地上的不同地点c1、c2处,获取半径5km~20km左右的范围内的远距离海面的图像。为了实现该半径5km~20km的远方海面的三维图像测量,在本实施方式中,使用基于立体观察的三维图像测量技术。各摄像机4a、4b间的距离即摄像机4a、4b间的基线长度b设为20米以上,各摄像机4a、4b的设置场所的海拔即摄像机4a、4b的设置高度h设为20米以上。
[0076]
为了使由海啸等引起的对第一摄像机组3a和第二摄像机组3b的损害为最小限度、且尽可能地抑制摄影时的干扰、且实现尽可能远的地方的照片拍摄,期望将第一摄像机组
3a和第二摄像机组3b设置在海岸处,且第一摄像机组3a和第二摄像机组3b越接近海洋越好。此外,在本实施方式中,使用第一摄像机组3a和第二摄像机组3b两组,也能够设为三组以上。
[0077]
期望的是,摄像机4a、4b使用具有8比特以上的灵敏度的高灵敏度摄像机,以能够提高测量精度并且在恶劣环境中也能够拍摄测量所需要的图像。另外,期望使用大口径的具有远摄功能的镜头或具有远摄功能的变焦镜头以拍摄远距离海面。
[0078]
客户端计算机2a、2b具有调节摄像机4a、4b的参数来对摄影进行控制的功能、取得由摄像机4a、4b拍摄到的图像的功能。服务器1具有与客户端计算机2a、2b通信的通信功能,并具有以下功能:对客户端计算机2a、2b进行控制,向客户端计算机2a、2b发送照片拍摄所需要的命令,并从客户端计算机2a、2b获取摄影图像。
[0079]
此外,在本实施方式的海面测量系统中,为了拍摄不同的海面,需要调节各摄像机4a、4b的视线。因此,摄像机固定调节部5a、5b能够进行平摇方向(水平方向)和俯仰方向(垂直方向)上的旋转。通过调节摄像机固定调节部5a、5b的在平摇方向和俯仰方向上的角度,来将设置于摄像机固定调节部5a、5b的摄像机4a、4b的摄影方向角度沿上下和左右调节,从而能够使多个摄像机4a、4b拍摄同一海面。
[0080]
根据客户端计算机2a、2b的命令来控制平摇方向和俯仰方向上的旋转角度。平摇方向和俯仰方向上的旋转角的具体的值根据测量海域的位置而不同,将该值从服务器1发送到对各摄像机组3a、3b进行控制的客户端计算机2a、2b。另外,通过根据电子外部同步信号或来自服务器1的软件同步命令,来使多个摄像机4a、4b同时拍摄,由此实现摄像机4a、4b的同步摄影。由客户端计算机2a、2b进行摄像机4a、4b的摄像机主体的曝光时间等参数控制、镜头的变焦控制、聚焦等控制。首先,由与各摄像机4a、4b连接的客户端计算机2a、2b取得由各摄像机4a、4b拍摄到的图像,然后根据服务器1的命令将该图像发送到服务器1,并在服务器1中进行图像处理。
[0081]
此外,在测量海啸等海面异常时需要24小时拍摄海洋。因此,在本实施方式的海面测量系统中,作为摄像机4a、4b,白天主要使用可见光摄像机,晚上主要使用远红外线摄像机等红外线摄像机。
[0082]
[雾影响减轻单元21]
[0083]
另外,服务器1具备用于减轻雾等的影响的雾影响减轻单元21。雾影响减轻单元21基于大气模型,根据太阳光的强度成分以及由摄像机4a、4b拍摄到的海洋的图像的纵向长度、即与图像的y坐标线性近似的海面的大气的透射率分布,减轻雾等对由所述摄像机4a、4b拍摄到的图像的不良影响。
[0084]
在图6中示出大气模型的影像图。在大气环境中存在大气辉光,其强度用a表示。进入到摄像机的光量包括直接进入的大气辉光成分和来自对象物的反射光成分这两个成分。大气模型能够通过式(1)来表示。
[0085]
[式1]
[0086]
i=j
·
t a
·
(1-t)
ꢀꢀꢀꢀ
(1)
[0087]
其中,i为进入到摄像机的光量即受到雾等环境影响的图像,a为大气辉光的强度成分,j为对象物的反射成分,t为从对象物到摄像机之间的大气透射率。j为没有雾等环境影响的情况下的图像,是在由雾影响减轻单元21进行雾等的去除处理时所要求取的目标成
分。
[0088]
透射率t能够通过式(2)来表示。
[0089]
[式2]
[0090]
t=e-βd
ꢀꢀꢀꢀꢀꢀ
(2)
[0091]
其中,β为大气扩散系数,d为从对象物到摄像机的距离。
[0092]
根据式(1),能够得到式(3)。
[0093]
[式3]
[0094][0095]
即,只要能够估计大气辉光的强度成分a和大气的透射率t,则能够通过雾影响减轻单元21去除雾的影响,来获得只具有单纯的反射成分j的图像。
[0096]
在参考文献1(robby t.tan,visibility in bad weather from a single image,ieee conference on computer vision and pattern recognition,pp.1-8,2008)中,tan氏提出了基于局部区域对比度最大化的方法。提出方法是指,在局部区域内t为常数且j<a的前提下,能够基于局部区域对比度最大化来估计该区域的t。该方法仅使用对比度,因此有时会过饱和。
[0097]
在参考文献2(r.fattal.single image dehazing[j].acm siggraph

08,pp.1-9,2008)中,提出了基于独立成分分析的雾去除方法。在局部区域内,由于对象物的反射成分j的模式与大气透射率在统计上是独立的,因此能够通过独立成分分析法来估计反射率。该方法的独立性的设定是基于图像的颜色信息的设定,因此针对颜色信息变化贫乏的海洋的照片而言,适应性低。
[0098]
在本实施方式中,大气辉光强度成分a的估计使用以往的方法,因此省略详细的说明。关于大气透射率t,提出如下一种计算方法:不进行繁杂的计算,而如下述的式(4)那样使大气透射率t以线性关系近似于拍摄到的海洋的图像的纵向长度、即图像的y坐标。
[0099]
[式4]
[0100]
t=t0 ky
ꢀꢀꢀꢀꢀꢀꢀ
(4)
[0101]
其中,t0为大气透射率的初始值,k为比例系数,y为由摄像机拍摄到的海洋的图像的y坐标。由此,能够大幅地缩短计算时间,并能够迅速地计算海啸信息。
[0102]
[雨影响减轻单元22]
[0103]
另外,服务器1具备用于减轻雨的影响的雨影响减轻单元22。图7是由雨影响减轻单元22进行的处理的流程图。如图7所示,雨影响减轻单元22使用连续拍摄到的多张图像,求出像素的颜色强度的最大值图像和最小值图像,应用将最小值图像作为引导的引导滤波器(日文:
ガイドフィルター
)来对最大值图像进行处理。
[0104]
首先,以固定的时间间隔连续拍摄五张图像(s301)。
[0105]
接着,针对五张输入图像in(n=1、2、3、4、5)的各像素,通过式(5)求出像素的颜色强度的最大值,通过式(6)求出像素的颜色强度的最小值。使用各像素的最小值生成最小值图像i
min
,使用各像素的最大值生成最大值图像i
max
(s302)。
[0106]
[式5]
[0107]imax
(i,j)=max(i1(i,j),...,in(i,j))
ꢀꢀꢀ
(5)
[0108]
[式6]
[0109]imin
(i,j)=min(i1(i,j),...,in(i,j))
ꢀꢀ
(6)
[0110]
其中,(i,j)为像素坐标,n=5。
[0111]
接着,针对关注像素k(i,j),在以k(i,j)为中心的小区域ωk内,将最小值图像i
min
设为引导图像,来如下述那样求出中间图像if(s303)。
[0112]
[式7]
[0113]
if(i,j)=aki
min
(i,j) bk,(i,j)∈ωkꢀꢀꢀꢀ
(7)
[0114]
[式8]
[0115][0116]
[式9]
[0117][0118]
其中,μk和σk分别是小区域ωk内的最小值图像i
max
的像素的颜色强度的平均及方差,|ω|为小区域ωk内的像素的数量,为小区域ωk内的输入图像的颜色强度的平均值。
[0119]
最后,基于输入图像in和中间图像if,如下述那样应用引导滤波器来对最大值图像进行处理,由此去除雨的影响,求出波浪的样态被强化的输出图像(s304)。
[0120]
[式10]
[0121]iout
=(i
n-if)a ifꢀꢀꢀꢀꢀꢀ
(10)
[0122]
其中,a为调整系数。
[0123]
[波浪提取单元11]
[0124]
波浪提取单元11从通过多个摄像机4a、4b同时拍摄到的各个图像中提取波浪。从海域的摄影图像中提取波浪、特别是从设为本实施方式中的海面测量系统的目标的5km~20km的远距离摄影图像中提取波浪是难以应用一般的图案提取法的。作为其原因,能够列举波浪形状的不确定性、由多个摄像机4a、4b拍摄到的多张图像中的波浪的特征的不同、全天候摄影所伴有的摄影环境的变化等。另外,要求检测的高速性。
[0125]
图8是波浪提取的影像图。如图8所示,首先,将拍摄到的图像分割为几个邻接的长方形或正方形的小区域。针对各小区域的左上角、右上角、左下角、右下角这四个角,在以角为中心的周围选择面积与小区域的面积相同的局部区域,求出用于对该局部区域内的图像进行二值化的阈值(小区域阈值)t。
[0126]
如式(11)那样,基于关注点p(i,j)所存在的小区域的左上角、右上角、左下角、右下角这四个角的小区域阈值t,乘以取决于到四个角的绝对距离的权重系数k1~k4来计算用于图像坐标(i,j)即关注点p(i,j)的二值化的阈值t(i,j)。
[0127]
[式11]
[0128]
t(i,j)=k1t(lui,luj) k2t(rui,ruj) k3t(ldi,ldj) k4t(rdi,rdj)
ꢀꢀꢀ
(11)
[0129]
[式12]
[0130][0131]
[式13]
[0132][0133]
[式14]
[0134][0135]
[式15]
[0136][0137]
[式16]
[0138][0139]
即,针对每个关注点,根据关注点与四个角之间的距离,使用不同的权重系数进行二值化。由此,关注小区域内的各像素的二值化结果能够平滑地与邻接的上下左右的各小区域内的像素的二值化结果衔接。
[0140]
如上述那样,将通过基于用于对关注点进行二值化的阈值动态地变化的动态阈值的二值化方法(动态阈值计算算法)提取出的成分设为图像中的波浪。
[0141]
[波浪对应单元12]
[0142]
波浪对应单元12针对由波浪提取单元11提取出的各波浪,从通过多个摄像机4a、4b同时拍摄到的各个图像中找出相同的波浪并进行对应。图9是示出使用两张图像进行的波浪的对应处理的流程的流程图。图像为三张以上的情况也能够使用相似的方法。
[0143]
首先,通过极线几何和二维图像的横纵坐标关系的分析,从通过设置于不同地点处的多个摄像机拍摄到的多张图像中求出对应波浪可能存在的区域、即对应候选区域(s401)。
[0144]
接着,对各图像的对应候选区域应用波浪的提取单元,从左图像提取m个波浪,从右图像提取n个波浪(s402)。另外,如图10所示,生成以所提取出的各波浪的图像的重心坐标、像素的数量、周长、圆形度、各方向上的长度等参数为该波浪的特征的二维特征向量。作为一例,能够如下述那样定义二维特征向量v。
[0145]
[式17]
[0146]
v={g,s,l,c,t}
ꢀꢀꢀꢀ
(17)
[0147]
其中,g为波浪的重心坐标,用g(x0,y0)表示。s为表示波浪的大小的参数,为所提取出的波浪中包括的像素的数量。l为波浪的周长,c为波浪的圆形度。t为表示波浪的纵向和横向的长度的参数,用t(h1,h2,w1,w2)表示。
[0148]
[式18]
[0149]
h1=y
0-y1ꢀꢀꢀꢀꢀ
(18)
[0150]
[式19]
[0151]
h2=y
2-y0ꢀꢀꢀꢀ
(19)
[0152]
[式20]
[0153]
w1=x
0-x1ꢀꢀꢀꢀꢀ
(20)
[0154]
[式21]
[0155]
w2=x
2-x0ꢀꢀꢀꢀꢀꢀ
(21)
[0156]
此外,p1(x1,y1)为包围波浪的长方形的左上角的像素坐标,p2(x2,y2)为右下角的像素坐标。
[0157]
接着,制作左图像的波浪的特征向量vl和右图像的波浪的特征向量vr(s403)。
[0158]
针对从左图像提取出的m个波浪,按照下述的式子求出第i个(i=1,2,

,m)波浪的特征参数vli,并生成左图像的波浪的向量vl。
[0159]
[式22]
[0160]
vl
t
={vl1,vl2,......,vl
m-1
,vlm}
ꢀꢀꢀꢀꢀ
(22)
[0161]
针对从右图像中提取出的n个波浪,按照下述的式子求出第j个(j=1,2,

,n)波浪的特征参数vrj,并生成右图像的波浪的向量vr。
[0162]
[式23]
[0163]
vr
t
={vr1,vr2,

,vr
n-1
,vrn}
ꢀꢀꢀꢀ
(23)
[0164]
接着,针对各图像的同一区域内的所有波浪,基于极线约束找出候选波浪,并求出其二维特征向量的相似度,按照相似度从高到低的顺序制作各波浪的对应波浪的候选波浪列表。具体地说,针对从左图像提取出的第i个波浪(i=1,2,

,m),从右图像提取处于极线附近的k个波浪,并将该k个波浪设为右图像中的与左图像的第i个波浪对应的对应波浪的k个候选波浪(s404、s405)。
[0165]
针对k个候选波浪,计算二维特征向量的相似度,按相似度从高到低的顺序,制作与从左图像提取出的第i个波浪对应的对应波浪的候选波浪列表(s406)。以将相似度最高的波浪设为第一个候选波浪、将相似度第二高的波浪设为第二个候选波浪的方式决定候选波浪的位次。
[0166]
基于关注波浪和候选波浪列表中的高位次波浪、即相似度高的波浪的重心的图像坐标,通过三维图像测量的方法求出关注波浪的三维世界坐标(x
关注
,y
关注
,z
关注
)和该候选波浪的三维世界坐标(x
候选
,y
候选
,z
候选
)(s407),判定关注波浪的三维世界坐标(x
关注
,y
关注
,z
关注
)与该候选波浪的三维世界坐标(x
候选
,y
候选
,z
候选
)之差是否小于事先决定的阈值(s408)。
[0167]
如果关注波浪的三维世界坐标(x
关注
,y
关注
,z
关注
)与该候选波浪的三维世界坐标(x
候选
,y
候选
,z
候选
)之差小于事先决定的阈值,则该波浪设为关注波浪的对应波浪(s409),否则验证下一个波浪(s410、s411)。
[0168]
[三维信息计算单元13]
[0169]
三维信息计算单元13通过对由波浪对应单元12对应起来的图像的各波浪进行三维分析,来计算海面的高度和波浪的高度。
[0170]
如果海面有波浪,则波浪的形状、尺寸随时变化。因此,在基于海面处的波浪的三维测量结果求出海面的高度时,难以决定将哪个波浪的哪个部分的高度设为海面的高度。
[0171]
另一方面,在海啸测量中不需要求出海面的特定点的高度,只要计算出某个区域的海面的高度来判定是否发生海啸即可。
[0172]
因此,在本实施方式的海面测量系统中,根据摄影图像计算某个小区域的海面的高度。图11是使用左右两台摄像机的摄影图像计算海面的高度的影像图,图12是示出使用
左右两台摄像机的摄影图像计算海面的高度的流程的流程图。
[0173]
如图12所示,首先,通过左右的摄像机4a、4b连续拍摄同一区域,获取多张连续时间序列图像(s501)。例如以1/30秒一张的摄影速度拍摄5秒钟,通过左右的摄像机各拍摄150张图像,获取共300张图像。
[0174]
接着,从各图像提取波浪,并进行对应(s502)。
[0175]
之后,针对对应起来的所有波浪,求出其重心、中心、上部、下部、左端、右端等特征点的三维世界坐标(x,y,z)(s503)。对于摄像机而言,x表示水平方向,y表示垂直方向,z表示深度方向。
[0176]
基于计算出的各波浪的三维世界坐标(x,y,z),将波浪分类为几个小区域(s504)。例如图11所示,将波浪分类为三个小区域。
[0177]
针对特定的小区域,如式(24)那样,针对所有时间序列图像的所有被对应起来的波浪,求出波浪的重心、除此以外的特征点的三维世界坐标的纵向即垂直于海面的方向上的y坐标值的平均值ya,将平均值ya乘以系数,来计算该小区域的海面的高度hs(s505)。
[0178]
[式24]
[0179]hs
=ksya h
s0
ꢀꢀꢀꢀꢀꢀꢀ
(24)
[0180]
其中,ks为用于计算的比例系数,h
s0
为用于调整计算海面的高度时的三维世界坐标的系数。
[0181]
最后,针对特定的小区域,如式(25)那样,求出所有时间序列图像的所有被对应起来的波浪的所述y坐标值的方差值v(y),将该方差值v(y)乘以系数,来计算该小区域的波浪的高度hv(s506)。
[0182]
[式25]
[0183]hv
=kvv(y) h
v0
ꢀꢀꢀꢀꢀꢀ
(25)
[0184]
其中,kv为用于计算的比例系数,h
v0
为用于调整计算波浪的高度时的世界坐标的系数。
[0185]
[远距离摄像机系统校准单元15]
[0186]
关于摄像机系统的校准,求出摄像机4a、4b的图像传感器的像素配置及镜头的焦距等内部参数、以及两台摄像机4a、4b间的距离即基线长度、摄像机4a、4b的位置及姿势即旋转角度等外部参数。摄像机系统的校准方法被提出很多种,当前特别是对张氏方法(zhang.z,a flexible new technique for camera calibration,ieee transactions on pattern analysis and machine intelligence.vol.22,no.11,2000,pp.1330-1334.)进行评价,还将该方法实际安装在opencv中。
[0187]
但是,需要进行远距离的照片拍摄以进行海啸测量,在远距离三维图像测量系统的校准中仍然存在需要大尺寸的校准用对象物、需要从不同视角拍摄对象物等问题点。
[0188]
因此,本实施方式中的远距离摄像机系统校准单元15设为融合了使用板状目标物的内部参数校准方法和使用圆柱状目标物的外部参数校准方法的、使用不同目标物的内外分离两级式校准。
[0189]
在本实施方式的远距离摄像机系统校准单元15中,关于外部参数校准用的目标物,使用两根以上的圆柱状的目标物或其它具有固定宽度的细长的目标物。图13是圆柱状目标物的一例。
[0190]
选择圆柱状的物体作为校准的对象物的理由是因为,如图13的(a)所示,即使从不同的视点观看,从图像测量所使用的测量点到圆柱状目标物的中心线的距离也不会变化。例如,当从视线1观看目标物时,在图像上观测到测量点1,从测量点1到圆柱的中心线的距离为r。即使从视线2观测,也与从视线1观测时相同,从观测点2到圆柱的中心线的距离为r。即,在圆柱状目标物中,不会由于观测点的变化而产生校准误差。
[0191]
此外,考虑是远距离校准,将圆柱状等细长的目标物的长度设定为1米以上。另外,为了易于识别目标物,设为目标物的表面具有白色和红色、或其它容易清楚地区分的两种颜色以上的颜色。
[0192]
将图13的(b)所示的特征点候选1和特征点候选3那样的上述不同颜色的边界线上的点、或者特征点候选2那样的处于颜色的中央的点、或者重心设为校准用的特征点。
[0193]
另外,将如图13的(c)所示的l1、l2、l3那样的上述目标物的不同颜色的区域的尺寸作为已知参数来在校准时使用。
[0194]
图14是示出由远距离摄像机系统校准单元15进行的校准的流程的流程图。
[0195]
首先,通过多个摄像机4a、4b来拍摄校准目标物的照片(s601)。
[0196]
接着,从各个图像中提取校准所需要的多个特征点(s602)。
[0197]
之后,进行所提取出的特征点的对应(s603)。
[0198]
最后,使用被对应起来的特征点,并使用上述的l1、l2、l3等已知参数,来计算摄像机的参数(s604)。
[0199]
[海面异常有无判定单元14]
[0200]
海面异常有无判定单元14基于由三维信息计算单元13计算出的海面的高度和波浪的高度,来判定是否发生海啸。图15是示出是否发生海啸的判定的流程的流程图。
[0201]
如图15所示,首先,为了判别是否发生海啸,将基于三维图像测量的技术计算出的几个特定区域的海面的高度与没有发生海啸的平常时的潮位进行比较(s701)。根据比较的结果,将这些特定区域分类为发生海啸的候选区域、可能发生海啸的区域、或没有发生海啸的区域这三个区域。
[0202]
在测量出的多个特定区域内的海面的高度与没有发生海啸的平常时的潮位之差均大于规定的阈值的情况下,将这些特定区域分类为发生海啸的候选区域。在混合存在上述的多个特定区域内的海面的高度与没有发生海啸的平常时的潮位之差大于规定的阈值的区域和该差小于规定的阈值的区域的情况下,将这些特定区域分类为可能发生海啸的区域。在上述的多个特定区域内的海面的高度与没有发生海啸的平常时的潮位之差都小于规定的阈值的情况下,将这些特定区域分类为没有发生海啸的区域。
[0203]
针对发生海啸的候选区域或可能发生海啸的区域,在该区域的周边进一步追加多个小区域作为测量区域(s702),将这些测量区域的海面的高度的测量结果与这些测量区域的平常潮位进行比较(s703)。基于这些测量区域的海面的高度与平常时的潮位的比较结果,使用噬菌体推测理论(日文:
ファージ
推測理論)来判断是否发生了海啸(s704),并将该特定区域判定为发生海啸的区域或没有发生海啸的区域。此外,关于除海啸以外的巨浪、满潮、急潮等海面异常,也能够通过同样的过程进行判定。
[0204]
[海啸行进速度计算单元16、海啸规模推测单元17、海啸到达时刻推测单元18]
[0205]
海啸行进速度计算单元16使用基于以往的海啸行进速度的计算方法计算出的海
啸的行进速度,对基于三维图像测量计算出的海啸行进速度进行修正,来计算更准确的海啸行进速度。
[0206]
图16是各时刻的海面的高度的三维图像测量的影像图。
[0207]
在作为海面异常而测量出发生海啸的情况下,随时测量海啸发生区域及其周边区域内的海面的高度,并随时检测海面的高度最大的区域。另外,通过摄像机4a、4b的姿势控制,来追踪与摄像机4a、4b最近的、海面的高度最大的区域、即海啸锋面,并获取到海啸锋面的距离、即三维测量结果的z值。
[0208]
如图16所示,通过时刻1、时刻2等各时刻的海面的高度的三维图像测量,获取海面的高度变化的四维信息(x,y,z,t(测量时刻)),来计算海啸锋面的海面高度的高度变化、行进方向、以及行进速度即基于图像测量计算出的海啸行进速度vm。
[0209]
另外,通过公式(26),基于以往的计算方法计算海啸的行进速度vc。
[0210]
[式26]
[0211][0212]
其中,h为测量区域的海洋的水深,g为重力加速度。
[0213]
使用基于以往的计算方法计算出的海啸的行进速度vc,对通过图像测量计算出的海啸的行进速度vm进行校正,来计算更准确的海啸行进速度v。
[0214]
[式27]
[0215]
v=k
mvm
(1-km)vcꢀꢀꢀꢀ
(27)
[0216]
其中,km为用于海啸速度校正的系数。
[0217]
海啸规模推测单元17基于海啸的高度来评价海啸的规模。例如,海啸规模推测单元17根据海啸的高度评价为-1、0、1、2、3、4六个等级。
[0218]
海啸到达时刻推测单元18根据通过三维图像测量而测量出的海啸锋面的位置及海啸的高度以及海啸行进速度,来计算海啸到达海岸的时刻。
[0219]
[存储单元20]
[0220]
首先,观察者制作在何时测量哪个海域那样的测量计划。服务器1按照该测量计划决定各摄像机的姿势,计算各摄像机的旋转角度并发送到各客户端计算机2a、2b。存储单元20保存从各客户端计算机2a、2b获取到的各摄像机4a、4b的摄影图像。服务器1通过对存储于存储单元20的这些图像进行处理,来计算并输出是否发生海啸等海面异常。
[0221]
[两级式控制单元23]
[0222]
并且,在本实施方式的海面测量系统中,摄影单元10具有融合了基于角度反馈的机械控制与基于图像摄影及图像处理的图像测量控制的、利用机械及图像的两级式控制单元23。
[0223]
在两级式控制单元23中,首先,使用低分辨率的角度传感器进行基于角度反馈的机械控制,来进行各摄像机4a、4b的视线的粗略控制以使多个摄像机4a、4b朝向大致相同的海域。作为粗略控制的目标,粗略地调节各摄像机4a、4b的视线以使各摄像机4a、4b均能够拍摄到测量目标海面处的波浪的一部分。
[0224]
在粗略控制之后,通过各摄像机4a、4b进行照片拍摄,对由各摄像机4a、4b拍摄到的图像进行反馈,来通过图像测量控制进行各摄像机4a、4b的视线的精密控制。精密控制的
目标是:即使存在机械控制的误差、摄像机的振动等,也能够通过各摄像机4a、4b最大限度地拍摄测量目标海面处的波浪。
[0225]
图17是摄像机的视线的机械控制的影像图。机械控制为第一级的粗略控制,是确保两台摄像机4a、4b的视线朝向同一海域的控制。关于粗略控制,不要求视线控制的精度,只要能够确保两台以上的摄像机4a、4b朝向同一海域即可,即,只要能够在由两台摄像机4a、4b拍摄到的图像中检测出几个相同的波浪即可。因此,能够在粗略控制中利用0.1
°
左右的低分辨率的角度传感器。
[0226]
图18是通过两台摄像机的图像测量控制进行的精密控制的影像图。摄像机的台数为三台以上的情况也能够应用同样的控制方式。图像测量控制为第二级的精密控制,通过图像反馈来对摄像机4a、4b的角度、图像进行微调,确保由两台摄像机4a、4b拍摄到同一目标。
[0227]
为了通过左右两台摄像机4a、4b拍摄同一目标,需要向左右的摄像机4a、4b发送各自的目标视线角度。将多个摄像机中的一个设为主摄像机,将另一个设为副摄像机,而在本实施方式中,将左摄像机4a设为主摄像机,将右摄像机4b设为副摄像机,根据左摄像机4a即主摄像机的视线角度来自动地计算右摄像机4b即副摄像机的视线角度。
[0228]
图像测量控制的关键点为:反馈由左右的摄像机4a、4b各自拍摄到的图像,与作为目标的摄影图像进行比较,基于它们的差,来进一步调节摄像机固定调节部5a、5b的平摇和俯仰的角度。由此,能够校正由摄像机4a、4b的振动等产生的误差即如图18所示的干扰那样的误差。
[0229]
图19是示出具有两台摄像机的情况下的通过图像测量控制进行的精密控制的流程的流程图。首先,在服务器1中,根据摄影目标区域计算用于调节摄像机4a、4b的姿势的视线角度、即平摇和俯仰方向上的旋转角度。之后,根据摄像机4a、4b间的空间关系和摄影目标区域,使用三角测量的理论计算右摄像机4b的视线角度、即平摇和俯仰方向上的旋转角度,并发送到对各个摄像机4a、4b进行控制的客户端计算机2a、2b(s801)。
[0230]
与左摄像机4a连接的客户端计算机2a将从服务器1接收到的左摄像机4a的平摇和俯仰的旋转角度设为视线角度的目标值,将从角度传感器获取的左摄像机4a的实际的视线角度与该视线角度的目标值进行比较,以使它们的差为最小的方式调节摄像机固定调节部5a的平摇和俯仰的角度,来进行机械控制(s802)。对右摄像机4b也同样地进行机械控制(s803)。
[0231]
之后,通过左摄像机4a和右摄像机4b分别进行照片拍摄,从左右的摄影图像中分别提取波浪(s804、s805)。
[0232]
接着,根据从左右的图像中提取出的波浪,求出在左右两张图像中都存在的区域即共通区域,提取处于共通区域的中央位置的一个波浪,并将该一个波浪设为摄像机姿势调整的目标波浪(s806)。
[0233]
期望目标波浪在左右两张图像中都位于图像的中央位置。目标波浪为左右摄像机4a、4b的姿势控制的目标、即摄影目标。针对左摄像机4a判定目标波浪是否进入到了图像的中央周边处的固定范围以内(s807),如果进入到范围内,则结束左摄像机4a的姿势调节,在没有进入到范围内的情况下,进行上述的机械控制,来调节平摇和俯仰的角度(s809)。
[0234]
接着,在机械控制后进行照片拍摄,从摄影图像中提取目标波浪,并且调查是否需
要图像调节(s811、s813)。如果目标波浪进入到了图像的中央的固定区域,则不需要图像调节,否则将图像沿水平方向和垂直方向在固定范围以内进行平移,来进行图像调节(s815)。
[0235]
通过上述的图像调节,调查目标波浪是否被调节到了左图像的中央,如果位于图像的中央区域则结束调节,否则返回到s809的机械控制,重复进行控制(s809~s815)。即,无论怎么进行图像的平移而目标波浪都不在图像的中心附近的位置的情况意味着图像测量控制能够控制的范围窄且无法达成控制目标。该情况、即图像测量控制所不能涵盖的部分通过机械控制来实现。对于右摄像机4b,也同样地调节摄像机4b的姿势。
[0236]
最后,调查左右两台摄像机4a、4b的调整是否均已结束,如果均已结束,则结束摄像机的控制,否则等待直到结束为止(s819、s820)。
[0237]
即,在上述图像测量控制中,具有由获取将摄像机4a、4b以能够调节姿势的方式固定的摄像机固定调节部5a、5b的平摇和俯仰的角度的角度传感器进行的角度反馈、以及从摄像机4a、4b的照片拍摄中提取波浪来确定摄像机4a、4b的姿势的图像反馈,将角度反馈设为内反馈,将图像反馈设为外反馈,图像反馈具有围绕角度反馈的双反馈结构,首先进行机械控制,机械控制不足的部分、即机械控制的误差通过图像测量控制来补充,由此图像测量控制所不能涵盖的部分通过机械控制来实现。
[0238]
产业上的可利用性
[0239]
本发明的海面测量系统、海面测量方法以及海面测量程序也能够应用于是否发生海啸的测量、在发生了海啸的情况下对海啸的高度的测量、海啸的行进速度的测量、海啸的到达时刻的推测等。另外,不只是海啸,也能够应用于由作为海面异常的台风产生的波浪的高度、位置、到达时刻的测量。除此之外,也能够应用于巨浪、满潮、急潮等海面异常的测量、海岸的监视。
再多了解一些

本文用于企业家、创业者技术爱好者查询,结果仅供参考。

发表评论 共有条评论
用户名: 密码:
验证码: 匿名发表

相关文献