CN111854621B - 机载分布式pos用光纤光栅传感器数据拟合方法和装置 - Google Patents

机载分布式pos用光纤光栅传感器数据拟合方法和装置 Download PDF

Info

Publication number
CN111854621B
CN111854621B CN202010505955.8A CN202010505955A CN111854621B CN 111854621 B CN111854621 B CN 111854621B CN 202010505955 A CN202010505955 A CN 202010505955A CN 111854621 B CN111854621 B CN 111854621B
Authority
CN
China
Prior art keywords
vibration
wavelength variation
wing
data
fitting
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010505955.8A
Other languages
English (en)
Other versions
CN111854621A (zh
Inventor
宫晓琳
孙一弘
刘刚
房建成
田珂珂
符倚伦
丁孝双
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202010505955.8A priority Critical patent/CN111854621B/zh
Publication of CN111854621A publication Critical patent/CN111854621A/zh
Application granted granted Critical
Publication of CN111854621B publication Critical patent/CN111854621B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/16—Measuring arrangements characterised by the use of optical techniques for measuring the deformation in a solid, e.g. optical strain gauge
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H9/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means
    • G01H9/004—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means using fibre optic sensors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Geometry (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Algebra (AREA)
  • Probability & Statistics with Applications (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本公开提供了机载分布式POS用光纤光栅传感器数据拟合方法，获取柔性基线上布设的光纤光栅传感器输出的波长变化量数据的振动分量，根据振动分量的瞬时振动频率和幅值将全程波长变化量数据分为若干小振动段。通过振动力学建立机翼振动时形变位移的数学模型，并基于此模型在分析波长变化量与机翼挠度之间的转换关系的基础上建立振动时波长变化量数据的数学模型。基于振动时波长变化量数据的数学模型结合机翼模态分析对波长变化量数据拟合，获得更符合机翼振动规律的更高精度的光纤光栅全程波长变化量数据，减小光纤光栅传感器分辨率有限导致的测量误差。本公开还提出机载分布式POS用光纤光栅传感器数据拟合装置。

Description

机载分布式POS用光纤光栅传感器数据拟合方法和装置
技术领域
本公开涉及航空航天技术领域，具体而言，涉及机载分布式POS用光纤光栅传感器数据拟合方法和装置。
背景技术
集成多任务遥感载荷的综合对地观测系统是目前机载对地观测的发展方向之一。对于装备多个或多种遥感载荷的高性能航空遥感系统，需要同时对位于机翼等柔性基线上的多个遥感载荷安装点的运动参数进行高精度测量。
机载分布式位置姿态测量系统(Position and Orientation System，POS)是目前用于多点高精度运动参数测量的有效手段。机载分布式POS由高精度POS系统(主节点)和多个安装于机翼下方的惯性测量单元(Inertial Measurement Unit，IMU)(子节点)组成。主节点由高精度IMU与卫星导航系统组成。由于子节点位于机翼下方，受安装方式及体积重量等限制，往往采用体积小、重量轻的中低精度IMU。因此，子IMU依赖于通过主POS的高精度导航信息以及主子节点间的高精度相对运动信息进行传递对准来获取所在节点处的高精度运动信息。然而，飞机机翼存在的复杂弹性变形导致各节点间的相对空间关系随时间变化。因此,实现分布式POS高精度传递对准的前提是实现主、子系统间柔性基线形变的高精度测量。
光纤光栅传感器是目前实现形变测量的主要手段，广泛用于桥梁等的形变测量。光纤光栅传感器输出的测量值为波长变化量(下文简称波长变化量)。在分布式POS的传递对准中，需要对数米长的机翼实现毫米级的形变测量精度。使用光纤光栅传感器测量机翼形变时，布设的相邻光纤光栅传感器之间波长间隔范围在2～4纳米左右，而光纤光栅解调仪得到的波长变化量数据的分辨率仅为皮米量级。由光纤光栅传感器分辨率不够而导致的形变测量误差显然不可忽略，将严重影响分布式POS传递对准的精度。
目前对于机翼结构的振动力学分析及光纤光栅传感器测量应变已有较多研究，但借助机翼结构的振动力学分析对光纤光栅测量的波长变化量数据进行进一步处理并用于分布式POS联合测量的研究或方法尚未见报道。
发明内容
为了解决现有技术中的技术问题，本公开实施例提供了机载分布式POS用光纤光栅传感器数据拟合方法和装置，能够减小由于光纤光栅传感器分辨率有限导致的测量误差，为机载分布式POS子节点的传递对准提供更高精度的柔性基线上各节点相对运动信息，具有操作的可行性与易用性。
第一方面，本公开实施例提供了机载分布式POS用光纤光栅传感器数据拟合方法，所述方法包括：滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值；根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定；根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度；建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
第二方面，本公开实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述的方法的步骤。
第三方面，本公开实施例提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述的方法的步骤。
第四方面，本公开实施例提供了机载分布式POS用光纤光栅传感器数据拟合装置，所述装置包括：滤除模块，用于滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量；计算模块，计算全程波长变化量中振动分量的瞬时振动频率和幅值，根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定；第一建立模块，用于根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度；第二建立模块，用于建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；拟合模块，用于使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
本发明提供的机载分布式POS用光纤光栅传感器数据拟合方法和装置，滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值；根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定；根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度；建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。该方法能够减小由于光纤光栅传感器分辨率有限导致的测量误差，为机载分布式POS子节点的传递对准提供更高精度的柔性基线上各节点相对运动信息，具有操作的可行性与易用性。
附图说明
为了更清楚地说明本公开实施例的技术方案，下面对实施例描述中所需要使用的附图作简单地介绍：
图1为本发明一个实施例中的机载分布式POS用光纤光栅传感器数据拟合方法的步骤流程示意图；
图2为本发明另一实施例中的机载分布式POS用光纤光栅传感器数据拟合方法的步骤流程示意图；
图3为本发明另一实施例中的机载分布式POS用光纤光栅传感器数据拟合方法中涉及的三角形示意图；
图4为本发明一个实施例中的机载分布式POS用光纤光栅传感器数据拟合装置的结构示意图；
图5为本发明一个实施例中的机载分布式POS用光纤光栅传感器数据拟合装置的硬件框图；
图6为本发明一个实施例中的计算机可读存储介质的示意图。
具体实施方式
下面结合附图和实施例对本申请进行进一步的详细介绍。
在下述介绍中，术语“第一”、“第二”仅为用于描述的目的，而不能理解为指示或暗示相对重要性。下述介绍提供了本公开的多个实施例，不同实施例之间可以替换或者合并组合，因此本申请也可认为包含所记载的相同和/或不同实施例的所有可能组合。因而，如果一个实施例包含特征A、B、C，另一个实施例包含特征B、D，那么本申请也应视为包括含有A、B、C、D的一个或多个所有其他可能的组合的实施例，尽管该实施例可能并未在以下内容中有明确的文字记载。
为了使本发明的目的、技术方案及优点更加清楚明白，以下通过实施例，并结合附图，对本发明机载分布式POS用光纤光栅传感器数据拟合方法和装置的具体实施方式进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
如图1所示，为一个实施例中的机载分布式POS用光纤光栅传感器数据拟合方法的流程示意图，具体包括以下步骤：
步骤11，滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值。
步骤12，根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定。
步骤13，根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度。
步骤14，建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作。
步骤15，使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
在本实施例中，滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值；针对机翼外力受力时变的特点，根据计算的全程波长变化量振动分量的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定；根据振动力学建立机翼振动时形变位移的数学模型，即挠度的数学模型；建立波长变化量与机翼挠度之间的转换关系，进而建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续操作的数据拟合；使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。该方法能够减小由于光纤光栅传感器分辨率有限导致的测量误差，为机载分布式POS子节点的传递对准提供更高精度的柔性基线上各节点相对运动信息，具有操作的可行性与易用性。
为了更加清晰与准确地理解与应用本公开所涉及的机载分布式POS用光纤光栅传感器数据拟合方法，进行以下示例。需要说明的是，本公开所保护的范围不限于以下示例。
如图2所示，为本发明另一实施例中的机载分布式POS用光纤光栅传感器数据拟合方法的步骤流程示意图。本公开涉及一种机载分布式POS用光纤光栅传感器数据拟合方法，可减小光纤光栅传感器分辨率有限导致的柔性基线形变测量误差，进而辅助光纤光栅形变测量系统得到更高精度的柔性基线形变数据，为机载分布式POS传递对准提供更高精度的柔性基线上各节点相对运动信息。
具体的，滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值。具体实施方式如下：
1)将光纤光栅传感器测量的波长变化量数据按照时间分为n段，每段100秒，逐段进行下述处理。
2)对于步骤1)分出的数据段xi(t,k)(i＝1,2,3…n)，这里的n指的是处理第n段时间段，k指的是执行2)至4)步骤的次数，k初值为零。通过求导确定xi(t,k)的所有极值点，并用三次样条插值法对极小值点形成下包络emini(t,k)，对极大值形成上包络emaxi(t,k)。
3)计算步骤2)中得到的上下包络的均值
Figure BDA0002526536140000071
从正在处理的信号段xi(t,k)中去掉均值分量，即di(t,k)＝xi(t,k)-mi(t,k)。
4)判断di(t,k)是否为内涵模态分量(Intrinsic Mode Functions,IMF)。判断标准为是否满足两个条件：在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个；在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零。若是，则从信号段xi(t,k)中去掉该分量并得到新信号xi(t,k+1)＝xi(t,k)-di(t,k)；若不是，则重复步骤2)至4)直至di(t,k)为IMF分量序列，计算残余新信号xi(t,k+1)＝xi(t,k)-di(t,k)。
5)对xi(t,k+1)重复步骤2)至4)，直至无法提取出IMF分量序列，最终得到缓变的残余分量xi(t,m)和处理的原始数据段xi(t,0)含有的m个IMF序列di(t,k),k＝1,2,3,…,m，其中，m为原始数据段xi(t,0)含有的IMF序列的个数。
6)对于按照时间划分的每个数据段xi(t,k)(i＝1,2,3…n)及对应的IMF序列di(t,k),k＝1,2,3,…,m，将IMF序列di(t,k),k＝1,2,3,…,m按照频率由高到低排列，并将di(t,k),k＝1,2,3,…,m中频率较低的几个IMF序列合并，由此获得每个数据段xi(t,k)(i＝1,2,3…n)的主要振动分量yi(t)(i＝1,2,3…n)。
7)对于每个数据段的主要振动分量yi(t)(i＝1,2,3…n)分别执行步骤8)至12)。
8)将yi(t)按照时间划分为p段信号yi(t,k)(i＝1,2,3…n；k＝1,2,3…p)，此处选定划分的时间间隔为50ms，即计算的瞬时频率和幅值为以50ms为单位的变化的瞬时频率和幅值。对于p段信号，分别执行步骤9)至12)。
9)选择一个小波函数基
Figure BDA0002526536140000081
其中ai,k为尺度因子，bi,k为伸缩因子，将
Figure BDA0002526536140000082
ai,k和bi,k作为初始值。这里的小波函数基是由同一小波母函数经过伸缩和平移得到的同一组函数序列；小波母函数是在有限时间范围内变化且平均值为零的函数。
10)计算
Figure BDA0002526536140000083
与yi(t,k)之间的相似程度，这里用小波系数CWTfi,k(ai,k,bi,k)(i＝1,2,3…n；k＝1,2,3…p)衡量该相似程度，小波系数的计算公式如下：积分区间为50ms。
Figure BDA0002526536140000084
11)判断CWTfi,k(ai,k,bi,k)是否达到极大值，若不满足，改变平移因子bi和尺度因子ai值，重复步骤4)。直至得到对于yi,k(t)而言CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k。
12)根据步骤11)得到的CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k计算yi,k(t)的瞬时振幅、瞬时频率。计算方法如下：CWTfi,k(ai,k,bi,k)可分解为实部SRi,k(t)和虚部SIi,k(t)两部分。信号的瞬时振幅可由下式确定：
Figure BDA0002526536140000085
信号的瞬时频率可由下式确定：
Figure BDA0002526536140000086
对于各通道的光纤光栅波长变化量yi(t)，均按照步骤8)按时间划分为p段信号yi(t,k)(i＝1,2,3…n；k＝1,2,3…p)，并对p段信号均执行步骤9)至12)，直至所有信号处理完成。之后对于每个通道，将得到的p段信号对应的瞬时频率和瞬时幅值按照时间连接，由此得到各通道光纤光栅波长变化量的全程瞬时频率和幅值。
进一步地，针对机翼外力受力时变的特点，根据步骤1计算的全程波长变化量振动分量的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出这些时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定的具体实施方式如下：
1)对前述计算得到的全程瞬时频率和幅值依次进行计算，若某时刻瞬时频率发生变化且变化幅值超过0.5Hz，且瞬时幅值绝对值大于2pm，则将振动开始标志位置1，并将该时刻的光纤光栅波长变化量数据的时间戳记录下来作为该段振动的开始时刻，进行步骤3.2；否则继续分析下一时刻数据的瞬时频率和幅值。
2)在振动开始标志位为1的前提下，对瞬时频率求导，若瞬时频率导数的符号发生变化或瞬时幅值小于2pm，则将振动结束标志位置一，并将该时刻的光纤光栅波长变化量数据的时间戳记录下来作为该段振动的结束时刻。之后继续返回进行下一时刻的全程瞬时频率和幅值，判断下一振动的开始时刻和结束时刻，直至该通道全程波长变化量数据分析完毕。
3)根据步骤2)计算得到的各振动段的开始时刻和结束时刻提取各振动段的波长变化量数据用于后续拟合。
更进一步地，根据振动力学建立机翼振动时形变位移即挠度的数学模型。具体实施方式包括机翼弯曲振动的运动微分方程的建立和求解两个步骤。
其中，机翼弯曲振动的运动微分方程的求解包括如下几步：
首先，使用振型叠加法将机翼弯曲振动的运动微分方程求解分解为各阶振型函数和广义坐标的求解；其次，使用脉冲响应法将机翼受复杂作用力的广义坐标微分方程的求解转化为多个单位脉冲激励下广义坐标微分方程的单位脉冲响应函数的解的叠加，并使用一元二次微分方程的求解方式对单位脉冲响应函数进行求解；再次，分析振动时刻各阶振型的求解，并由此通过振型叠加法得到机翼弯曲振动的运动微分方程的解即机翼振动时挠度的数学模型。
前述步骤中建立机翼振动时的形变位移即挠度的数学模型的步骤如下：
1)将大展弦比机翼视作等截面悬臂梁，对于机翼挠曲变形明显的一维弯曲振动，考虑外阻尼力Fd(x,t)以及材料应变的粘性阻尼力σd(x,t)两种阻尼力，建立机翼弯曲振动的运动微分方程模型。
近似于悬臂梁的机翼的弯曲振动示意图如图1所示。通过建立机翼弯曲振动的运动方程确立机翼上某点的弯曲位移随时间的变化规律。机翼弯曲振动的运动方程建立如下：
设梁长为l，单位长度的梁的质量ρ及抗弯刚度EI均为常数，建立图1所示的坐标系。在梁上距左端x处取微元段dx，在任一瞬时t，按其受力情况，微元段沿y(x,t)方向的运动方程可表示为：
Figure BDA0002526536140000101
即
Figure BDA0002526536140000102
其中，Q(x,t)为剪力，p(x,t)为分布干扰力，y(x,t)为此微元段横向位移即弯曲产生的位移。
由微元右端截面形心的力矩平衡方程可以得到
Figure BDA0002526536140000103
其中，M(x,t)为弯矩，移去二阶微量，则有
Figure BDA0002526536140000104
由材料力学中的平面假设条件，可得到弯矩与挠度曲线的关系：
Figure BDA0002526536140000111
将式(4)和式(5)代入(2)，得
Figure BDA0002526536140000112
上式即为梁弯曲振动的运动微分方程。
考虑到实际飞行过程中存在与绝对速度有关的外阻尼力Fd(x,t)以及材料应变的粘性阻尼力σd(x,t)两种阻尼力，因此结合梁弯曲振动运动微分方程(6)，建立有阻尼时的梁弯曲振动的运动微分方程如下：
当梁各点以y(x,t)横向位移做弯曲振动时，设分布的粘性外阻尼系数为c(x)，则分布的外阻尼力Fd(x,t)为
Figure BDA0002526536140000113
假定σd(x,t)与梁截面的应变速度成比例。用cs表示应变速度内阻尼系数，则材料应变的粘性阻尼力σd(x,t)可表示为
Figure BDA0002526536140000114
其中，ε(x,t)为应变，则梁的横截面上此分布阻尼力对中性轴的矩为
Figure BDA0002526536140000115
其中，z为微段长度，A为梁截面面积。根据材料力学梁的弯曲理论，应变为
Figure BDA0002526536140000116
将其代入式(7)，则得
Figure BDA0002526536140000117
由此将式(2)写为
Figure BDA0002526536140000118
则式(5)变成
Figure BDA0002526536140000121
将式(9)对x求导后，代入式(8)，得到有阻尼时的梁弯曲振动的运动微分方程如下：
Figure BDA0002526536140000125
即
Figure BDA0002526536140000122
2)对步骤1)中考虑阻尼力的机翼弯曲振动运动微分方程进行求解，具体步骤如下：
首先使用振型叠加法将机翼弯曲振动的运动微分方程求解分解为各阶振型函数和广义坐标的求解，之后使用脉冲响应法对广义坐标的微分方程进行求解，最后结合各振型函数的特点得到机翼弯曲振动的运动微分方程的解即机翼振动时挠度的数学模型。
a)使用振型叠加法将对考虑阻尼力的机翼弯曲振动运动微分方程的求解转换为各振型的运动微分方程和广义坐标乘积的叠加。进而将考虑阻尼力的机翼弯曲振动运动微分方程的求解分解为对各振型的运动微分方程和广义坐标的求解。步骤如下：
采用广义坐标表示位移曲线，将梁的位移表示为
Figure BDA0002526536140000123
式中，Yi(x)为梁的第i个主振型函数，qi(t)为第i个广义坐标。
将式(11)代入考虑阻尼力的梁弯曲振动运动微分方程式(10)，得
Figure BDA0002526536140000124
对式(12)两边各项乘以Yi(x)，并沿梁的全长进行积分，利用主振型函数对质量和刚度的正交性，即各振型之间互不影响的特性，可得
Figure BDA0002526536140000131
式中，
Figure BDA0002526536140000132
Figure BDA0002526536140000133
Figure BDA0002526536140000134
Figure BDA0002526536140000135
其中，
Figure BDA0002526536140000136
分别为第j振型对应的广义质量、广义刚度和广义载荷。对固定和自由端的梁，广义刚度
Figure BDA0002526536140000137
可表示为：
Figure BDA0002526536140000138
假设阻尼作用与质量、刚度成正比，即
c＝aρ,cs＝bE (14)
式中，a、b为比例常数，E为弹性模量，ρ为梁单位体积的质量。将式(14)代入式(13)，并利用振型的正交条件即不同阶次的振型互不影响的条件，可得
Figure BDA0002526536140000139
令
Figure BDA00025265361400001310
则ζj代表第j振型的阻尼比，可得考虑阻尼力时机翼弯曲振动的挠度微分方程j振型对应的运动方程如下：
Figure BDA00025265361400001311
b)对于a)中的各阶振型函数(17)使用脉冲响应法进行求解，即将大小随时间变化作用一段时间后停止作用的非周期激振力
Figure BDA00025265361400001312
分解为一系列强度为
Figure BDA0002526536140000141
的脉冲，先求得每一个冲量对系统激励的响应，然后利用叠加原理，对所有脉冲引起的响应进行叠加，从而得到整个非周期力
Figure BDA0002526536140000142
对系统激励的响应。
首先对单位脉冲函数及脉冲响应说明如下并求解机翼弯曲振动的挠度微分方程中第j振型对应的运动微分方程的单位脉冲响应如下：
单位脉冲函数即δ-函数，其定义为
Figure BDA0002526536140000143
考虑实际工程实现不存在无穷大的冲量，因此对单位脉冲函数做如下假设：假设某δ∈(t-τ)函数如图2所示，该函数定义为
Figure BDA0002526536140000144
当∈趋近于0时，δ∈(t-τ)函数趋近于δ-函数
Figure BDA0002526536140000145
δ-函数的单位为1/s。
在t＝τ时，非周期激振力
Figure BDA0002526536140000146
产生的冲量Iε为
Figure BDA0002526536140000147
其中，ε为此时刻的应变，当ε趋近于0时，有关系式
Figure BDA0002526536140000148
由上述单位脉冲函数的定义对振动时挠度的运动微分方程(17)求解，从而获取对应于单位脉冲作用的挠度的振动函数模型。具体步骤如下：
系统在t＝τ＝0的初始时刻受到单位冲量I＝1的作用，此时系统的运动方程为
Figure BDA0002526536140000151
其中，ζj代表第j振型的阻尼比，
Figure BDA0002526536140000152
分别为第j振型对应的广义质量、广义刚度和广义载荷，qj(t)为第j振型的广义坐标表示的梁横向振动对应的弯曲位移。
此时，系统相当于受到单位脉冲δ(t)的激励并在此脉冲作用结束之后，产生自由振动。根据动量定理可得：
Figure BDA0002526536140000153
由
Figure BDA0002526536140000154
得知系统在受到冲量I激励后的t＝0+时刻获得的速度为
Figure BDA0002526536140000155
此时系统位移尚来不及变化。故系统在受到冲量I的冲击结束后的运动方程为
Figure BDA0002526536140000156
由此可得：
Figure BDA0002526536140000157
其中，ζj代表第j振型的阻尼比，
Figure BDA0002526536140000158
分别为第j振型对应的广义质量、广义刚度和广义载荷，qj(t)为第j振型的广义坐标表示的梁横向振动对应的弯曲位移。
上述微分方程(20)即单位脉冲作用下机翼挠度的响应符合的运动微分方程，使用一元二次微分方程的求解方式对单位脉冲响应函数(20)求解。由此获得单位脉冲作用下机翼挠度的响应表达式；具体步骤如下：
方程式(20)具有如下形式的通解：
q(x,t)＝Xest
其中，X和s为待定常数，X为实数，s为复数。将上式代入式(20)，得到特征方程：
Figure BDA0002526536140000161
由此解得一对特征根如下：
Figure BDA0002526536140000162
实际情况中，0＜ζj＜1，特征根为一对共轭复根
Figure BDA0002526536140000163
式中
Figure BDA0002526536140000164
此时(20)的通解为
Figure BDA0002526536140000165
分别以X1＝C1+C2，X2＝(C1-C2)i代入上式，得到自由振动方程的解
Figure BDA0002526536140000166
其中，X1和X2为积分常数。以式(21)中的积分常数X1和X2为直角边，作图3所示直角三角形。利用该三角形可将(21)化简为
Figure BDA0002526536140000167
式中，X为振幅；ωd为振动频率，
Figure BDA0002526536140000168
为相位角，
Figure BDA0002526536140000169
为初相位。令t＝0时，系统的初位移为x0，初速度为v0，则可得系统的振幅和初相位为
Figure BDA00025265361400001610
Figure BDA00025265361400001611
由式(19)可得，对于单位脉冲响应函数，
Figure BDA00025265361400001612
所以由式(22)得系统的响应为
Figure BDA0002526536140000171
对于单位冲量I＝1，得到在单位脉冲δt的激励下的响应，记为
Figure BDA0002526536140000172
其中，h(t)称为t＝0时刻的单位脉冲响应函数。
在任意的t＝τ时刻，单位脉冲δ(t-τ)的冲击响应滞后时间τ，此时的单位脉冲响应函数为
Figure BDA0002526536140000173
在前述单位脉冲响应函数求解结果的基础上，对于线性系统，利用叠加原理，第j振型任意激励函数
Figure BDA0002526536140000174
引起的响应等于时区1≤τ≤t上所有脉冲响应的总和，即系统对任意激励
Figure BDA0002526536140000175
的响应为
Figure BDA0002526536140000176
其中，τ为积分变量，对τ积分时将t视为常量。
由此可知，各阶振型振动时广义坐标表示的机翼横向振动对应的纵向弯曲位移具有如下形式：
Figure BDA0002526536140000177
其中，参数项Cij可通过拟合确定。
c)通过前述得到非周期力作用下机翼横向振动时各阶振型微分方程的响应函数获取机翼横向振动时挠度的数学模型：
用振型叠加法求解时机翼横向振动时的挠度表示为各振型的运动微分方程和广义坐标乘积的叠加，公式如下：
Figure BDA0002526536140000181
式中，Yi(x)为梁的第i个主振型函数，qi(t)为第i个广义坐标；
光纤光栅传感器测量输出的波长变化量与挠度相关，因此需计算得到物理坐标表示的机翼横向振动对应的纵向弯曲位移。对于近似为悬臂梁的机翼，每个振型的振型函数仅与位置有关，与时间无关，因此可将机翼上位于确定位置的某个点的每阶振型的振型函数看作常值。由此根据式(25)和式(24)可知，物理坐标表示的机翼横向振动对应的纵向弯曲位移具有如下形式：
Figure BDA0002526536140000182
式中，y(t)为t时刻物理坐标表示的机翼横向振动对应的纵向弯曲位移，Dij i＝1,2,3,…；j＝1,2,3,4对于步骤2提取的每个振动段而言为待求常值。
由此得到每个小振动段物理坐标表示的机翼横向振动对应的纵向弯曲位移的模型(26)，由于光纤光栅传感器测量输出的波长变化量与挠度相关，因此可通过式(26)通过步骤4转换得到光纤光栅传感器测量输出的波长变化量在机翼横向振动时的表达式用于步骤5的数据拟合。
更进一步地，建立波长变化量与机翼挠度之间的转换关系，进而基于前述建立机翼振动时间段的光纤光栅传感器波长变化量的数学模型用于数据拟合，具体实施方式如下：
光纤光栅传感器直接提供的物理量为波长变化量，需将前述步骤3建立的每个小振动段物理坐标表示的机翼横向振动对应的纵向弯曲位移模型转换为光纤光栅波长变化量的模型用于具体实施方式5的数据拟合。具体实施过程如下：
首先使用光纤光栅形变测量中成熟普遍的Ko位移法建立光纤光栅传感器波长变化量与机翼弯曲振动的振动位移之间的关系。假设对于大展弦比柔性机翼，将其看做悬臂梁，将该梁按照距离分段，每个小段内变形为小量，并且满足线性完全弹性假设，即物体在任一瞬时的形状完全取决于它在这一瞬时受到的力，与它过去的受力情况无关。由此依据下列步骤建立梁变形后的变形位移与应变之间的关系：
结构力学中，对于承受弯曲载荷的梁结构，其微分方程为
Figure BDA0002526536140000191
式中x为结构长度的方向坐标，y(x,t)为测点处的垂直变形量即挠度，M(x,t)为结构承受的弯曲载荷，E为材料的弹性模量(杨氏模量)，I为横截面的惯性矩。x位置处的应力σ(x,t)与载荷M(x,t)之间的关系可以表示为
Figure BDA0002526536140000192
式中，c表示结构厚度的
Figure BDA0002526536140000193
根据胡克定律可得应力σ(x,t)与应变ε(x,t)之间的关系为
Figure BDA0002526536140000194
根据式(27)和(28)，结构承受的载荷M(x,t)可以用结构应变函数ε(x,t)表示为
Figure BDA0002526536140000195
则结构的微分方程可以写为
Figure BDA0002526536140000196
将大展弦比机翼视作悬臂梁，在长度为l的悬臂梁结构上，沿其长度方向等间距地布置n+1个应变传感器。在两个相邻的传感器之间的区域{xi-1,xi}(i＝1,2,3…n)内，结构的厚度函数c(x,t)和应变函数ε(x,t)可表示为线性函数
Figure BDA0002526536140000197
Figure BDA0002526536140000198
式中，{ci-1-ci}和{εi-1-εi}分别对应于{xi-1,xi}位置处的结构厚度和应变大小。{xi-1,xi}区域内，x位置处结构的斜率函数tan(θ(x))可通过式(29)积分计算得到，即
Figure BDA0002526536140000201
式中，tanθi-1为xi-1处的斜率值。
在两个相邻的传感器之间的区域{xi-1,xi}内，x位置处的结构位移函数y(x,t)可以表示为
Figure BDA0002526536140000202
其中，yi-1为xi-1处的挠度。
由此，基于悬臂梁的边界条件可知，起始段的斜率函数tan(θ0)和挠度y0均为零，可在已知起始段的斜率函数tan(θ0)和挠度y0的情况下可通过积分得到结构位移函数，即y(x,t)与ε(x,t)之间为积分关系；由于忽略温度影响时光纤布拉格光栅波长变化量与应变呈线性关系，因此光纤光栅波长变化量与物理坐标表示的机翼横向振动对应的纵向弯曲位移之间的关系可表示为成比例的积分关系。由式(26)可得，光纤布拉格光栅波长变化量λ(t)在每个机翼横向振动小振动段的模型如下：
Figure BDA0002526536140000203
Eij i＝1,2,3,…；j＝1,2,3,4为由步骤5中拟合确定的常值。
更进一步地，使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，从而获得明显振动时间段内更高精度的波长变化量数据的具体实施方式为对于前述计算得到的r个振动段，依次执行下列步骤：
1)对需要处理的柔性基线结构使用ANSYS软件进行模态分析，确定(33)中待拟合参数Eij i＝1,2,3,…；j＝1,2,3,4的初值范围；其中，每个模态的振动频率初值设定为由模态分析得到的模态频率，振动幅值的初值设定为振动挠度的最大值，由此确定几组不同的初值。
2)根据步骤1)选取的初值，对前述提取的各振动时间段的波长变化量振动分量数据与建立的振动时波长变化量的数学模型(33)使用最小二乘拟合进行拟合，使用拟合的参数结果作为下一次最小二乘拟合的初值重新计算，重复迭代，直至本次拟合的结果与前一次拟合结果的残差的差值在10-4内或达到了设定的最大迭代次数40次。
3)计算步骤2)迭代结束后的残差是否低于阈值10-4，若低于阈值，则本振动段的处理结束，进行下一振动段的处理。若残差大于阈值，则更换使用另一组步骤1)中确定的拟合参数初值，重复步骤2)步骤3)，直至步骤2)迭代结束后的残差在阈值范围内之后进行下一振动段的处理。直至所有振动段处理完毕结束。
由此得到更符合机翼振动规律的更高精度的光纤光栅传感器波长变化量数据，该数据可辅助光纤光栅形变测量系统计算得到更高精度的柔性基线形变信息，进而辅助针对多遥感载荷的机载分布式POS进行传递对准和信息融合。
综上可知，首先获取柔性基线上布设的光纤光栅传感器输出的波长变化量数据的振动分量，根据振动分量的瞬时振动频率和幅值将全程波长变化量数据分为若干小振动段。之后通过振动力学建立机翼振动时的形变位移即挠度的数学模型，并基于此模型建立振动时波长变化量数据的数学模型。最后基于振动时波长变化量数据的数学模型结合机翼模态分析对波长变化量数据进行拟合，从而获得更符合机翼振动规律的更高精度的光纤光栅全程波长变化量数据，减小光纤光栅传感器分辨率有限导致的测量误差。该数据可辅助光纤光栅形变测量系统得到更高精度的柔性基线形变信息，进而辅助针对多遥感载荷的机载分布式POS进行传递对准和信息融合。
本实施例中，针对多遥感载荷对地观测系统用机载分布式POS对柔性基线高精度相对运动信息的需求，结合机翼振动力学规律和光纤光栅应变传感器的测量特性，建立机翼振动时光纤光栅传感器输出的波长变化量的数学模型，并通过模态分析和最小二乘拟合对波长变化量进行回归拟合，得到更高精度的光纤光栅传感器波长变化量信息，从而克服了光纤光栅传感器分辨率低导致的测量误差，即提高了光纤光栅传感器的测量精度，进而辅助光纤光栅形变测量系统获取更高精度的柔性基线形变位移信息，为机载分布式POS提供更高精度的柔性基线各节点相对运动信息。
基于同一发明构思，还提供了机载分布式POS用光纤光栅传感器数据拟合装置。由于此装置解决问题的原理与前述机载分布式POS用光纤光栅传感器数据拟合方法相似，因此，该装置的实施可以按照前述方法的具体步骤实现，重复之处不再赘述。
如图4所示，为一个实施例中的机载分布式POS用光纤光栅传感器数据拟合装置的结构示意图。该机载分布式POS用光纤光栅传感器数据拟合装置10包括：滤除模块100、计算模块200、第一建立模块300、第二建立模块400和拟合模块500。
其中，滤除模块100用于滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量；计算模块200计算全程波长变化量中振动分量的瞬时振动频率和幅值，根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力近似恒定；第一建立模块300用于根据振动力学建立机翼振动时形变位移，以完成挠度的数学模型建立；第二建立模块400用于建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；拟合模块500用于使用结合机翼模态分析的最小二乘拟合方法基于建立的振动时波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
图5是图示根据本公开的实施例的机载分布式POS用光纤光栅传感器数据拟合装置的硬件框图。如图5所示，根据本公开实施例的机载分布式POS用光纤光栅传感器数据拟合装置50包括存储器501和处理器502。机载分布式POS用光纤光栅传感器数据拟合装置50中的各组件通过总线系统和/或其它形式的连接机构(未示出)互连。
存储器501用于存储非暂时性计算机可读指令。具体地，存储器501可以包括一个或多个计算机程序产品，计算机程序产品可以包括各种形式的计算机可读存储介质，例如易失性存储器和/或非易失性存储器。易失性存储器例如可以包括随机存取存储器(RAM)和/或高速缓冲存储器(cache)等。非易失性存储器例如可以包括只读存储器(ROM)、硬盘、闪存等。
处理器502可以是中央处理单元(CPU)或者具有数据处理能力和/或指令执行能力的其它形式的处理单元，并且可以控制机载分布式POS用光纤光栅传感器数据拟合装置50中的其它组件以执行期望的功能。在本公开的一个实施例中，所述处理器502用于运行存储器501中存储的计算机可读指令，使得机载分布式POS用光纤光栅传感器数据拟合装置50执行上述机载分布式POS用光纤光栅传感器数据拟合方法。机载分布式POS用光纤光栅传感器数据拟合装置与上述机载分布式POS用光纤光栅传感器数据拟合方法描述的实施例相同，在此将省略其重复描述。
图6是图示根据本公开的实施例的计算机可读存储介质的示意图。如图6所示，根据本公开实施例的计算机可读存储介质600其上存储有非暂时性计算机可读指令601。当所述非暂时性计算机可读指令601由处理器运行时，执行参照上述描述的根据本公开实施例的机载分布式POS用光纤光栅传感器数据拟合方法。
以上，根据本公开实施例的机载分布式POS用光纤光栅传感器数据拟合方法和装置，以及计算机可读存储介质，能够减小由于光纤光栅传感器分辨率有限导致的测量误差，为机载分布式POS子节点的传递对准提供更高精度的柔性基线上各节点相对运动信息，具有操作的可行性与易用性的有益效果。
以上结合具体实施例描述了本公开的基本原理，但是，需要指出的是，在本公开中提及的优点、优势、效果等仅是示例而非限制，不能认为这些优点、优势、效果等是本公开的各个实施例必须具备的。另外，上述公开的具体细节仅是为了示例的作用和便于理解的作用，而非限制，上述细节并不限制本公开为必须采用上述具体的细节来实现。
本公开中涉及的器件、装置、设备、系统的方框图仅作为例示性的例子并且不意图要求或暗示必须按照方框图示出的方式进行连接、布置、配置。如本领域技术人员将认识到的，可以按任意方式连接、布置、配置这些器件、装置、设备、系统。诸如“包括”、“包含”、“具有”等等的词语是开放性词汇，指“包括但不限于”，且可与其互换使用。这里所使用的词汇“或”和“和”指词汇“和/或”，且可与其互换使用，除非上下文明确指示不是如此。这里所使用的词汇“诸如”指词组“诸如但不限于”，且可与其互换使用。
另外，如在此使用的，在以“至少一个”开始的项的列举中使用的“或”指示分离的列举，以便例如“A、B或C的至少一个”的列举意味着A或B或C，或AB或AC或BC，或ABC(即A和B和C)。此外，措辞“示例的”不意味着描述的例子是优选的或者比其他例子更好。
还需要指出的是，在本公开的系统和方法中，各部件或各步骤是可以分解和/或重新组合的。这些分解和/或重新组合应视为本公开的等效方案。
可以不脱离由所附权利要求定义的教导的技术而进行对在此所述的技术的各种改变、替换和更改。此外，本公开的权利要求的范围不限于以上所述的处理、机器、制造、事件的组成、手段、方法和动作的具体方面。可以利用与在此所述的相应方面进行基本相同的功能或者实现基本相同的结果的当前存在的或者稍后要开发的处理、机器、制造、事件的组成、手段、方法或动作。因而，所附权利要求包括在其范围内的这样的处理、机器、制造、事件的组成、手段、方法或动作。
提供所公开的方面的以上描述以使本领域的任何技术人员能够做出或者使用本公开。对这些方面的各种修改对于本领域技术人员而言是非常显而易见的，并且在此定义的一般原理可以应用于其他方面而不脱离本公开的范围。因此，本公开不意图被限制到在此示出的方面，而是按照与在此公开的原理和新颖的特征一致的最宽范围。
为了例示和描述的目的已经给出了以上描述。此外，此描述不意图将本公开的实施例限制到在此公开的形式。尽管以上已经讨论了多个示例方面和实施例，但是本领域技术人员将认识到其某些变型、修改、改变、添加和子组合。

