研究文章|开放获取
韩少林,木下武那 “非零前进速度下船舶非线性横摇阻尼力矩的随机逆辨识“,工程数学问题那 卷。2012那 文章ID.769385那 22 页面那 2012. https://doi.org/10.1155/2012/769385
非零前进速度下船舶非线性横摇阻尼力矩的随机逆辨识
摘要
船舶横摇运动的非线性响应以横摇阻尼力矩为特征,是造船和海洋工程的重要研究方向。非线性阻尼力矩的建模和识别是将固有非线性纳入船舶设计、分析和控制的关键。文献中提出了一种用于识别一般机械系统非线性阻尼的随机非参数方法(Han and Kinoshits 2012)。该方法还被应用于零前进速度下船舶非线性阻尼力矩的识别(Han和Kinoshits 2013)。船舶在有前进速度的情况下,由于升力的作用,横摇阻尼力矩特性发生了显著的变化。本文应用随机逆方法对非零前进速度下船舶的非线性阻尼矩进行了辨识。在控制条件下进行了室内试验,验证了该方法的可行性和有效性。在试验中考虑了两种不同类型的船舶横摇运动:时变运动和频变周期运动。结果表明,该方法可以估计阻尼力矩的固有非线性,并对其进行可靠性分析。
1.介绍
船舶的运动是由船舶在海上所能经历的六自由度来定义的。其中,横摇运动是指船舶在纵轴上的旋转运动,由于较大的横摇运动可能会对船舶的安全造成严重威胁,如船舶倾覆或结构破坏,因此近年来引起了广泛的研究关注。因此,一个合适的横摇运动模型对于准确预测给定海况下船舶的横摇响应至关重要。
关于滚动运动的建模已有大量的研究[1-20.].船舶的横摇运动可以通过分析惯性矩、阻尼矩和恢复矩等力矩的分量来表征。其中,由于横摇阻尼力矩的数值难以确定,因此认为横摇阻尼力矩是最重要的分量。困难的原因是由于流体粘度的非线性和对前进速度的依赖。船舶横摇阻尼矩的识别通常采用参数法[8.-14,它预先假定为非线性的形式。然后,通过将滚动衰减试验每个周期的势能损失与作为阻尼的能量耗散联系起来,估计出规定形式的非线性系数。这种方法得到的结果相当准确,但当采用不准确的非线性模型时,由于不同类型的非线性模型产生的运动响应不同,有时会失败[15].
大量的研究集中在上述参数识别上。相比之下,很少有涉及非参数识别的研究,其中阻尼矩的非线性的先验知识是不必要的。例如,Haddara和Hinchey [16]提出了一种基于神经网络技术与标准参数识别的组合的方法,用于从自由滚衰减曲线建模非线性阻尼力矩。他们还将该方法应用于海洋试验期间测量的随机卷响应,以研究海上全尺寸船舶的卷特性[17].罗伯茨和瓦斯拉[18]提出了一种利用响应能量包络的马尔可夫模型估计阻尼矩和激励谱的随机方法。
近年来,逆方法被提出[1那2那19那20.,用于船舶阻尼矩的非参数识别。从概念上讲,这些方法是基于反形式主义的,源于对一个原始非线性运动方程的转换,用于船舶在确定性[19那20.或随机方式[1那2].基于反问题公式,从系统响应的测量中识别出零前进速度下船舶的横摇阻尼力矩[2那20.],这在意义上没有提出,即输入数据的小变化可能导致错误的逆解决方案。与确定性逆方法相比[19那20.],即随机逆方法[1那2]是稳健和可靠的,因为它还可以对已确定的结果进行可靠性分析。
在文献中[2[介绍了对零前进速度的情况的应用研究的应用。成功识别了在零前进速度下测试船的阻尼片段的固有非线性,包括其可靠的间隔作为识别可靠性的指标。然而,没有保证该方法可以应用于非零前进速度的情况,因为由于在前进速度存在下的提升效应,辊阻尼的非线性特性发生了显着变化。众所周知,船的卷筒阻尼力矩也强烈依赖于前进速度[14].
在本文中,注意力集中在随机逆方法的实际实际应用中[1]来识别以“非零”前进速度运动的船舶的非线性阻尼矩。为此,我们首先通过将非线性阻尼矩定义为一个非确定性参数,即多元随机变量,推导出一个随机逆模型。然后将这种随机逆建模方法应用于室内试验,以评估其可行性和实用性。在不同的试验条件下,考虑了瞬态运动和强迫周期运动两种不同的运动。该方法的独特特点可以概括为以下几个方面:第一,它是非参数的,即与传统参数方法不同,它不需要规定形式的非线性。其次,该方法基于随机逆模型,给出了给定噪声数据下辨识解置信度的量化方法。
本文的提纲如下。本节推导了船舶非线性阻尼的随机逆模型2.在一节中解释了实验设置3..部分4.给出了实验数据的分析结果。本节作结束语5..
2.问题公式化
2.1.船舶横摇运动控制方程
假设滚动运动一艘以前进速度行驶的船,如图所示1,由以下非线性运动微分方程控制: 在哪里是实际质量转动惯量,是附加质量惯性矩,为非线性横摇阻尼力矩,是恢复的时刻吗为滚激力矩。在上式中,上点表示对时间的微分。
横摇阻尼力矩表示为横摇角和角速度的正非线性函数。横摇阻尼力矩也受到船舶前进速度的影响,因为由于前进速度的存在而产生升力效应。恢复的时刻,表示为横摇角的反对称函数,它是由处于平衡状态的流体在重力作用下施加的静压引起的。如果横摇角的幅值足够小,则常表示为船舶位移与稳心高度与重心之间的距离的乘积,也就是说,.在本研究中,我们的注意力局限于小幅度运动。
2.2。非线性阻尼矩的逆设置
通过基于参数变分概念的数学运算,得到以下关系式[1那2那19那20.]: 在哪里和,分别满足以下方程: 的左边2.2)是由 在哪里和.
研究的目的是对给定的响应数据进行非线性阻尼矩的反辨识.获取一组测得的横摇响应数据在,对于未知非线性阻尼,给出如下系统那: 在哪里
非线性阻尼的辨识可以通过将矩阵系统(2.5)从可观察参数.该识别过程可视为逆问题[21-23,因为它的目的是从物理可观测数据的影响中找到原因。
值得注意的是,系统(2.5)通过离散化第一类积分算子(2.2).根据逆问题理论[21-23,这种系统的逆过程对数据的微小变化非常敏感。反问题的数据驱动性质使反分析更加困难,因为在反问题中,应该使用包含各种噪声源产生的随机误差的测量数据。对于这种常被称为不适定的系统,通常的逆过程会产生物理上没有意义的错误解。
2.3.随机反演
为保证求解过程稳定,未知阻尼力矩这里会被建模成一个随机变量序列吗那,在那里是否每个可能的结果都是一个任意的非空样本空间.然后,随机变量是否与直接可观测量相关根据贝叶斯公式[23]: 在哪里的可能性,先验概率密度函数,和为后验概率密度函数。
使用适当的概率模型,概率表达式(2.8)可以更清楚地加以说明。当测量误差为均值和标准差均为零的独立加性高斯随机噪声时,可能的形式: 在哪里是指欧几里得范数为测量的次数。利用两两马尔可夫随机场模型[23-25),前可以写成 的矩阵被定义为 在哪里这个点的邻居数是多少,意味着和是相邻的。因此,(2.8)可指定为 的参数和因为概率模型也是不确定性的,而且很难被先验所知。在随机反演中,这些不确定性通过扩展到层次模型中自然得到了解决[23-25]: 在哪里一对伽马分布对吗,是先验的一对逆分布吗.基于层次模型(2.13)时,可以同时确定稳定的逆解和超参数和通过马尔可夫链蒙特卡罗等数值采样技术[1那23-25].
2.4.密度模拟
为了从所构造的随机逆模型中提取阻尼矩信息,需要采用马尔可夫链蒙特卡罗等仿真技术,其目的是根据目标密度得到一组相同的独立分布样本。层次模型(2.13)通过以下混合算法可以探索给定的可观察量,该混合算法是通过在GIBBS采样器中混合Metropolis-Hastings步骤来设计的[1那23-26].(我)初始化 那,.(2)为了样本从样本从 样本从样本从均匀分布上样本从如果那其他的样本从均匀分布上样本从如果,其他的.
采样的一套实现然后可以用来近似目标密度的统计吗经过 在哪里是蒙特卡罗模拟的数量和是狄拉克函数。
综上所述,船舶在非零前进速度运动时的阻尼力矩可以通过以下步骤进行识别。(1)推导逆模型(2.5)将期望阻尼力矩与可观测参数联系起来,可观测参数是横摇响应的函数;(2)测量船舶在非零前进速度下的动力横摇响应;(3)建立随机逆模型(2.8)给出测得的横摇响应;(4)利用所设计的MCMC算法识别阻尼矩。
3.实验装置
通过对试验模型的阻尼力矩辨识,验证了该方法的可操作性和准确性。相关实验是在东京大学海洋工程(OE)盆地进行的。水池通常被称为拖曳水池,是用船模进行水动力试验的物理水池。数字2显示用于实验分析的测试模型和实验设置的概述。测试模型的方案如图所示3..表格1总结了测试模型的细节。在实验过程中,试验模型在除滚转运动外的所有自由度中进行刚性夹紧,以尽量减少其他模态运动的影响。
|
||||||||||||||||||||
(一)
(b)
数字4.显示了OE盆地的布局,长50 m,宽10 m,深5 m。盆内装有牵引车,牵引车在盆两侧的两条轨道上以0 - 2米/秒的速度运行。拖车还配备了计算机和设备来测量或控制拖曳拖车的速度。试验模型首先放置在水池的一侧,然后由马车拖到另一侧,速度从0到2米/秒不等。横摇运动由附在试验模型重心上的电位器记录下来。
该模型首先进行了测试,没有任何附件,如舭龙骨,以考虑船体阻尼特性。之后,为评估其可操作性,在试验模型的舭部转弯处,两侧分别安装两根长约1 m的舭龙骨(BK),如图所示5..BK的安装产生了完全不同的横摇特性,因为BK增加了横摇运动的水动力阻力,使船舶的横摇更少。即使在零前进速度下,BK的影响也很容易观察到,而在非零前进速度存在时,BK的影响会变得更大。
在实验应用中,考虑了两种不同类型的运动:时变的瞬态运动和频变的周期运动。首先,考虑了由初始横摇角引起的瞬态运动。首先通过静态方法对试验模型施加一个外部力矩来给出一个初始横摇角。然后消去力矩,测量滚动衰减运动。其次,考虑了周期激励引起的强迫运动。单频正弦滚转运动首先由图中力振荡装置产生的垂直力施加6..由此产生的垂直力由测压元件记录下来,并通过与力臂相乘转化为激励力矩。
(一)
(b)
应该注意的是,恢复力矩的线性近似仅对足够小的横摇角有效。如前所述,在本研究中,我们将注意力局限于小振幅运动,从而使恢复力矩可以用线性形式表示。
4.实验数据分析
4.1.瞬态运动:自由衰减滚动
作为第一个应用,考虑了由初始横摇角引起的瞬态运动。这种运动称为自由衰减滚动运动。试验中,初始滚转角为当试验模型被马车以前进的速度拖着时.然后消除静力矩,并用测量设备记录合成响应。这种自由衰减滚动运动由下列初值问题控制(2.1)自: 平分(4.1),我们可以获得 在哪里和.
记录的滚转响应如图所示7.对于没有BK和有BK的测试模型,这里定义弗劳德数(Fr)为,在那里为试验模型的速度,是重力加速度,和是船的长度。可以观察到反应测试模型与BK衰减到零的速度比没有BK测试模型。此外,减少横摇角变大的速度前进速度增加时的两种情况下没有和BK模型。这清楚地表明,阻尼是依赖于前进速度。这一事实是众所周知的,并在各种论文中进行了描述[14那27].
(a)没有bk
与BK (b)
值得注意的是,物理系数对于自由衰减滚动是不需要的。应用这种方法所需要的唯一东西就是固有角频率。试验模型的固有频率如表所示1.由于添加舱底龙骨而导致的辊自然频率的变化非常大。在本研究中,为简单起见,最大傅里叶组分的频率被认为是BK模型的自然频率。
通过对图中实验数据的应用来说明该方法7.,选择具有FR = 0.2的滚角数据的特定情况。作为识别的第一阶段,滚角数据集()转换为一组可观察量使用(2.7).如果从激励力矩开始记录横摇角,这显然是直接的对于自由衰减滚动运动是零。因此,可观测量可以用结果如图所示8..利用这个量,我们可以构造一个随机逆模型这两种情况下的实验。
(a)没有bk
与BK (b)
在说明MCMC仿真的结果之前,对逆解决方案的最小二乘估计(2.5)首先出现在图中9..结果清楚地表明(的标准逆解的困难2.5),即溶液缺乏稳定性。在解决数据驱动的逆问题时,经常会遇到符号频繁变化的不稳定性,因为在采集和分析被测信号时,无论多么小的误差都是不可避免的,这是不稳定性的主要原因。
(a)没有bk
与BK (b)
现在我们考虑随机逆模型并对其进行了MCMC仿真,得到了稳定可靠的解。为了模拟的目的,那,和0向量用于算法的初值。一个较小的值是为参数对选择的吗和.获得的模型模拟,生成了样本,最后样本用于估计的统计数据,如平均值和标准偏差.MCMC结果如图所示10.上面和下面的虚线表示95%可信区间,量化给出的测量结果的不确定性程度。与图中的结果相比,后验密度的均值估计是相当稳定的9..图中的结果10描述多元随机变量的相对可能性,其中包含未知非线性阻尼矩的信息。后验均值表示随机逆模型的最可能值.置信区间可以解释为以实测数据为条件的估计值的可靠性指标。
(a)没有bk
与BK (b)
值得注意的是,MCMC模拟需要一段时间来正确地采样目标分布。在本研究中,我们使用组件的演化来检查链条是否正常工作。数字11显示了用BK试验MCMC结果的边际分布的示踪图示例。各组分的示踪图明显处于燃烧阶段,即静止状态。这意味着后密度成功地利用所设计的MCMC算法进行了探索。
(一)
(b)
(c)
(d)
一旦概率密度函数,则非线性阻尼矩可由.识别出的非线性横摇阻尼力矩如图所示12.值得注意的是,本方法完全非参数。因此,首先基于估计函数的最可能值的一组数据揭示非线性阻尼矩的信息.然后用散点集指定阻尼矩的非线性模型。在这里,形式用于拟合已识别的数据。系数的单位为[s]-1]和[S]的.
(a)没有bk
与BK (b)
检查识别模型的准确性也很重要。为此,利用辨识出的阻尼力矩重新模拟了横摇运动。数字13给出了模拟结果与实测结果的对比。试验结果表明,试验结果与实测结果吻合较好。
(a)没有bk
与BK (b)
最后,将上述方法应用于自由衰减滚动运动的其它实验数据。结果汇总于表中2.结果清楚地说明了提高前进速度对横摇阻尼力矩的影响。此外,有BK的试验比没有BK的试验得到的阻尼矩要大得多。
|
|||||||||||||||||||||||||||||||||||||||||||||||||
根据表中确定的结果,可以自然地得出结论2,可以为非线性阻尼矩作为形式开发一个实证公式: 在哪里是零前进速度的线性贡献。的价值,可以通过拟合表中识别的结果来确定2结果证明是对于非bk模型和BK模型见图14.可以发现,在BK模型中,前进速度的影响要大得多。
4.2.周期运动:强迫滚动运动
作为第二个应用,考虑了由周期激励引起的周期强迫运动。船舶以前进速度运动时,产生了单频周期运动: 相应的卷筒令人兴奋的时刻然后用测压元件测量,如第3..值得指出的是,在将所提方法应用于强迫滚动运动时,与自由滚动衰减试验不同,首先需要知道滚动惯性矩。这些值是通过相关测试先验确定的。
对于强迫振荡试验,一般通过傅里叶分析得到横摇惯性矩的频率相关系数: 在哪里是令人兴奋的时刻的幅度,横摇角和激振力矩的相位差,和是令人兴奋的频率。所得系数在图中示出15.
(a)没有bk
与BK (b)
现在我们准备把目前的方法应用于实测数据。为了说明这一问题,本文选取了一个特定的实验数据作为受迫滚动运动的识别实例。数字16给出了试验模型在Fr = 1.0和激励频率下的横摇响应和激励力矩.数字17给出了从测量数据识别横摇阻尼力矩的分步结果。测量数据首先转换为可观测参数通过 (2.7).转换后的量,如图所示(17日),可构造随机逆模型,如(2.13).在不失一般性的前提下,我们使用相同的条件来探索所构造的MCMC算法的随机逆模型,详见本节4.1.MCMC结果显示在图中17(b).接下来,用识别的逆解决方案重新诱发辊响应。结果如图所示17(c).结果表明,模拟结果与实测结果吻合较好。这证明了识别的横摇阻尼力矩的准确性。
(一)
(b)
(a)g
(b)获得的结果
(c)再模拟侧倾响应
这里值得注意的是,对于自由衰减滚动运动,滚动衰减对滚动幅值的影响主要是由于阻尼[28,但在强迫滚动运动的情况下就不是这样了。可以预料,恢复过程中的非线性会起作用,估计的阻尼力矩并不仅仅依赖于横摇角速度。因此,与本节中描述的自由衰减滚动运动的情况不同,很难建立识别的非线性阻尼矩的解析模型4.1.
相反,我们在这里进行了一个定量分析的结果,通过应用上述程序的所有其他实验数据的强迫滚动运动。这背后的基本原理是,由于施加了单频周期运动,相关的结果可以被认为与激励频率是周期性的。在这种情况下,用频率来说明和比较结果是方便的。数据18和19显示从目前的随机识别过程中识别的解的峰值。在所有情况下,运动的幅度是相同的。应该注意的是函数包括阻尼力矩的信息。也就是figure中的结果18和19可以用阻尼力矩取决于前进速度和激励频率这一事实来解释。
(a)确定解决方案与前进速度
(b)确定了溶液与激励频率的关系
(a)确定解决方案与前进速度
(b)确定了溶液与激励频率的关系
5.结论
本文研究了一种用于识别非零前进速度船舶横摇特性的随机逆方法。滚动运动被视为一个单自由度非线性运动方程,与其他运动不耦合。在此基础上,导出了非线性阻尼矩的随机逆模型,该模型以可观测参数为横摇角测量值和激励值的函数的概率表达式表示。随机逆模型以多元随机变量的形式包含非线性阻尼贡献的信息。根据实测数据,采用设计的马尔可夫链蒙特卡罗算法识别非线性阻尼矩。
为保证该方法的适用性,将该方法应用于船舶横摇运动的瞬态运动和强迫周期运动两种不同力学现象的实验数据。在非线性系统辨识中,辨识结果的好坏取决于辨识的目的,因而很难确定辨识结果的好坏。本研究的目的是寻找非线性系统模型,可以再现测量系统的响应。由此可见,该方法能够准确地识别非零前进速度下船舶阻尼的非线性。
预设随机逆方法有以下局限性。首先,该方法基于滚动运动幅值较小的假设,使恢复力矩可以近似为线性形式。在现实中,恢复力矩也是非线性的。非线性贡献不可忽视。具有非线性恢复的系统的推广在数学上并不困难。然而,为了识别目的,新的公式需要恢复非线性的额外信息。其次,实验设置没有严格反映船舶的实际运动情况。横摇运动通常与其他运动,如摇摆和垂荡相结合。在现实中,耦合效应也应该被考虑在内。但是,为简化实验应用,本文在匀速运动时,对试验模型除横摇外的所有运动都进行了限制。
致谢
作者要感谢编辑和匿名审稿人提供的宝贵意见、建议和建设性的批评,这对论文的实质性改进非常有帮助。
参考文献
- 韩少林,“基于随机逆方法的非线性动力系统的非线性阻尼辨识,”工程数学问题文章编号574291,2012。视图:出版商的网站|谷歌学者
- 韩少林,“船舶非线性横摇阻尼力矩的随机逆建模,”,T. Kinoshita,应用海洋研究, vol. 39, pp. 11-19, 2013。视图:谷歌学者
- r . Bhattacharyya,船舶动力学,约翰瓦利和儿子,纽约,纽约,美国,1978年。
- s . n . Blagoveshchensky船舶运动理论。我和二世,多佛出版物,纽约,纽约,美国,1962年。
- J.O. de Kat和J. R. Paulling,“船舶运动的仿真和严重的海洋倾覆”sname的交易,第117卷,第127-135页,1989。视图:谷歌学者
- A. Zborowski和M. Taylan,“对共振条件下的小型船舶横摇运动稳定性储备的评估”,刊于第14艘船舶技术和研究明星研讨会1989年的诉讼程序与Sname Spring春会结合1989年4月,美国路易斯安那州新奥尔良市。视图:谷歌学者
- J.A.Witz,C. B. Ablett,以及J. H. Harrison,“具有非线性恢复时刻特征的半手会的滚动响应”应用海洋研究,卷。11,不。3,PP。153-166,1989。视图:谷歌学者
- j·f·达尔泽尔,《关于船舶减摇形式的说明》船舶研究杂志第22卷第2期第3页,178 - 185,1978。视图:谷歌学者
- J. B. Roberts,“从自由衰变数据估计非线性船舶卷起”,“船舶研究杂志,第29卷,第2期2,页127-138,1985。视图:谷歌学者
- J. B. Mathisen和W. G. Price,“船舶横摇阻尼系数的估计”,RINA事务,第127卷,第169-186页,1985。视图:谷歌学者
- D. W. Bass和M. R. Haddara,《船舶横摇阻尼的非线性模型》,国际造船进度,卷。35,不。401,pp。5-24,1988。视图:谷歌学者
- J. R. Spouge,“大幅度滚动实验的非线性分析”,国际造船进度,卷。35,不。403,pp。271-320,1988。视图:谷歌学者
- 黄文龙,“船舶大振幅横摇运动的非线性阻尼系数估计,”应用海洋研究,第十七卷,第二期第4页,217-224页,1995。视图:谷歌学者
- Y. Himeno,“船舶横摇阻尼的预测——最先进的状态”,技术代表239,美国密歇根大学造船与海洋工程系,安阿伯,密歇根州,1981。视图:谷歌学者
- M. Taylan,“船舶横摇中的非线性阻尼和恢复效应”,海洋工程第27卷第2期9,页921-932,2000。视图:出版商的网站|谷歌学者
- M. R. Haddara和M. Hinchey,“关于使用神经网络技术在自由滚衰减曲线分析中,”国际造船进度,卷。42,pp。166-178,1995。视图:谷歌学者
- M. R. Haddara和M. Wishahy,“海上两艘全尺寸船的横摇特性调查”,海洋工程,第29卷,第2期6,页651-666,2002。视图:出版商的网站|谷歌学者
- J. B. Roberts和M. Vasta,“随机波浪中非线性船舶滚动的马尔科夫模型和随机识别”,英国皇家学会哲学汇刊A,第358卷,第358号1771页,1917-1941页,2000年。视图:出版商的网站|谷歌学者|Zentralblatt数学|MathSciNet
- 张天顺,崔海生,韩淑玲,“从瞬态数据中检测非线性振动系统的非线性阻尼和恢复力的新方法”,非线性力学学报,第44卷,第5期。7,第801-808页,2009。视图:出版商的网站|谷歌学者
- “基于瞬态响应的船舶非线性横摇阻尼力矩非参数识别的数值研究”,开放海洋工程杂志,第3卷,第100-107页,2010。视图:谷歌学者
- a·基尔希反问题数学理论导论,卷。120,Springer,纽约,纽约,美国,1996。视图:出版商的网站|MathSciNet
- c . w . Groetsch数学科学中的逆问题, Informatica International, 1993。视图:MathSciNet
- J. Kaipio和E. Somersalo,统计和计算逆问题,春天,纽约,纽约,美国,2005年。视图:MathSciNet
- J. Wang和N. Zabaras,“热传导逆问题的层次贝叶斯模型”,逆问题第21卷第2期1,页183-206,2005。视图:出版商的网站|谷歌学者|Zentralblatt数学|MathSciNet
- 王建军,“多孔介质流场中污染源识别的马尔可夫随场论模型,”,N. Zabaras,国际传热与传质学报,第49卷,第49期。5-6,第939-950页,2006。视图:出版商的网站|谷歌学者
- C. Andriu,N. de Freitas,A. Doucet和M. I. Jordan,“MCMC为机器学习介绍”,机器学习,第50卷,第5期。1-2,页5-43,2003。视图:出版商的网站|谷歌学者
- 张士生,“三艘小型渔船前倾速度对横摇阻尼的影响”,海洋力学与北极工程学报,卷。116,没有。2,pp。102-108,1994。视图:谷歌学者
- G. Bulian,“用衰减时间历程的解析近似解估计非线性横摇衰减参数”,国际造船进度,第51卷,第5-32页,2004。视图:谷歌学者
版权
版权所有©2012韩少林、木下武。这是一篇发布在创意公共归因许可证,允许在任何媒介上不受限制地使用、传播和复制,但必须正确引用原作。