CN112326162A - 一种机载分布式pos用机翼弹性变形测量方法 - Google Patents

Abstract

本发明涉及一种机载分布式POS用机翼弹性变形测量方法，包括步骤：建立柔性基线上主节点和子节点的加速度关系模型；对主节点高精度陀螺仪数据以及主、子节点加速度计数据依次进行低通和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段；分别计算机翼挠曲变形时间段和机翼振动时间段中各子节点机翼对应的挠曲角；建立机翼总体挠度曲线方程，进而得到各子节点机翼对应的挠度；根据挠曲角及挠度计算机翼轴向位移变化量。本发明计算的上述信息可用于辅助机载分布式POS通过传递对准获取多个子节点的高精度运动信息，进而辅助阵列天线SAR等多任务机载对地观测遥感载荷进行高精度成像。

Description

一种机载分布式POS用机翼弹性变形测量方法

技术领域

本发明涉及弹性变形测量技术领域，特别是涉及一种机载分布式POS用机翼弹性变形测量方法，可计算机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度和挠曲角，以及机翼弯曲时不可忽略的轴向位移，辅助机载对地观测多任务遥感载荷柔性基线上多节点高精度运动信息的获取，进而辅助多任务遥感载荷对地观测进行高精度成像。

背景技术

多任务遥感载荷是目前机载对地观测的重要发展方向之一，如集成高分辨率测绘相机、成像光谱仪、大视场红外扫描仪、合成孔径雷达(Synthetic Aperture Radar，SAR)于同一载机的多任务载荷，机载分布式阵列天线SAR和柔性多基线干涉SAR等。对于装备多任务遥感载荷的高性能航空遥感系统，需要对位于机腹的主节点和机翼上分布安装的多个载荷点即子节点的运动参数实现高精度测量。

对于主节点，往往使用高精度位置姿态测量系统(Position and OrientationSystem，POS)进行绝对运动信息的高精度测量。POS是目前机载对地观测中遥感载荷获取位置、速度、姿态等运动参数的主要手段，其主要组成包括高精度惯性测量单元(InertialMeasurement Unit，IMU)、全球导航卫星系统(Global Navigation Satellite System，GNSS)接收机、导航计算机和后处理软件。但对于子节点，受体积、重量、成本等限制，通常为体积小、轻量化的中低精度IMU。因此需要依靠主节点POS的高精度位置、速度、姿态等运动参数对子节点中低精度IMU进行传递对准以实现所在节点处运动信息的精确测量。由于阵风、湍流、发动机振动和机翼内部载荷变化等影响，飞机机翼存在复杂弹性变形，导致分布式POS各节点间的空间距离(基线)呈柔性变化，因此必须对机翼弹性变形进行高精度测量才能保证传递对准的精度。

目前国内外对弹性变形信息的测量、建模和估计的方法主要分为三类。第一类方法为被动式建模，代表性方法之一是将弹性变形角理想化为马尔科夫过程，并将其增广为卡尔曼滤波的状态变量进行估计(例如李端昌,钟麦英,郭丁飞.分布式POS传递对准中的误差检测与补偿[C].第25届中国控制与决策会议论文集.2013:4194-4199.)，但机翼内部负载及外部受力复杂且时变，因此模型参数难以准确确定。另一种代表性方法是将机翼等效为悬臂梁，通过振动力学规律建模，之后将机体弹性变形角增广为卡尔曼滤波的状态变量进行估计，例如公开号为CN102621565A的专利采用ANSYS辅助建模的方法模拟机翼弹性变形。但该方法所建立的模型随飞机材质的变化而变化，并且没有考虑时变的外部气动载荷和自身载荷如油量等对飞机弹性变形中参数的影响。第二类方法是直接利用挠曲变形测量传感器测量各子节点处的弹性变形，例如视觉测量或光纤光栅形变测量系统。但这些传感器的数据更新频率较低，无法满足遥感载荷对高频实时运动参数数据的需求，且测量精度易受环境的影响。例如，光纤光栅传感器受温度影响较大，而且光纤脆弱、安装复杂；视觉测量易受遮挡和环境光强的影响。第三类方法是通过惯性测量单元计算挠曲变形信息，即通过主节点、子节点的陀螺仪或加速度计信息计算弹性变形信息。例如公开号为CN106679612A的专利将主节点、子节点的加速度计测量值之差和陀螺仪测量值之差作为量测，建立系统的非线性系统测量方程进而估计子节点处的挠曲变形和挠曲角，但该方法不可避免地用到子节点陀螺仪数据，对于分布式POS而言，子节点通常为中低精度IMU，其中陀螺仪漂移较大，可达0.1°/h至10°/h量级，而机载对地观测成像时间长达数小时，因此上述方法的测量精度无法保证。考虑到对于中低精度IMU，尽管陀螺仪漂移严重，但加速度计精度与高精度POS中的加速度计在一个数量级或者仅相差一个数量级，因此部分学者提出仅使用加速度计进行机体弹性变形估计的方法，例如公开号为CN104655132A的专利将主节点、子节点的加速度计测量值之差作为量测，建立系统的非线性系统量测方程，进而估计出机体弹性变形信息。但该方法建立的机体变形模型基于二阶马尔科夫过程，模型参数的选择多凭经验确定。在实际机载应用环境中，机体受到的内外载荷如大气扰动等复杂时变，导致模型参数难以准确确定和更新，上述问题严重影响了该方法的计算精度和实用性。由此可知，现有形变测量方法应用于机载对地观测环境时均存在各自不足。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提出一种机载分布式POS用机翼弹性变形测量方法，该方法可在仅安装分布式POS的情况下，通过主节点高精度POS和子节点加速度计计算机翼变形信息，包括机翼主要发生形变的横向弯曲对应的挠度和挠曲角，以及机翼弯曲时不可忽略的轴向位移。上述机翼形变信息可用于辅助机载分布式POS获取多个子节点的高精度运动信息，进而辅助阵列天线SAR等机载对地观测多任务遥感载荷进行高精度成像。