Claims (10)

1.机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述方法包括：
滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值；
根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并基于此从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力恒定；
根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度；
建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；
使用结合机翼模态分析的最小二乘拟合方法，基于建立的机翼振动时光纤光栅传感器波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
2.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量，并计算全程波长变化量中振动分量的瞬时振动频率和幅值包括：
将光纤光栅传感器测量的波长变化量数据按照时间分为n段，每段100秒，逐段进行处理；
对于已分出的数据段xi(t,k)，i＝1,2,3…n，通过求导确定xi(t,k)的所有极值点，并用三次样条插值法对极小值点形成下包络emini(t,k)，对极大值形成上包络emaxi(t,k)；其中，n指的是处理第n段时间段，k指的是执行处理的次数，k的初值为零；
计算上述操作中得到的上下包络的均值
Figure FDA0003223773150000021
从正在处理的数据段xi(t,k)中去掉均值分量，即di(t,k)＝xi(t,k)-mi(t,k)；
判断di(t,k)是否为内涵模态分量，内涵模态分量简记为IMF，其中，判断标准为是否满足两个条件：在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个；在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零；若是，则从数据段xi(t,k)中去掉该分量并得到新数据段xi(t,k+1)＝xi(t,k)-di(t,k)；若不是，则重复上述操作直至di(t,k)为IMF分量序列，计算残余新数据段xi(t,k+1)＝xi(t,k)-di(t,k)；
对xi(t,k+1)重复上述操作，直至无法提取出IMF分量序列，最终得到缓变的残余分量xi(t,m)和处理的原始数据段xi(t,0)含有的所有IMF序列di(t,k),k＝1,2,3,…,m，其中，m为原始数据段xi(t,0)含有的IMF序列的个数；
对于按照时间划分的每个数据段xi(t,k)及对应的IMF序列di(t,k),k＝1,2,3,…,m，将di(t,k),k＝1,2,3,…,m按照频率由高到低排列，并将di(t,k),k＝1,2,3,…,m中频率较低的几个IMF序列合并，即为每个数据段的主要振动分量yi(t)；
对于每个数据段的主要振动分量yi(t)分别执行后续操作；
将yi(t)按照时间划分为p段信号yi(t,k)，k＝1,2,3…p，此处选定为50ms，即计算的瞬时频率和幅值为以50ms为单位的变化的瞬时频率和幅值；对于p段信号，分别执行后续操作；
选择一个小波函数基
Figure FDA0003223773150000022
其中ai,k为尺度因子，bi,k为伸缩因子，将
Figure FDA0003223773150000031
ai,k和bi,k作为初始值；这里的小波函数基是由同一小波母函数经过伸缩和平移得到的同一组函数序列；小波母函数是在有限时间范围内变化且平均值为零的函数；
计算
Figure FDA0003223773150000032
与yi(t,k)之间的相似程度，这里用小波系数CWTfi,k(ai,k,bi,k)，i＝1,2,3…n；k＝1,2,3…p衡量该相似程度，小波系数的计算公式如下：积分区间为50ms；
Figure FDA0003223773150000033
判断CWTfi,k(ai,k,bi,k)是否达到极大值，若不满足，改变伸缩因子bi,k和尺度因子ai,k值，重复判断CWTfi,k(ai,k,bi,k)是否达到极大值，直至得到对于yi,k(t)而言CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k；
根据上述操作得到的CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k计算yi,k(t)的瞬时振幅、瞬时频率；计算方法如下：CWTfi,k(ai,k,bi,k)可分解为实部SRi,k(t)和虚部SIi,k(t)两部分；信号的瞬时振幅可由下式确定：
Figure FDA0003223773150000034
信号的瞬时频率可由下式确定：
Figure FDA0003223773150000035
对于各通道的光纤光栅波长变化量yi(t)，均按照将yi(t)按照时间划分为p段信号时间划分为p段信号yi(t,k)，i＝1,2,3…n；k＝1,2,3…p，并对p段信号均执行选择一个小波函数基、计算
Figure FDA0003223773150000036
与yi(t,k)之间的相似程度，判断CWTfi,k(ai,k,bi,k)是否达到极大值以及根据判断结果得到的CWTfi,k(ai,k,bi,k)达到极大值的尺度因子ai,k和伸缩因子bi,k计算yi,k(t)的瞬时振幅、瞬时频率，直至所有信号处理完成。
3.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，还包括：对于每个通道，将得到的p段信号对应的瞬时频率和瞬时振幅按照时间连接，以便得到各通道光纤光栅波长变化量的全程瞬时频率和幅值。
4.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述将光纤光栅波长变化量数据分为不同的振动时间段包括：
对得到的全程瞬时频率和幅值依次进行判断，若某时刻瞬时频率发生变化且变化幅值超过0.5Hz，且瞬时振幅绝对值大于2pm，则定义振动开始，将振动开始标志位置一，并将该时刻的光纤光栅波长变化量数据的时间戳记录下来作为该段振动的开始时刻，进行后续操作；否则继续分析下一时刻数据的瞬时频率和幅值；
在振动开始标志位为一的前提下，对瞬时频率求导，若瞬时频率导数的符号发生变化或瞬时振幅小于2pm，则将振动结束标志位置一，并将该时刻的光纤光栅波长变化量数据的时间戳记录下来作为该段振动的结束时刻；
继续返回执行对下一时间段的瞬时频率和幅值的判断，直至通道全程波长变化量数据分析完毕；
根据各振动时间段的开始时刻和结束时刻提取各振动时间段的波长变化量数据用于后续拟合。
5.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述根据振动力学建立机翼振动时形变位移的数学模型包括：机翼弯曲振动的运动微分方程的建立和求解两阶段操作；
将大展弦比机翼视作等截面悬臂梁，对于机翼挠曲变形明显的一维弯曲振动，考虑外阻尼力Fd(x,t)以及材料应变的粘性阻尼力σd(x,t)两种阻尼力，建立机翼弯曲振动的运动微分方程模型；
通过建立机翼弯曲振动的运动微分方程确立机翼上某点的弯曲位移随时间的变化规律，机翼弯曲振动的运动微分方程建立如下：
设梁长为l，单位长度的梁的质量ρ及抗弯刚度EI均为常数，在梁上距左端x处取微元段dx，在任一瞬时t，按其受力情况，微元段沿y(x,t)方向的运动方程可表示为：
Figure FDA0003223773150000051
即
Figure FDA0003223773150000052
其中，Q(x,t)为剪力，p(x,t)为分布干扰力，y(x,t)为此微元段横向位移即弯曲产生的位移；
由微元右端截面形心的力矩平衡方程可以得到
Figure FDA0003223773150000053
其中，M(x,t)为弯矩，移去二阶微量，则有
Figure FDA0003223773150000054
由材料力学中的平面假设条件，可得到弯矩与挠度曲线的关系：
Figure FDA0003223773150000055
将力矩平衡方程和弯矩与挠度曲线的关系表达式代入微元段沿y(x,t)方向的运动方程，得
Figure FDA0003223773150000056
上式即为机翼弯曲振动的运动微分方程；
考虑到实际飞行过程中存在与绝对速度有关的外阻尼力Fd(x,t)以及材料应变的粘性阻尼力σd(x,t)两种阻尼力，因此结合机翼弯曲振动的运动微分方程，建立有阻尼时的机翼弯曲振动的运动微分方程如下：
当梁各点以y(x,t)横向位移做弯曲振动时，设分布的粘性外阻尼系数为c(x)，则分布的外阻尼力Fd(x,t)为
Figure FDA0003223773150000061
假定σd(x,t)与梁截面的应变速度成比例；用cs表示应变速度内阻尼系数，则材料应变的粘性阻尼力σd(x,t)可表示为
Figure FDA0003223773150000062
其中，ε(x,t)为应变，则梁的横截面上此分布阻尼力对中性轴的矩为
Figure FDA0003223773150000063
其中，z为微段长度，A为梁截面面积；根据材料力学梁的弯曲理论，应变为
Figure FDA0003223773150000064
则得
Figure FDA0003223773150000065
由此微元段沿y(x,t)方向的运动方程可写为
Figure FDA0003223773150000066
则弯矩与挠度曲线的关系表达式变为
Figure FDA0003223773150000067
将上式对x求导后，代入微元段沿y(x,t)方向的运动方程，得到有阻尼时的机翼弯曲振动的运动微分方程如下：
Figure FDA0003223773150000068
即
Figure FDA0003223773150000069
其中，机翼弯曲振动的运动微分方程的求解得到对应的纵向位移模型的过程包括：
使用振型叠加法将机翼弯曲振动的运动微分方程求解分解为各阶振型函数和广义坐标的求解；
使用脉冲响应法将机翼受复杂作用力的广义坐标的求解转化为多个单位脉冲激励下单位脉冲响应函数的解的叠加，并使用一元二次微分方程的求解方式对单位脉冲响应函数进行求解；
分析各阶振型、广义坐标与挠度之间的关系进而获得机翼弯曲振动时挠度的解，即机翼振动时的挠度的数学模型如下：
Figure FDA0003223773150000071
式中，y(t)为t时刻物理坐标表示的机翼横向振动对应的纵向弯曲位移，Dij对于提取的每个振动时间段而言为待求常值 ， i ＝ 1,2,3, …； j ＝ 1,2,3,4 。
6.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述建立波长变化量与机翼挠度之间的转换关系，进而基于建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作包括：
使用光纤光栅形变测量中成熟普遍的Ko位移法建立光纤光栅传感器波长变化量与机翼弯曲振动的振动位移之间的关系，进而建立光纤光栅传感器波长变化量的振动模型如下：
Figure FDA0003223773150000072
其中，Eij为由需通过拟合确定的常值，其中，Eij中的i＝1,2,3……；j＝1,2,3,4。
7.根据权利要求1所述的机载分布式POS用光纤光栅传感器数据拟合方法，其特征在于，所述使用结合机翼模态分析的最小二乘拟合方法，基于建立的机翼振动时光纤光栅传感器波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据包括：
对需要处理的柔性基线结构使用ANSYS软件进行模态分析，确定预设数学模型中待拟合参数Eij的初值范围，其中，Eij中的i＝1,2,3……；j＝1,2,3,4；其中，每个模态的振动频率初值设定为由模态分析得到的模态频率，振动幅值的初值设定为振动挠度的最大值，由此确定几组不同的初值；
根据上述操作选取的初值，对于前述步骤提取的各振动时间段与建立的机翼振动时光纤光栅传感器波长变化量的数学模型使用最小二乘拟合进行拟合，使用拟合的参数结果作为下一次最小二乘拟合的初值重新计算，重复迭代，直至本次拟合的结果与前一次拟合结果的残差的差值在10-4内或达到了设定的最大迭代次数40次；
计算上述操作迭代结束后的残差是否低于阈值10-4，若低于阈值，则本振动时间段的处理结束，进行下一振动时间段的处理。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现所述权利要求1-7中任一项所述方法的步骤。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现所述权利要求1-7中任一项所述方法的步骤。
10.机载分布式POS用光纤光栅传感器数据拟合装置，其特征在于，所述装置包括：
滤除模块，用于滤除光纤光栅传感器全程输出的波长变化量数据中的缓变分量，获取全程波长变化量中的振动分量；
计算模块，计算全程波长变化量中振动分量的瞬时振动频率和幅值，根据计算的瞬时振动频率和幅值，将光纤光栅波长变化量数据分为不同的振动时间段，并从全程波长变化量振动分量中提取出前述时间段的数据；其中，定义每个振动时间段内机翼受到的气动外力恒定；
第一建立模块，用于根据振动力学建立机翼振动时形变位移的数学模型，其中，机翼振动时形变位移即挠度；
第二建立模块，用于建立波长变化量与机翼挠度之间的转换关系，进而基于前述模型建立机翼振动时光纤光栅传感器波长变化量的数学模型用于后续数据拟合操作；
拟合模块，用于使用结合机翼模态分析的最小二乘拟合方法，基于建立的机翼振动时光纤光栅传感器波长变化量的数学模型对提取的各个振动时间段的波长变化量振动分量进行拟合，以便获得更高精度的波长变化量数据。
CN202010505955.8A 2020-06-05 2020-06-05 机载分布式pos用光纤光栅传感器数据拟合方法和装置 Active CN111854621B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010505955.8A CN111854621B (zh) 2020-06-05 2020-06-05 机载分布式pos用光纤光栅传感器数据拟合方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010505955.8A CN111854621B (zh) 2020-06-05 2020-06-05 机载分布式pos用光纤光栅传感器数据拟合方法和装置

