本发明涉及数字近景摄影测量技术领域，具体而言，涉及一种振动测量方法，尤其是一种基于结构边缘亚像素骨架提取的无标记的振动测量方法，适用于具有显著边缘特征的结构或对象的动态测量。
背景技术：
传统的模态测试大多采用加速度、速度或位移传感器来获取振动响应，一个传感器只能获得单点特定方向的响应输出，对于密集型空间采样需安装多个传感器，易出现安装布线困难、测试效率低和附加质量等问题。
薄壁结构一般是指由薄板、薄壳和细长杆件组成的结构，能以较小的重量和较少的材料承受较大的载荷，且包含显著的边缘特征，在风力发电和航空航天等诸多工程领域应用广泛。例如，飞机广泛采用的薄壁结构是指由薄型板件和加强构件组成的结构，板件有蒙皮、腹板、隔板等，加强构件有桁条和梁、肋、框和缘条等；风力机叶片也是一种薄壁结构的典型应用。该类结构件尺寸长、面积大，传统的测量方法已不能满足测试需求。
专利号为201310656022.9的发明专利申请公开了一种利用hough变换和骨架提取特征点的扩展目标跟踪方法，采用高斯平滑滤波对待处理图像进行预处理以去除噪声对后续算法的影响，采用模糊c均值聚类算法fcm(fuzzyc-meanscluster)对上述平滑后的图像进行分割，获得二值图像，然后对分割后获得的二值目标图像进行hough变换，检测出目标上直线特征明显的部分，而对于目标上线特征不明显的部分采用骨架提取获得特征点，再对骨架上的特征点进行直线拟合，获得目标轴线，最后将上述获得的轴线所在直线的交点作为最后的跟踪点，从而实现对姿态变化较大情况下的稳定跟踪。但是该方法在进行飞机机身和机翼的定位跟踪时，只将轴线所在直线的交点作为最后的跟踪点，跟踪的精度不高，同时，该方法只能进行目标点的跟踪，难以获得某些特殊边缘的振动响应，因而，仍然不能很好满足某些薄壁结构(如风力机叶片)的振动测试需求。
摄影测量技术具有全场、非接触等优点，如何应用于振动测试领域是近年来学者们的研究热点。现有的方法大多采用在物体表面布设、粘贴标志或者纹理的方式实现特征点的动态追踪。然而，在一些特殊工况或高温、高压等恶劣环境中人力难以在工件表面布置标记。
因此，如何实现了无人工合作标志情况下，仅通过被测物体本身的边缘信息来实现动态跟踪，进行振动测量，成为本领域技术人员亟需解决的技术问题。
技术实现要素：
鉴于此，本发明提出了一种基于结构边缘亚像素骨架提取的无标记振动测量方法，基于双目动态摄影立体成像技术可实现虚拟标志点的三维空间振动响应测量。
本发明提出的一种振动测量方法，包括如下步骤：
s10：采用双摄像机组成的立体摄站系统对待测对象进行动态连续拍摄；
s20：对拍摄获得的图像进行去噪处理；
s30：对去噪处理后的图像进行阈值分割，获得包含边缘特征的二值图像；
s40：对二值图像进行对象边缘的粗定位；
s50：获取对象边缘的像素骨架；
s60：确定对象边缘的像素骨架的起始点和末端点，并将起始点与末端点之间的骨架线段分割成n份，将用于分割骨架线段的分割点确定为虚拟目标点；
s70：获取连续拍摄的图像的像素骨架，并对像素骨架采用与步骤s60中相同分割方式，对相同的虚拟目标点进行动态追踪、匹配；
s80：对立体摄站系统进行系统标定，获得立体摄站系统的内外参数；
s90：计算虚拟目标点的空间振动响应。
上述的振动测量方法，优选地，在步骤s20中，采用自适应维纳滤波方法，对拍摄获得的图像进行去噪处理。
上述的振动测量方法，优选地，在步骤s30中，采用otsu双阈值法，对去噪处理后的图像进行阈值分割。
上述的振动测量方法，优选地，在步骤s40中，通过形态学方法，对二值图像进行对象边缘的粗定位。
上述的振动测量方法，优选地，步骤s50包括：
s51：获得边缘的单像素骨架；
s52：求取单像素骨架上各点的法向方向，并计算法向方向上对象边缘的光学中心，获得对象边缘的亚像素骨架。
上述的振动测量方法，优选地，所述k3m顺序迭代法包括如下步骤：
s511：提取最新的目标轮廓，并记录所有轮廓点；
s512：依次检测所有轮廓点的8像素邻域，判断是否只含有3连通像素，如果有，把只含有3连通像素的轮廓点删除，并在目标图像中删除对应点；
s513：依次检测s512中剩余轮廓点的8像素邻域，判断是否只含有3或4连通像素，如果有，把只含有3或4连通像素的轮廓点删除，并在目标图像中删除对应点；
s514：依次检测s513中剩余轮廓点的8像素邻域，判断是否只含有3或4或5连通像素，如果有，把只含有3或4或5连通像素的轮廓点删除，并在目标图像中删除对应点；
s515：依次检测s514中剩余轮廓点的8像素邻域，判断是否只含有3或4或5或6连通像素，如果有，把只含有3或4或5或6连通像素的轮廓点删除，并在目标图像中删除对应点；
s516：依次检测s515中剩余轮廓点的8像素邻域，判断是否只含有3或4或5或6或7连通像素，如果有，把只含有3或4或5或6或7连通像素的轮廓点删除，并在目标图像中删除对应点；如果在这步骤中仍有轮廓点被删除，则进入第7步，否则迭代终止，得到最终骨架；
s517：依次检测目标图像中像素骨架的8像素邻域，判断是否只含有2或3或4或5或6或7连通像素，如果有，把只含有2或3或4或5或6或7连通像素的轮廓点删除，得到最终的骨架。
上述的振动测量方法，优选地，步骤s52包括：
s521：针对像素骨架上的轮廓点f(x,y)，按sobel算子求取梯度向量g[fx,fy]；
s522：则点f(x,y)处的法线方向夹角
s523：获得法线方向后，以f(x,y)为基点，针对该法线方向上的每一个轮廓点f(x',y')，其中x'＝x+l·nx，,y'＝y+l·ny，以临近的三个骨架点pi-1,pi-2,f(x',y')，求对应曲率，取最小曲率时的f(x',y')为新找到的骨架点pi；
s524：连接所有的pi即可得到亚像素骨架。
上述的振动测量方法，优选地，在步骤s51中，通过k3m顺序迭代法获得边缘的单像素骨架；在步骤s52中：采用能量加权平均法计算法向方向上对象边缘的光学中心。
上述的振动测量方法，优选地，在步骤s60中，确定对象边缘的像素骨架的起始点和末端点，并计算起始点和末端点之间的骨架线段的像素长度，采用等分方式将起始点与末端点之间的骨架线段分割成n份。
上述的振动测量方法，优选地，在步骤s90中，根据双目立体摄影成像原理求解计算虚拟目标点的空间振动响应，空间振动响应包振动加速度、振动速度和振动位移。
本发明提出的振动测量方法，采用双摄像机组成的立体摄站系统对待测对象进行动态连续拍摄；对拍摄获得的图像进行去噪处理；对去噪处理后的图像进行阈值分割，获得包含边缘特征的二值图像；对二值图像进行对象边缘的粗定位；获取对象边缘的像素骨架；确定对象边缘的像素骨架的起始点和末端点，并将起始点与末端点之间的骨架线段分割成n份，将用于分割骨架线段的分割点确定为虚拟目标点；获取连续拍摄的图像的像素骨架，并对像素骨架采用与前述步骤中相同分割方式，对相同的虚拟目标点进行动态追踪、匹配；对立体摄站系统进行系统标定，获得立体摄站系统的内外参数；计算虚拟目标点的空间振动响应。本发明的振动测量方法，无需提前在被测物体上粘贴或布置合作标志即可实现振动测量，不会造成结构表面的污染和损伤，特别适用于高温、高压等特殊工况中人力难以布置标记的场合，可以节约大量的人力物力。
在进一步的方案中，对振动中的运动物体拍摄后，采用自适应维纳滤波方法消除图像退化的影响，与传统的维纳滤波相比，本方法可自适应地选取信噪比参数，具有更好的适应能力。而且为了克服滤波后边缘位置的模糊化问题，本方法单独对图像进行了sobel算子提取边缘图像，并叠加到滤波后图像上，获得的图像更为平滑，边缘更清晰。
在进一步的方案中，采用两步法获得高精度骨架，即首先采用k3m顺序迭代算法提取单像素骨架；再在骨架各点的梯度方向上求取新的骨架点，使得新骨架线曲率变化最小。该方法综合考虑几何中心和光学中心，可获得连续光滑的亚像素骨架线。
在进一步的方案中，对提取的骨架线等分，并可采用多帧同名点匹配算法，获得振动物体上若干点的连续振动位移曲线，实现了无标记振动监测，避免了粘贴标记的繁琐和识别标记的困难等问题。
附图说明
图1为本发明具体实施例提出的一种振动测量方法的流程图；
图2为具体实施例中在一种叶片边缘结构拍摄原图；
图3为阈值分割后的图像；
图4为腐蚀运算后的图像；
图5为单像素骨架提取的结果；
图6为骨架等分求交点示意图；
图7为实施例中9个虚拟目标点的空间运动轨迹图；
图8为实施例中某个虚拟目标点在y方向上的振动响应图；
图9为双目立体摄影成像原理示意图。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互组合。
需要说明的是，本发明中的“和/或”，指既可以同时具备，也可以选择其中之一，如方案a和/或方案b，包括方案a、方案b、方案a且方案b此三种情况。
如图1-8所示，本发明提出的一种振动测量方法，针对具有显著边缘特征的结构进行振动测量，无需提前布设合作标志，该方法基于边缘的亚像素骨架提取和交比不变特性构造虚拟目标点，通过对相同的虚拟目标点进行追踪匹配、摄影成像实现待测对象的动态测量。包括如下步骤：
s10：采用双摄像机组成的立体摄站系统对待测对象进行动态连续拍摄；
待测对象是需要进行振动测量的结构或部件，具体可以根据需要进行选择，如图1中所示，以某型风力机叶片作为待测对象，叶片根部固定在测试平板上，叶尖端自由，叶片水平放置。根据仿真分析可知叶片第三阶频率不超过40hz，故选用120fps的拍摄频率，用两个工业摄像机组成的立体摄站系统进行摄影测量。图像采集时，打开摄像机，施加外力载荷使叶片处于运动状态，保证左右摄像机同时触发拍照并将所采集的连续图像帧同步传输到pc端上。
s20：对拍摄获得的图像进行去噪处理；
如图2所示，对步骤s10中拍摄的图像，由于采集的图像本身存在噪声，因此需要对读入的图像进行去噪处理，去噪处理可以通过多种方式实现，如高斯平滑滤波，从而得到平滑的图像，在一种优选的方式中，采用自适应维纳滤波方法进行滤波去噪处理，其原理如下：
含噪声图像模型为：g(x,y)＝h[f(x,y)]+n(x,y)
其中h为退化函数，f(x,y)为原图像，n(x,y)为噪声，g(x,y)为含噪声图像；
要从g(x,y)恢复原始图像f(x,y)，可在频率域按如下维纳滤波公式进行：
其中g(u,v)是退化图像的傅立叶变换，h(u,v)是退化函数的傅立叶变换，pn(u,v)是噪声功率谱，pf(u,v)是原始图像的功率谱。信号功率谱pf(u,v)常难以估计，通常可设为可变参数
按以上公式获得了f(u,v)之后，再进行逆傅立叶变换即可得到滤波后的图像f1(x,y)，虽然维纳滤波方法去噪效果好，但有时参数k选择不当会造成边缘信息丢失。因此，可根据实际情况，采用小波模极大值法(或sobel算子等其他方法)提取图像边缘，得到边缘信息图像f2(x,y)，则新的滤波去噪后且含边缘信息的图像为f(x,y)＝f1(x,y)+f2(x,y)。
上述方案中，对振动中的运动物体拍摄后，采用自适应维纳滤波方法消除图像退化的影响，与传统的维纳滤波相比，本方法可自适应地选取信噪比参数，具有更好的适应能力。而且为了克服滤波后边缘位置的模糊化问题，本方法单独对图像进行了sobel算子提取边缘图像，并叠加到滤波后图像上，获得的图像更为平滑，边缘更清晰。
s30：对去噪处理后的图像进行阈值分割，获得包含边缘特征的二值图像；
该步骤可以使图像背景与目标分割开，阈值分割可以通过多种方式实现，比如模糊c均值聚类算法(fcm,fuzzyc-meanscluster),在一种优选的方式中，如图3所示，优选采用otsu双阈值方法进行图像的阈值分割，具体方法如下：
c0＝{0～t1}c1＝{t1+1～t2}c2＝{t2+1～255}
设c0出现的概率为ω0，其灰度平均值为m0；c1出现的概率为ω1，其灰度平均值为m1。c2出现的概率为ω2，其灰度平均值为m2，灰度值为i的像素数为ni，总像素数为n，则：
各类灰度值出现的概率pi为：
pi＝ni/n
整幅图像的灰度平均值m为：
进一步得到：
s40：对二值图像进行对象边缘的粗定位；
该步骤在于去除区域中不连续的局部噪声干扰，只保留目标区域，优选通过形态学方法，分割出来的目标区域(二值图像)进行形态学运算，如图4所示。需要说明的是，形态学方法可以包括腐蚀、膨胀、开运算或闭运算等方法中的一种或多种，具体采用何种形态学方法，可以根据具体情况进行选择，比如在粗线条时，通过形态学腐蚀将线条变细，形成粗定位。
s50：获取对象边缘的像素骨架；
该步骤中，可以采用现有方式获得对象边缘的像素骨架，如中轴变换(medialaxistransform，mat)、迭代替换算法等等，具体可以参考现有技术，像素骨架也可以包括多种形式，如单像素骨架、亚像素骨架等等，优选地，步骤s50包括：
s51：获得边缘的单像素骨架；
s52：求取单像素骨架上各点的法向方向，并计算法向方向上对象边缘的光学中心，获得对象边缘的亚像素骨架。
进一步优选地，如图5所示，在步骤s51中，通过k3m顺序迭代法获得边缘的单像素骨架；具体地，针对目标图像，确定目标外围轮廓，然后从外轮廓开始，不断腐蚀目标图像边界(这是一次迭代过程)，直至腐蚀到不能再腐蚀为止。步骤如下：
s511：提取最新的目标轮廓(初始为原目标轮廓)，并记录所有轮廓点；
s512：依次检测所有轮廓点的8像素邻域，判断是否只含有3连通像素，如果有，把只含有3连通像素的轮廓点删除，并在目标图像中删除(腐蚀掉)对应点(即只含有3连通像素的轮廓点)；
s513：依次检测s512中剩余轮廓点的8像素邻域，判断是否只含有3或4连通像素，如果有，把只含有3或4连通像素的轮廓点删除，并在目标图像中删除(腐蚀掉)对应点(即只含有3或4连通像素的轮廓点)；
s514：依次检测s513中剩余轮廓点的8像素邻域，判断是否只含有3或4或5连通像素，如果有，把只含有3或4或5连通像素的轮廓点删除，并在目标图像中删除(腐蚀掉)对应点(即只含有3或4或5连通像素的轮廓点)；
s515：依次检测s514中剩余轮廓点的8像素邻域，判断是否只含有3或4或5或6连通像素，如果有，把只含有3或4或5或6连通像素的轮廓点删除，并在目标图像中删除(腐蚀掉)对应点(即只含有3或4或5或6连通像素的轮廓点)；
s516：依次检测s515中剩余轮廓点的8像素邻域，判断是否只含有3或4或5或6或7连通像素，如果有，把只含有3或4或5或6或7连通像素的轮廓点删除，并在目标图像中删除(腐蚀掉)对应点(只含有3或4或5或6或7连通像素的轮廓点)；
如果在这步骤中仍有轮廓点(即像素)被删除(即腐蚀)，则进入第7步，否则迭代终止，得到最终骨架(即单像素骨架)；
s517：依次检测目标图像中像素骨架的8像素邻域，判断是否只含有2或3或4或5或6或7连通像素，如果有，把只含有2或3或4或5或6或7连通像素的轮廓点删除(腐蚀)，从而得到最终的骨架，即单像素骨架。
同时，在步骤s52中：采用能量加权平均法计算法向方向上对象边缘的光学中心。具体地，步骤如下：
s521：针对骨架上的轮廓点f(x,y)，按sobel算子求取梯度向量g[fx,fy]，sobel算子是在prewitt算子的基础上改进的，在中心系数上使用一个权值2，相比较prewitt算子，sobel模板能够较好的抑制(平滑)噪声；
s522：则点f(x,y)处的法线方向夹角
s523：获得法线方向后，以f(x,y)为基点，针对该法线方向上的每一个轮廓点f(x',y')，其中x'＝x+l·nx，,y'＝y+l·ny(其中l为步长参数，取决于图像细化程度)，以临近的三个骨架点pi-1,pi-2,f(x',y')，求对应曲率，取最小曲率时的f(x',y')为新找到的骨架点pi，由于l·nx非整数，因此我们最后得到的pi具有亚像素精度；
s524：连接所有的pi即可得到具有亚像素精度的骨架，即亚像素骨架。
该方式中，通过采用两步法获得高精度骨架，即首先采用k3m顺序迭代算法提取单像素骨架；再在骨架各点的梯度方向上求取新的骨架点，使得新骨架线曲率变化最小。该方法综合考虑几何中心和光学中心，可获得连续光滑的亚像素骨架线。
s60：确定对象边缘的像素骨架的起始点和末端点，并将起始点与末端点之间的骨架线段分割成n份，将用于分割骨架线段的分割点确定为虚拟目标点；
该步骤中，可以通过差分、随机等方式设置分割点，n可以为5个、6个、7个、8个、9个、10个、15个、20个等。优选地，在步骤s60中，确定对象边缘的像素骨架的起始点和末端点，并计算起始点和末端点之间的骨架线段的像素长度，采用等分方式将起始点与末端点之间的骨架线段分割成n份。如图6所示，如将骨架线等分为8份，则相应生成9个交点(分割点，分割点可以包括起始点和末端点)，交点对应被测对象边缘上8等分的相应点，此方法构造了无标记场合下的9个虚拟目标点。
该步骤中，对提取的骨架线等分，并可采用多帧同名点匹配算法，获得振动物体上若干点的连续振动位移曲线，实现了无标记振动监测，避免了粘贴标记的繁琐和识别标记的困难等问题。
s70：获取连续拍摄的图像的像素骨架，并对像素骨架采用与步骤s60中相同分割方式，对相同的虚拟目标点进行动态追踪、匹配；
需要说明的是，该步骤中，可以采用步骤s20-s50相同的方法获取连续拍摄的图像的像素骨架，并采用步骤s60的方法针对像素骨架进行分割，从而获得虚拟目标点，并对虚拟目标点进行动态追踪、匹配，通过连续拍摄图像，实现连续追踪和匹配。该步骤中，追踪主要指虚拟目标点的位置追踪，匹配包括两层含义，第一，对于同一摄像机连续拍摄的图像，均需要采用相同的分割方法进行分割，从而得到相同数目的虚拟目标点(在上述的分割方式中，如虚拟目标点1、虚拟目标点2……虚拟目标点9)，对于不同时间拍摄的图像(即连续拍摄的不同图像)，相同的虚拟目标点需要匹配，如第1个摄像机在第1时刻拍摄图像的虚拟目标点3与第1个摄像机在其它时刻(第2-n，n为最后时刻)拍摄图像的虚拟目标点3均需要匹配；第二，对于两个不同的摄像机同一时刻拍摄的图像，相同的虚拟目标点需要匹配，比如第1时刻、第1个摄像机拍摄图像的虚拟目标点3与第1时刻、第2个摄像机拍摄图像的虚拟目标点3，需要匹配。
s80：对立体摄站系统进行系统标定，获得立体摄站系统的内外参数；
立体摄站系统包括两个工业相机，内外参数具体指每个相机的内参和外参，内参包括焦距、主点坐标、畸变系数、比例因子等，外参是两个摄像机坐标系的空间变换关系(r、t矩阵)。
s90：计算虚拟目标点的空间振动响应。
计算方式可以采用现有方法，振动响应可以根据需要确定，振动响应具体可以是振动加速度、振动速度和振动位移等几种变量形式。优选地，根据双目立体摄影成像原理求解计算虚拟目标点的空间振动响应，如图7示出了9个虚拟目标点的空间运动轨迹图，图8示出了某个虚拟目标点在y方向上的振动响应图，当然根据附图7，还可以得出该虚拟目标点在x方向和z方向的振动响应图，也可以得出其它虚拟目标点在x、y和z方向的振动响应图，根据振动响应图可以得出振动的位移变化量，根据位移与时间的关系，可以计算得到速度、加速度响应。双目立体摄影成像原理可以参考附图9，该系统中的一对工业相机包括左相机和右相机，设左相机像空间坐标系o-xyz与标记点的坐标系重合，像平面坐标系为ol-ulvl，有效焦距为fl，右相机像空间坐标系or-xryrzr，像平面坐标系为or-urvr，有效焦距为fr；设标记点p在o-xyz中的坐标为(x,y,z)，其在左照片中对应的像点pl在o-xyz中的坐标为(x,y,-fl)，标记点p在右照片中对应的像点pr在or-xryrzr中的坐标为(xr,yr,-fr)；
对于左照片，有o、pl和p三点共线有：
即：
对于右照片，由or、pr和p三点共线有：
式中，(x',y',z')为点p在or-xryrzr中的坐标，即：
o-xyz坐标系与or-xryrzr坐标系之间的相互位置关系通过旋转和平移实现转换，设
将(5)式代入(4)式得：
根据相机焦距参数、待测空间点在左、右照片中的图像坐标、旋转矩阵r和平移矩阵t，联合(2)式与(6)式计算标记点的三维坐标，通过该方法可以得到虚拟目标点的坐标，从而得到虚拟目标点的位移、速度、加速度等。
本发明的振动测量方法，无需提前在被测物体上粘贴或布置合作标志即可实现振动测量，不会造成结构表面的污染和损伤，特别适用于高温、高压等特殊工况中人力难以布置标记的场合，可以节约大量的人力物力。与现有技术相比的有益效果在于：
(1)该方法在提取边缘的单像素骨架基础上结合sobel梯度算子和能量加权平均法进一步获得边缘的亚像素骨架，为后续振动测量计算提供了更精确的数据。
(2)对振动中的运动物体拍摄后，采用自适应维纳滤波方法消除图像退化的影响，与传统的维纳滤波相比，本方法可自适应地选取信噪比参数，具有更好的适应能力。而且为了克服滤波后边缘位置的模糊化问题，本方法单独对图像进行了sobel算子提取边缘图像，并叠加到滤波后图像上，获得的图像更为平滑，边缘更清晰。
(3)采用两步法获得高精度骨架，即首先采用k3m顺序迭代算法提取单像素骨架；再在骨架各点的梯度方向上求取新的骨架点，使得新骨架线曲率变化最小。该方法综合考虑几何中心和光学中心，可获得连续光滑的亚像素骨架线。
(4)对提取的骨架线等分，并采用多帧同名点匹配算法，获得振动物体上若干点的连续振动位移曲线，实现了无标记振动监测，避免了粘贴标记的繁琐和识别标记的困难等问题。
以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。