为解决上述技术问题，本发明采取如下的技术方案：

一种机载分布式POS用机翼弹性变形测量方法，所述方法包括以下步骤：

建立柔性基线上主节点和子节点的加速度关系模型；

对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段；

基于所述加速度关系模型分别计算机翼挠曲变形时间段和机翼振动时间段中各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角；

基于机翼结构分析和材料力学中的叠加原理建立机翼总体挠度曲线方程，并利用计算得到的所述挠曲角计算所述机翼总体挠度曲线方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度；

根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度计算机翼轴向位移变化量。

本发明与现有技术相比的优点在于：

针对机翼变形严重影响机载分布式POS子节点传递对准精度的问题，本发明提出了一种机载分布式POS用机翼弹性变形测量方法。该方法可在仅使用主节点高精度POS的角速度和加速度数据和子节点加速度计数据的情况下，结合机翼结构分析，计算得到机翼弹性变形信息，包括机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度和挠曲角，以及机翼轴向变形对应的形变位移。该方法直接根据机翼形变导致的加速度分量计算机翼变形信息，相比现有方法，有如下几个优势：一是无需计算机载环境下复杂时变的机翼内外载荷信息，且无需建立复杂的微分方程模型并进行时变参数的更新；二是无需利用子节点的陀螺仪，避免了子节点陀螺仪精度较低、漂移较大导致的形变测量误差，尤其适用于长航时的机载对地观测实验；三是与现有光纤光栅形变传感器和视觉测量等形变测量方式相比，该发明具有受环境影响小、可靠性高和安装简单的优点。

附图说明

为了更清楚地说明本公开实施例的技术方案，下面对实施例描述中所需要使用的附图作简单地介绍：

图1为本发明一个实施例中的一种机载分布式POS用机翼弹性变形测量方法的流程图；

图2为本发明一个实施例中的一种机载分布式POS用机翼弹性变形测量方法的右侧机翼挠曲变形示意图；

图3为本发明一个实施例中的一种机载分布式POS用机翼弹性变形测量方法的悬臂梁加载示意图；

图4为本发明一个实施例中的一种机载分布式POS用机翼弹性变形测量方法的右侧机翼轴向位移示意图。

具体实施方式

下面将结合附图和较佳实施例对本发明的技术方案做进一步的详细介绍。

在下述介绍中，术语“第一”、“第二”仅为用于描述的目的，而不能理解为指示或暗示相对重要性。下述介绍提供了本公开的多个实施例，不同实施例之间可以替换或者合并组合，因此本发明也可认为包含所记载的相同和/或不同实施例的所有可能组合。因而，如果一个实施例包含特征A、B、C，另一个实施例包含特征B、D，那么本发明也应视为包括含有A、B、C、D的一个或多个所有其他可能的组合的实施例，尽管该实施例可能并未在以下内容中有明确的文字记载。

为了使本发明的目的、技术方案及优点更加清楚明白，以下通过实施例，并结合附图，对本发明一种机载分布式POS用机翼弹性变形测量方法和装置的具体实施方式进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，为一个实施例中的一种机载分布式POS用机翼弹性变形测量方法的流程示意图，具体包括以下步骤：

步骤11，建立柔性基线上主节点和子节点的加速度关系模型；

步骤12，对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段；

步骤13，基于加速度关系模型分别计算机翼挠曲变形时间段和机翼振动时间段中各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角；

步骤14，基于机翼结构分析和材料力学中的叠加原理建立机翼总体挠度曲线方程，并利用计算得到的挠曲角计算机翼总体挠度曲线方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度；

步骤15，根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度计算机翼轴向位移变化量。

在本实施例中，提出一种机载分布式POS用机翼弹性变形测量方法，该方法可在仅安装分布式POS的情况下，通过主节点高精度POS和子节点加速度计计算机翼变形信息，包括机翼主要发生形变的横向弯曲对应的挠度和挠曲角，以及机翼弯曲时不可忽略的轴向位移。首先建立柔性基线上主节点和子节点的加速度关系模型。之后利用低通滤波器和高通滤波器对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段。之后基于前述建立的柔性基线上主节点和子节点的加速度关系模型，分别计算机翼挠曲变形时间段和机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角。之后基于机翼结构分析和材料力学中的叠加原理，建立机翼受垂直于轴线的横向力作用发生弯曲时位移即机翼总体挠度曲线方程，并基于前述计算得到的挠曲角通过拟合等方式计算该方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度。最后根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度，计算机翼轴向位移变化量。由此最终确定机翼弹性变形信息，其包括机翼主要发生形变的横向弯曲对应的挠度和挠曲角，以及机翼弯曲时不可忽略的轴向位移。上述机翼形变信息可用于辅助机载分布式POS获取多个子节点的高精度运动信息，进而辅助阵列天线SAR等机载对地观测多任务遥感载荷进行高精度成像。

本实施例提出的一种机载分布式POS用机翼弹性变形测量方法直接根据机翼形变导致的加速度分量计算机翼变形信息，相比现有方法，有如下几个优势：一是无需计算机载环境下复杂时变的机翼内外载荷信息，且无需建立复杂的微分方程模型并进行时变参数的更新；二是无需利用子节点的陀螺仪，避免了子节点陀螺仪精度较低、漂移较大导致的形变测量误差，尤其适用于长航时的机载对地观测实验；三是与现有光纤光栅形变传感器和视觉测量等形变测量方式相比，该发明具有受环境影响小、可靠性高和安装简单的优点。

为了更加清晰与准确地理解与应用本发明所涉及的一种机载分布式POS用机翼弹性变形测量方法，进行以下示例。需要说明的是，本发明所保护的范围不限于以下示例。

作为一种具体的实施方式，步骤11，建立柔性基线上主节点和子节点的加速度关系模型具体包括以下步骤：

首先对本发明使用的几个假设说明如下：

①由于分布式POS主要应用于大展弦比柔性机翼，翼型近似于变截面悬臂梁，因此本发明将柔性机翼视作悬臂梁；

②将分布式POS主节点、子节点视为质点；

③考虑到机载环境中，分布式POS各子节点安装于机翼上，柔性基线的变化主要体现为机翼受垂直于轴线的横向力作用发生的弯曲变形，因此仅考虑机翼受垂直于轴线的横向力作用发生弯曲对应的一维挠曲角、挠度和机翼变形时不可忽略的轴向位移，在下文中分别简称为挠曲角、挠度和轴向位移。

(1)坐标系和向量定义

基于上述假设，考虑到左右机翼对称，因此以右侧机翼为例建立坐标系及所需向量，如图2所示。

图2中建立的坐标系说明如下：O_mX_mY_mZ_m为主节点载体坐标系，坐标原点O_m为分布式POS主节点，垂直右侧机翼平面向上为Z_m轴，飞机机体主轴向前为Y_m轴，机体主轴沿右侧机翼方向为X_m轴；O_sX_sY_sZ_s为子节点载体坐标系，坐标原点O_s为分布式POS子节点，垂直右侧机翼平面往上为Z_s轴，机体主轴沿右侧机翼远离机舱方向为X_s轴，Y_s轴根据右手定则确定。O_iX_iY_iZ_i为地心惯性坐标系，坐标原点O_i为地心，X_i轴和Y_i轴在地球赤道平面内，X_i轴指向春分点，Z_i轴指向地球极轴，由右手定则确定Y_i轴指向。

图2中使用的向量说明如下：l为初始时刻机翼未发生弯曲变形时子节点相对于主节点的位移矢量在主节点载体坐标系下的投影，设l＝[l 0 0]^T；r_ms为机翼存在弹性变形时子节点相对于主节点的位移矢量在主节点载体坐标系下的投影；△l为r_ms与l的差值即机翼形变量，△l＝[0 0 l_z(x)]^T，其中l_z(x)为挠度；λ(x)为子节点处机翼挠曲角，绕子节点载体坐标系Y_s轴逆时针为正，x为机翼上子节点距离固定端的距离；r_m为主节点在地心惯性坐标系下的坐标矢量，r_s为子节点在地心惯性坐标系下的坐标矢量。

(2)基于前述定义的坐标系和向量，确定柔性基线上主节点和子节点的加速度关系模型的具体过程如下：

由上述定义可知：

r_ms＝r_s-r_m

对上式求导并结合哥氏定理，可得：

其中，

为主节点陀螺仪输出，

为子节点相对于主节点的位置矢量在主节点载体坐标系下的投影，

其中l_z(x)＝l_zchange(x)+l_zstatic(x)为挠度，挠度包括缓变分量l_zstatic(x)和振动分量l_zchange(x)，

为子节点在地心惯性坐标系下的坐标矢量的二阶导数，

为主节点在地心惯性坐标系下的坐标矢量的二阶导数，上式在主节点载体坐标系的投影为：

其中，

为子节点载体坐标系到主节点载体坐标系的转换矩阵，

为子节点在子节点载体坐标系下的坐标矢量，

为主节点在主节点载体坐标系下的坐标矢量；

由主节点、子节点三轴正交安装的加速度计测量公式可知：

其中，

为主节点加速度计输出，

为子节点加速度计输出，

为从地球坐标系到主节点载体坐标系的转换矩阵，

为从地球坐标系到子节点载体坐标系的转换矩阵，g为地球坐标系下表示的重力加速度矢量。其中，地球坐标系的定义说明如下：原点位于地心，x轴穿越本初子午线与赤道的交点，z轴穿越地球北极点，y轴穿越东经90°子午线与赤道的交点。该坐标系与地球固联。

将式(2)和式(3)代入式(1)可得：

其中，

为子节点载体坐标系到主节点载体坐标系的转换矩阵，

和

分别为主节点、子节点加速度计输出，

为主节点陀螺仪输出，

为子节点相对于主节点的位置矢量在主节点载体坐标系下的投影，

其中l_z(x)＝l_zchange(x)+l_zstatic(x)为挠度，挠度包括缓变分量l_zstatic(x)和振动分量l_zchange(x)。

考虑到机载对地观测时成像段是平稳飞行直线段，且挠曲角振动分量较小，因此

为小量，可忽略。其中，

和

分别为主节点陀螺仪输出在主节点载体坐标系下的三维矢量分量。

为主节点陀螺仪输出在主节点载体坐标系的一阶导数。而对于悬臂梁，实验结果表明，相同负载情况下，挠曲角振动分量较小时该振动分量与挠度振动分量近似呈线性比例关系。设挠曲角λ(x)＝λ_change(x)+λ_static(x)，λ_static(x)为挠曲角缓变分量，λ_change(x)为挠曲角振动分量。因此可设子节点处的挠曲角振动分量λ_change(x)与挠度振动分量l_zchange的比例为m(x)，即l_z＝-mλ_change(x)+l_zstatic(x)。此时，由于

且

因此根据前述假设可将式(4)化简并展开如下：

式(5)-(7)即为柔性基线上主节点和子节点的加速度关系模型。

作为一种具体的实施方式，步骤12，对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段包括：

(1)使用巴特沃斯低通滤波器对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据使用低通滤波去噪，可选地，频率阈值设置为30Hz，以减小高频量化噪声对原信号的影响。巴特沃斯低通滤波器的设计步骤如下：

A)确定通带边界频率Ω_p、阻带边界频率Ω_s、通带衰减a_p和阻带衰减a_s。可选地，使用巴特沃斯低通滤波器进行低通滤波去噪时，此处设定通带截止频率Ω_p＝0.01Hz、阻带开始频率Ω_s＝30Hz、通带衰减a_p＝3dB和阻带衰减a_s＝40dB。

B)确定滤波器的阶数N，阶数N为满足下列条件的最小整数：

C)利用阶数N查表求得归一化传输函数H(s)，使用

代替归一化传输函数H(s)中的s，可得实际滤波器传输函数Ha(s)。

D)使用双线性变换法由实际滤波器传输函数Ha(s)计算得到数字滤波器传递函数Ha(z)。将

代入实际滤波器传输函数Ha(s)计算即可得到数字滤波器传递函数Ha(z)，其中T为采样间隔。

E)由数字滤波器传递函数Ha(z)即可得主节点高精度陀螺仪数据以及主节点、子节点加速度计数据的差分方程，即通过前若干时刻数据计算当前时刻滤波结果的公式，并基于此差分方程计算主节点高精度陀螺仪数据以及主节点、子节点加速度计数据的低通滤波结果。

(2)使用巴特沃斯高通滤波器对前述步骤(1)得到的低通滤波去噪后的子节点X轴加速度计数据进行高通滤波，可选地，频率阈值设置为0.5Hz，以去除缓变分量，仅保留均值近似为零的振动分量。巴特沃斯高通滤波器的设计具体步骤如下：

a)确定通带边界频率Ω_p、阻带边界频率Ω_s、通带衰减a_p和阻带衰减a_s。可选地，使用巴特沃斯高通滤波器进行高通滤波去噪时，此处设定通带截止频率Ω_p＝0.5Hz、阻带开始频率Ω_s＝30Hz、通带衰减a_p＝3dB和阻带衰减a_s＝40dB。

b)确定滤波器的阶数N，阶数N为满足下列条件的最小整数：

c)利用阶数N查表求得归一化传输函数H(s)，使用

代替归一化传输函数H(s)中的s，可得实际滤波器传输函数Ha(s)。

d)使用双线性变换法由实际滤波器传输函数Ha(s)计算得到数字滤波器传递函数Ha(z)。将

代入实际滤波器传输函数Ha(s)计算即可得到数字滤波器传递函数Ha(z)，其中T为采样间隔。

e)由数字滤波器传递函数Ha(z)即可得子节点X轴加速度计数据的差分方程，即通过前若干时刻数据计算当前时刻滤波结果的公式，并基于此差分方程计算子节点X轴加速度计数据的高通滤波结果。

(3)对前述步骤(2)得到的X轴加速度计数据振动分量进行分析，根据该分析结果对机翼挠曲变形时间段和机翼振动时间段进行划分，即根据高通滤波去噪后的子节点X轴加速度计数据划分机翼挠曲变形时间段和机翼振动时间段。具体划分方法为将高通滤波去噪后的子节点X轴加速度计数据振动分量中幅值大于0.5m/s₂且持续时间超过0.3s的时间段划分为机翼振动时间段，剩余时间段认为是缓慢挠曲变形的时间段，即剩余时间段划分为机翼挠曲变形时间段。

作为一种具体的实施方式，步骤13，基于加速度关系模型分别计算机翼挠曲变形时间段和机翼振动时间段中各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角包括：

(1)机翼挠曲变形时间段中挠曲角的计算

式(5)-式(7)为复杂微分方程，求解困难。考虑到机翼缓慢变形时，

均为小量，可以忽略，因此可将式(5)-式(7)化简如下：

上述式(8)和式(10)中都含有挠曲角λ(x)。由于加速度计量化噪声及飞机平飞直行的运动特征导致

易出现极小值，若采用式(10)(10)计算挠曲角λ(x)，其计算结果易出现较大偏差。因此，本发明将式(8)(8)作为机翼挠曲变形时间段的挠曲角计算式，即根据式(8)计算机翼挠曲变形时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角。

(2)机翼振动时间段中挠曲角的计算

由于机翼振动时，

不可忽略，因此机翼振动时间段挠曲角的计算与机翼挠曲变形时间段不同。首先通过对挠曲角λ(x)的小角度近似可将式(5)-式(7)简化为：

上述三式均含有挠曲角λ(x)。与机翼挠曲变形时间段挠曲角的计算类似，加速度计量化噪声及飞机平飞直行的运动特征会导致

和

易出现接近于零的极小值，因此采用式(11)或式(12)计算挠曲角时，结果易出现较大偏差。

事实上，简化后得到的式(13)仍很难求解，在此根据成像段飞机平飞直行的运动特征将式(13)进一步化简并通过积分来计算挠曲角λ(x)。在载机平飞直行时，

较小且λ(x)通常为10度以内，因此

为小量，可以忽略，由此得到下式：

对式(14)进行积分，即可得到机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角λ(x)。

进一步地，根据式(14)即可通过积分获取机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角λ(x)的结果，具体计算步骤如下：

首先确定积分初值，具体步骤如下：取机翼振动结束之后的0.05秒时间段内的挠曲角均值作为积分初值。其中振动结束后的时间段属于机翼挠曲变形时间段，其挠曲角已通过上述步骤13的具体实施方式计算得出。

基于该积分初值，通过对机翼振动时的挠曲角速度

进行反向积分即可得到机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角计算结果。其中挠曲角速度

的具体计算方式如下：

设初始时刻挠曲角速度

为零，将对式(14)的

进行积分获得振动时挠曲角速度

并使用频率阈值为0.5Hz的高通滤波滤除挠曲角速度

中的低频分量。其中高通滤波仍使用巴特沃斯高通滤波器，巴特沃斯高通滤波器的设计步骤与前述具体实施方式中巴特沃斯高通滤波器的设计步骤类似，具体如下：

a)确定通带边界频率Ω_p、阻带边界频率Ω_s、通带衰减a_p和阻带衰减a_s。此处设定通带截止频率Ω_p＝0.01Hz、阻带开始频率Ω_s＝0.5Hz、通带衰减a_p＝3dB和阻带衰减a_s＝40dB。

b)确定滤波器的阶数N，阶数N为满足下列条件的最小整数：

c)利用阶数N查表求得归一化传输函数H(s)，使用

代替归一化传输函数H(s)中的s可得实际滤波器传输函数Ha(s)。

d)使用双线性变换法由实际滤波器传输函数Ha(s)计算得到数字滤波器传递函数Ha(z)。将

代入实际滤波器传输函数Ha(s)计算即可得到数字滤波器传递函数Ha(z)，其中T为采样间隔。

e)由数字滤波器传递函数Ha(z)即可得挠曲角速度

的差分方程，即通过前若干时刻数据计算当前时刻滤波结果的公式，并基于此差分方程计算得到挠曲角速度

的高通滤波结果。

上述机翼振动时间段连同机翼挠曲变形时间段的挠曲角计算结果即为整个测量实验过程中的挠曲角计算结果。此外，加速度计常值偏差会导致挠曲角计算结果存在一个近似为常值的误差，该误差可在实验开始时的地面静止段采用经纬仪进行标定补偿。

作为一种具体的实施方式，步骤14，基于机翼结构分析和材料力学中的叠加原理建立机翼总体挠度曲线方程，并利用计算得到的挠曲角计算机翼总体挠度曲线方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度包括：

(1)根据机翼结构分析和材料力学中的叠加原理建立机翼一维挠度曲线方程

首先根据实际机翼受力情况，将机翼受到的外力归纳为三类：第一类气动力属于分布力，合力作用在压力中心线上，大小与弦长成正比；第二类机翼质量力属于分布力，合力作用在重心线上，大小与弦长成正比；第三类部件质量力属于集中力，与部件质量和过载成正比合力作用于部件重心处。因此可将机翼受力简化为图3所示。

图3中，a为集中力作用点距离机翼固定端的距离，l为机翼展长即悬臂梁长度。C₁为距离机翼固定端的距离为x处的集中力，

为与展长成正比的分布力，F_d2(x)＝C₃为均匀分布力。C₁、C₂和C₃为待求参数。

根据材料力学中的叠加原理可知，机翼受力可看作一个集中力及两个分布力作用的叠加，因此可由结构力学计算得到的机翼总体挠度曲线方程为：

其中，EI为抗弯刚度，其中E是弹性模量，I是截面惯性矩。a为集中力作用点距离机翼固定端的距离。，l为机翼展长，即悬臂梁长度，C₁、C₂和C₃为待求参数。

(2)利用计算得到的各子节点的挠曲角对机翼总体挠度曲线方程进行求解

由于对方程(15)(15)求导即为机翼总体挠曲角方程，因此在认为机翼固定端不发生弯曲位移即式(15)常数项为零的情况下，机翼总体挠度曲线方程与机翼总体挠曲角方程的待求参数相同。因此由式(15)可知，可通过三个或三个以上子节点的挠曲角结果确定出式(15)的待求参数，也可通过至少四个子节点的挠曲角结果进行拟合来确定式(15)的待求参数，由此即可得到各子节点的挠度。此处选择采用至少四个子节点的挠曲角进行最小二乘法拟合的方法求解待求参数，进而求解子节点的挠度。最小二乘法拟合多项式的具体步骤如下：

设待拟合的子节点挠曲角的目标多项式如下：

其中，λ(x_i)为挠曲角，x_i∈R是i时刻的输入x的观测值，此处x_i∈R为机翼上各子节点距离机翼固定端的距离，R为实数集。a₀,a₁,a₂,a₃,…,a_M是M+1个待求参数。

对于给定数据点(x_i,y_i),1≤i≤N，定义目标损失函数如下：

其中，N为给定的数据点个数。

由目标损失函数最小即对目标损失函数求导并使导数等于零可得：

通过克莱姆法则求解上述方程组即可求得M+1个待求参数a₀,a₁,a₂,a₃,…,a_M。求得的待求参数a₀,a₁,a₂,a₃,…,a_M与下述挠度方程的待求参数一致，因此可根据下式确定各子节点的挠度：

其中，w(x_i)为x_i处子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度，此处x_i∈R为机翼上各子节点距离机翼固定端的距离。

作为一种具体的实施方式，步骤15，根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度计算机翼轴向位移变化量包括：

将机翼近似为悬臂梁，如图4所示。从中性面内沿x方向取出线段AB，长度为dx，位移后成为A'B'，由于中性面位移u₀引起的伸长为

则由于挠度引起的伸长为：

其中，w(x)为x处子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度。

根据极限近似公式：在x无穷小时，近似等式

成立。由此可得：

在挠曲角较小时，设机翼上距离机腹距离为x处的挠曲角为λ(x)，则对于机翼该处的微元段下式近似成立：

此时式(17)变为：

由式(18)可知，对于机翼上每个小段，其微元段长度近似公式为

轴向位移变化量近似公式为

由此将整个机翼按照0.1mm的距离间隔划分为小结构段，并从机翼(悬臂梁)固定端开始记为第一段、第二段、第三段……。此处忽略机翼形变导致的展长变化。依次通过积分确定机翼上各子节点的轴向位移，具体步骤如下：

i)计算第k段微小结构段的轴向位移变化量l_x(k)及自身微段长度l_all(k)，单位为毫米，k的初始值为1。公式如下：

其中，λ(x)＝a₀x+a₁x²+a₂x³+a₃x⁴，a_i(i＝0,1,2,3)为通过步骤14计算得到的机翼总体挠度曲线方程的待求参数，这里仅以采用四个子节点的挠曲角进行最小二乘法拟合的方法求解待求参数。

ii)计算前k段的微小结构段的轴向位移变化量L_x(k)及自身微段长度L_all(k)，公式如下：

L_x(k)＝L_x(k-1)+l_x(k)

L_all(k)＝L_all(k-1)+l_all(k)

iii)判断前k段的微小结构段的自身微段长度L_all(k)是否近似等于所求节点距离机翼固定点的长度m。若|L_all(k)-m|≤0.05，则停止计算，将此时的前k段的微小结构段的自身微段长度L_x(k)作为柔性基线上该节点的轴向位移变化量。若|L_all(k)-m|>0.05，则令k加1后重复步骤i)和步骤ii)直至前k段的微小结构段的自身微段长度L_x(k)满足|L_all(k)-m|≤0.05。

由此，可实现在仅安装分布式POS的情况下，通过主节点高精度POS和子节点加速度计计算机翼变形信息，包括机翼主要发生形变的横向弯曲对应的挠度和挠曲角，以及机翼弯曲时不可忽略的轴向位移。上述机翼形变信息可用于辅助机载分布式POS获取多个子节点的高精度运动信息，进而辅助阵列天线SAR等机载对地观测多任务遥感载荷进行高精度成像。

综上所述，针对机翼变形严重影响机载分布式POS子节点传递对准精度的问题，基于分布式POS中的主节点高精度POS和子节点加速度计数据进行传递对准所需的机翼等柔性基线变形信息的测量。首先建立柔性基线上主节点、子节点的加速度关系模型；之后对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据进行去噪，并将实验全程划分为机翼挠曲变形时间段和机翼振动时间段；其次，分别计算两种时间段中机翼受垂直于轴线的横向力作用发生弯曲时导致的挠曲角，并基于此计算各子节点受垂直于轴线的横向力作用发生弯曲时对应的挠度；最后基于前述数据计算机翼轴向位移变化量。该方法无需计算机翼内外载荷信息和建立复杂时变的微分方程模型；无需利用子节点的漂移较大陀螺仪；且与现有光纤光栅形变传感器和视觉测量等形变测量方式相比，该发明具有受环境影响小、可靠性高和安装简单的优点。通过本发明计算得到的机翼变形信息可用于辅助机载分布式POS进行高精度的传递对准进而获取机翼上多个子节点的高精度运动信息，上述运动信息可用于辅助阵列天线SAR等多任务机载对地观测遥感载荷进行高精度成像。

以上，根据本发明实施例的一种机载分布式POS用机翼弹性变形测量方法，能够基于分布式POS中的主节点高精度POS和子节点加速度计数据实现机翼等柔性基线的变形测量，从而辅助机载分布式POS传递对准获取机翼上多个子节点的高精度运动信息。

以上结合具体实施例描述了本发明的基本原理，但是，需要指出的是，在本发明中提及的优点、优势、效果等仅是示例而非限制，不能认为这些优点、优势、效果等是本发明的各个实施例必须具备的。另外，上述发明的具体细节仅是为了示例的作用和便于理解的作用，而非限制，上述细节并不限制本发明为必须采用上述具体的细节来实现。

本发明中涉及的器件、装置、设备、系统的方框图仅作为例示性的例子并且不意图要求或暗示必须按照方框图示出的方式进行连接、布置、配置。如本领域技术人员将认识到的，可以按任意方式连接、布置、配置这些器件、装置、设备、系统。诸如“包括”、“包含”、“具有”等等的词语是开放性词汇，指“包括但不限于”，且可与其互换使用。这里所使用的词汇“或”和“和”指词汇“和/或”，且可与其互换使用，除非上下文明确指示不是如此。这里所使用的词汇“诸如”指词组“诸如但不限于”，且可与其互换使用。

另外，如在此使用的，在以“至少一个”开始的项的列举中使用的“或”指示分离的列举，以便例如“A、B或C的至少一个”的列举意味着A或B或C，或AB或AC或BC，或ABC(即A和B和C)。此外，措辞“示例的”不意味着描述的例子是优选的或者比其他例子更好。

还需要指出的是，在本发明的方法中，各步骤是可以分解和/或重新组合的。这些分解和/或重新组合应视为本发明的等效方案。

可以不脱离由所附权利要求定义的教导的技术而进行对在此所述的技术的各种改变、替换和更改。此外，本发明的权利要求的范围不限于以上所述的处理、机器、制造、事件的组成、手段、方法和动作的具体方面。可以利用与在此所述的相应方面进行基本相同的功能或者实现基本相同的结果的当前存在的或者稍后要开发的处理、机器、制造、事件的组成、手段、方法或动作。因而，所附权利要求包括在其范围内的这样的处理、机器、制造、事件的组成、手段、方法或动作。

提供所公开的方面的以上描述以使本领域的任何技术人员能够做出或者使用本发明。对这些方面的各种修改对于本领域技术人员而言是非常显而易见的，并且在此定义的一般原理可以应用于其他方面而不脱离本发明的范围。因此，本发明不意图被限制到在此示出的方面，而是按照与在此公开的原理和新颖的特征一致的最宽范围。

为了例示和描述的目的已经给出了以上描述。此外，此描述不意图将本公开的实施例限制到在此公开的形式。尽管以上已经讨论了多个示例方面和实施例，但是本领域技术人员将认识到其某些变型、修改、改变、添加和子组合。

Claims (10)

1.一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述方法包括以下步骤：

建立柔性基线上主节点和子节点的加速度关系模型；

对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段；

基于所述加速度关系模型分别计算机翼挠曲变形时间段和机翼振动时间段中各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角；

基于机翼结构分析和材料力学中的叠加原理建立机翼总体挠度曲线方程，并利用计算得到的所述挠曲角计算所述机翼总体挠度曲线方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度；

根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度计算机翼轴向位移变化量。

2.根据权利要求1所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述建立柔性基线上主节点和子节点的加速度关系模型包括：

(1)坐标系和向量定义

考虑到左右机翼对称，因此以右侧机翼为例建立坐标系及所需向量，建立的坐标系说明如下：O_mX_mY_mZ_m为主节点载体坐标系，坐标原点O_m为分布式POS主节点，垂直右侧机翼平面向上为Z_m轴，飞机机体主轴向前为Y_m轴，机体主轴沿右侧机翼方向为X_m轴；O_sX_sY_sZ_s为子节点载体坐标系，坐标原点O_s为分布式POS子节点，垂直右侧机翼平面往上为Z_s轴，机体主轴沿右侧机翼远离机舱方向为X_s轴，Y_s轴根据右手定则确定；O_iX_iY_iZ_i为地心惯性坐标系，坐标原点O_i为地心，X_i轴和Y_i轴在地球赤道平面内，X_i轴指向春分点，Z_i轴指向地球极轴，由右手定则确定Y_i轴指向；

向量说明如下：l为初始时刻机翼未发生弯曲变形时子节点相对于主节点的位移矢量在主节点载体坐标系下的投影，设l＝[l 0 0]^T；r_ms为机翼存在弹性变形时子节点相对于主节点的位移矢量在主节点载体坐标系下的投影；△l为r_ms与l的差值即机翼形变量，△l＝[00 l_z(x)]^T，其中l_z(x)为挠度；λ(x)为子节点处机翼挠曲角，绕子节点载体坐标系Y_s轴逆时针为正，x为机翼上子节点距离固定端的距离；r_m为主节点在地心惯性坐标系下的坐标矢量，r_s为子节点在地心惯性坐标系下的坐标矢量；

(2)基于前述定义的坐标系和向量，确定柔性基线上主节点和子节点的加速度关系模型的具体过程如下：

由上述定义可知：

r_ms＝r_s-r_m

对上式求导并结合哥氏定理，可得：

其中，

为主节点陀螺仪输出，

为子节点相对于主节点的位置矢量在主节点载体坐标系下的投影，

其中l_z(x)＝l_zchange(x)+l_zstatic(x)为挠度，挠度包括缓变分量l_zstatic(x)和振动分量l_zchange(x)，

为子节点在地心惯性坐标系下的坐标矢量的二阶导数，

为主节点在地心惯性坐标系下的坐标矢量的二阶导数，上式在主节点载体坐标系的投影为：

其中，

为子节点载体坐标系到主节点载体坐标系的转换矩阵，

为子节点在子节点载体坐标系下的坐标矢量，

为主节点在主节点载体坐标系下的坐标矢量；

由主节点、子节点三轴正交安装的加速度计测量公式可知：

其中，

为主节点加速度计输出，

为子节点加速度计输出，

为从地球坐标系到主节点载体坐标系的转换矩阵，

为从地球坐标系到子节点载体坐标系的转换矩阵，g为地球坐标系下表示的重力加速度矢量；其中，地球坐标系的定义为：原点位于地心，x轴穿越本初子午线与赤道的交点，z轴穿越地球北极点，y轴穿越东经90°子午线与赤道的交点，地球坐标系与地球固联；

将式(2)和式(3)代入式(1)可得：

考虑到机载对地观测时成像段是平稳飞行直线段，因此

为小量，可忽略，其中

和

分别为

在主节点载体坐标系下的三维分量，即主节点正交安装的三轴陀螺仪的输出；

为主节点陀螺仪输出的一阶导数；此外，由于挠曲角振动分量较小，因此挠曲角振动分量与挠度振动分量近似呈线性比例关系，设挠曲角λ(x)＝λ_change(x)+λ_static(x)，λ_static(x)为挠曲角缓变分量，λ_change(x)为挠曲角振动分量，则可设子节点处的挠曲角振动分量λ_change(x)与挠度振动分量l_zchange的比例为m(x)，即l_z(x)＝-m(x)λ_change(x)+l_zstatic(x)，此时，由于

且

因此可将式(4)化简并展开如下：

式(5)-(7)即为柔性基线上主节点和子节点的加速度关系模型。

3.根据权利要求1或2所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据依次进行低通滤波和高通滤波去噪，并根据去噪后的数据划分机翼挠曲变形时间段和机翼振动时间段包括：

(1)使用巴特沃斯低通滤波器对主节点高精度陀螺仪数据以及主节点、子节点加速度计数据进行低通滤波去噪，所述巴特沃斯低通滤波器的设计步骤如下：

A)确定通带边界频率Ω_p、阻带边界频率Ω_s、通带衰减a_p和阻带衰减a_s；

B)确定滤波器的阶数N；

C)利用阶数N查表求得归一化传输函数H(s)，使用

代替归一化传输函数H(s)中的s，可得实际滤波器传输函数Ha(s)；

D)使用双线性变换法由实际滤波器传输函数Ha(s)计算得到数字滤波器传递函数Ha(z)；

E)根据数字滤波器传递函数Ha(z)可得主节点高精度陀螺仪数据以及主节点、子节点加速度计数据的差分方程，并基于此差分方程计算主节点高精度陀螺仪数据以及主节点、子节点加速度计数据的低通滤波结果；

(2)使用巴特沃斯高通滤波器对低通滤波去噪后的子节点X轴加速度计数据进行高通滤波去噪，所述巴特沃斯高通滤波器的设计步骤如下：

a)确定通带边界频率Ω_p、阻带边界频率Ω_s、通带衰减a_p和阻带衰减a_s；b)确定滤波器的阶数N；

c)利用阶数N查表求得归一化传输函数H(s)，使用

代替归一化传输函数H(s)中的s，可得实际滤波器传输函数Ha(s)；

d)使用双线性变换法由实际滤波器传输函数Ha(s)计算得到数字滤波器传递函数Ha(z)；

e)根据数字滤波器传递函数Ha(z)可得子节点X轴加速度计数据的差分方程，并基于此差分方程计算子节点X轴加速度计数据的高通滤波结果；

(3)根据高通滤波去噪后的子节点X轴加速度计数据划分机翼挠曲变形时间段和机翼振动时间段。

4.根据权利要求3所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，

使用巴特沃斯低通滤波器进行低通滤波去噪时，巴特沃斯低通滤波器的频率阈值设置为30Hz，通带截止频率Ω_p＝0.01Hz，阻带开始频率Ω_s＝30Hz，通带衰减a_p＝3dB，阻带衰减a_s＝40dB。

5.根据权利要求3所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，

使用巴特沃斯高通滤波器进行高通滤波去噪时，巴特沃斯高通滤波器的频率阈值设置为0.5Hz，通带截止频率Ω_p＝0.5Hz，阻带开始频率Ω_s＝30Hz，通带衰减a_p＝3dB，阻带衰减a_s＝40dB。

6.根据权利要求3所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，

将高通滤波去噪后的子节点X轴加速度计数据振动分量中幅值大于0.5m/s²且持续时间超过0.3s的时间段划分为机翼振动时间段，剩余时间段划分为机翼挠曲变形时间段。

7.根据权利要求1或2所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述基于所述加速度关系模型分别计算机翼挠曲变形时间段和机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角包括：

(1)机翼挠曲变形时间段中挠曲角的计算

考虑到机翼缓慢变形时，

均为小量忽略，将式(5)-式(7)化简如下：

根据式(8)计算机翼挠曲变形时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角；

(2)机翼振动时间段中挠曲角的计算

通过对挠曲角λ(x)的小角度近似可将式(5)-式(7)简化为：

根据成像段飞机平飞直行的运动特征将式(13)进一步化简为：

对式(14)进行积分，得到机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角。

8.根据权利要求7所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，对式(14)进行积分，得到机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角的具体计算步骤如下：

取机翼振动结束之后的0.05秒时间段内的挠曲角均值作为积分初值；

基于所述积分初值，通过对机翼振动时的挠曲角速度

进行反向积分，得到机翼振动时间段中机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角，其中挠曲角速度

的具体计算方式如下：设初始时刻挠曲角速度为零，对式(14)进行积分获得振动时的挠曲角速度

并使用频率阈值为0.5Hz的高通滤波滤除挠曲角速度

中低频分量。

9.根据权利要求1或2所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述基于机翼结构分析和材料力学中的叠加原理建立机翼总体挠度曲线方程，并利用计算得到的所述挠曲角计算所述机翼总体挠度曲线方程的待求参数，进而得到各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度包括：

(1)根据机翼结构分析和材料力学中的叠加原理建立机翼一维挠度曲线方程

根据实际机翼受力情况，将机翼受到的外力归纳为三类：第一类气动力属于分布力，合力作用在压力中心线上，大小与弦长成正比；第二类机翼质量力属于分布力，合力作用在重心线上，大小与弦长成正比；第三类部件质量力属于集中力，与部件质量和过载成正比合力作用于部件重心处；

根据材料力学中的叠加原理可知，机翼受力可看作一个集中力及两个分布力作用的叠加，因此可由结构力学计算得到的机翼总体挠度曲线方程为：

其中，EI为抗弯刚度，a为集中力作用点距离机翼固定端的距离，l为机翼展长即悬臂梁长度，C₁、C₂和C₃为待求参数；

(2)利用计算得到的各子节点的挠曲角对机翼总体挠度曲线方程进行求解由于对方程(15)(15)求导即为机翼总体挠曲角方程，因此在认为机翼固定端不发生弯曲位移即式(15)常数项为零的情况下，机翼总体挠度曲线方程与机翼总体挠曲角方程的待求参数相同，此处选择采用至少四个子节点的挠曲角进行最小二乘法拟合的方法求解待求参数；最小二乘法拟合多项式的具体步骤如下：

设待拟合的子节点挠曲角的目标多项式如下：

其中，λ(x_i)为挠曲角，x_i∈R是i时刻的输入x的观测值，x_i∈R为机翼上各子节点距离机翼固定端的距离，R为实数集，a₀,a₁,a₂,a₃,…,a_M是M+1个待求参数；

对于给定数据点(x_i,y_i),1≤i≤N，定义目标损失函数如下：

其中，N为给定数据的个数；

由目标损失函数最小即对目标损失函数求导并使导数等于零可得：

通过克莱姆法则求解上述方程组即可求得M+1个待求参数a₀,a₁,a₂,a₃,…,a_M；求得的待求参数a₀,a₁,a₂,a₃,…,a_M与下述挠度方程的待求参数一致，因此可根据下式确定各子节点的挠度：

其中，w(x_i)为x_i处子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度，x_i∈R为机翼上各子节点距离机翼固定端的距离。

10.根据权利要求1或2所述的一种机载分布式POS用机翼弹性变形测量方法，其特征在于，所述根据各子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠曲角及挠度计算机翼轴向位移变化量包括：

将机翼近似为悬臂梁，从中性面内沿x方向取出线段AB，长度为dx，位移后成为A'B'，由于中性面位移u₀引起的伸长为

则由于挠度引起的伸长为：

其中，w(x)为x处子节点机翼受垂直于轴线的横向力作用发生弯曲时对应的挠度；

根据极限近似公式：在x无穷小时，近似等式

成立，由此可得：

在挠曲角较小时，设机翼上距离机腹距离为x处的挠曲角为λ(x)，则对于机翼该处的微元段下式近似成立：

此时式(17)变为：

由式(18)可知，对于机翼上每个小段，其微元段长度近似公式为

轴向位移变化量近似公式为

由此将整个机翼按照0.1mm的距离间隔划分为小结构段，并从机翼固定端开始记为第一段、第二段、第三段……，依次通过积分确定机翼上各子节点的轴向位移，具体步骤如下：

i)计算第k段微小结构段的轴向位移变化量l_x(k)及自身微段长度l_all(k)，单位为毫米，k的初始值为1，公式如下：

其中，λ(x)＝a₀x+a₁x²+a₂x³+a₃x⁴，a_i(i＝0,1,2,3)为计算得到的机翼总体挠度曲线方程的待求参数；

ii)计算前k段的微小结构段的轴向位移变化量L_x(k)及自身微段长度L_all(k)，公式如下：

L_x(k)＝L_x(k-1)+l_x(k)

L_all(k)＝L_all(k-1)+l_all(k)

iii)判断前k段的微小结构段的自身微段长度L_all(k)是否近似等于所求子节点距离机翼固定点的长度m，若|L_all(k)-m|≤0.05，则停止计算，将此时的前k段的微小结构段的轴向位移变化量L_x(k)作为柔性基线上该子节点的轴向位移变化量；若|L_all(k)-m|>0.05，则令k加1后重复步骤i)和步骤ii)直至|L_all(k)-m|≤0.05。

Images