Publications (2)

Publication Number Publication Date
CN111854621A CN111854621A (zh) 2020-10-30
CN111854621B true CN111854621B (zh) 2021-10-15

Family

ID=72986066

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010505955.8A Active CN111854621B (zh) 2020-06-05 2020-06-05 机载分布式pos用光纤光栅传感器数据拟合方法和装置

Country Status (1)

Country Link
CN (1) CN111854621B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116628878B (zh) * 2023-05-22 2024-01-16 深圳大学 基于主动振动控制的减振方法、结构及计算机系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103398800A (zh) * 2013-07-20 2013-11-20 北京航空航天大学 一种用于大型结构体准分布式光纤光栅温度应变测量系统
CN105115437A (zh) * 2015-08-06 2015-12-02 中国电子科技集团公司第三十八研究所 一种机载雷达一体化天线的变形实时测量系统及方法
US9588001B2 (en) * 2014-10-17 2017-03-07 National Kaohsiung University Of Applied Sciences Pressure detecting apparatus made by 3D printing technologies being able to be used in dangerous areas
CN108413887A (zh) * 2018-02-22 2018-08-17 北京航空航天大学 光纤光栅辅助分布式pos的机翼形变测量方法、装置和平台
CN108801166A (zh) * 2018-05-29 2018-11-13 北京航空航天大学 基于悬臂梁理论的光纤光栅机翼形变测量建模及标定方法
CN108896078A (zh) * 2018-04-04 2018-11-27 天津大学 基于探测器时域响应的光纤布拉格光栅弱信号解调方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104007175B (zh) * 2014-05-09 2017-04-05 华南理工大学 一种悬臂柔性梁多裂缝损伤识别装置及方法
US10143523B2 (en) * 2015-03-31 2018-12-04 7D Surgical Inc. Systems, methods and devices for tracking and calibration of flexible instruments
CN109323659B (zh) * 2018-09-29 2024-03-29 株洲菲斯罗克光电科技股份有限公司 一种机载合成孔径雷达基线长度测量方法及装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103398800A (zh) * 2013-07-20 2013-11-20 北京航空航天大学 一种用于大型结构体准分布式光纤光栅温度应变测量系统
US9588001B2 (en) * 2014-10-17 2017-03-07 National Kaohsiung University Of Applied Sciences Pressure detecting apparatus made by 3D printing technologies being able to be used in dangerous areas
CN105115437A (zh) * 2015-08-06 2015-12-02 中国电子科技集团公司第三十八研究所 一种机载雷达一体化天线的变形实时测量系统及方法
CN108413887A (zh) * 2018-02-22 2018-08-17 北京航空航天大学 光纤光栅辅助分布式pos的机翼形变测量方法、装置和平台
CN108896078A (zh) * 2018-04-04 2018-11-27 天津大学 基于探测器时域响应的光纤布拉格光栅弱信号解调方法
CN108801166A (zh) * 2018-05-29 2018-11-13 北京航空航天大学 基于悬臂梁理论的光纤光栅机翼形变测量建模及标定方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Strain Monitoring of Combustible Gas Implosion Test Based on Fiber Bragg Grating;Gao, Lei等;《SHOCK AND VIBRATION》;20191231;全文 *
基于光纤光栅的结构应变模态测试及参数识别;朱超杰 等;《半导体光电》;20190215(第1期);第123-128页 *

Also Published As

Publication number Publication date
CN111854621A (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
Raveh CFD-based models of aerodynamic gust response
Buning et al. CFD approaches for simulation of wing-body stage separation
Zeng et al. GVT-based ground flutter test without wind tunnel
CN103970964A (zh) 一种挠性卫星模态参数在轨辨识方法
CN101131311A (zh) 一种智能化机载导弹动基座对准及标定方法
Chwalowski et al. Preliminary computational analysis of the (hirenasd) configuration in preparation for the aeroelastic prediction workshop
Bazilevs et al. Adjoint-based control of fluid-structure interaction for computational steering applications
CN109086250B (zh) 适用于带斜置光纤陀螺的mems惯组的数据融合方法
CN111854621B (zh) 机载分布式pos用光纤光栅传感器数据拟合方法和装置
Valente et al. A doublet-lattice method correction approach for high fidelity gust loads analysis
Song et al. Experimental determination of unsteady aerodynamic coefficients and flutter behavior of a rigid wing
Ritter et al. Comparison of nonlinear aeroelastic methods for maneuver simulation of very flexible aircraft
Grauer et al. Real-time parameter estimation for flexible aircraft
Kiviaho et al. Application of a time-accurate aeroelastic coupling framework to flutter-constrained design optimization
Raveh et al. Wind-tunnel study of the ARMA flutter prediction method
Suh et al. Modal filtering for control of flexible aircraft
CN112326162B (zh) 一种机载分布式pos用机翼弹性变形测量方法
Zhu et al. Research on a modeling method of wing deformation under the influence of separation and compound multi-source disturbance
Pak et al. Acceleration and Velocity Sensing from Measured Strain
Mocsányi et al. Grid and polytopic LPV modeling of aeroelastic aircraft for co-design
He et al. Geometrically nonlinear analysis for elastic beam using point interpolation meshless method
Mallik et al. Computationally efficient gust load analysis at high angles of attack
Gahan et al. Uncertainty in Wind Estimates, Part 1: Analysis using Generalized Polynomial Chaos
Zeng et al. Ground vibration test identified structure model for flutter envelope prediction
Schirmer et al. Design, implementation and experimental testing of an inertial sensor system to quantify wing deflection

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant