基于系统输出的时变特征参数辨识
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
T d Γr
中的一个质点处施加白噪声激励, 而在另一个质点 处记录振动加速度, 采样率= 1 000 H z。为保证可控 可观性, 激励及测量点取在两个端部质点处。 对响应 数据使用式 ( 2) , 式 ( 6)～ ( 8) 以及阻尼修正和滤波措 施可得到图 1 ～ 图 4 的估计结果, 估计中使用了 8 阶 的参数模型, 每段40 点数据。 图1 表明频率估计与精 确值比较吻合, 而虚假特征频率的分布具有随机性。 按式 ( 9)～ 12) 得到的CE I 分布说明了虚假特征的能
( 14)
d iag ( [ 1 0. 5 1 ] ) , 阻尼阵 = d iag [ 10 10 10 ], 刚度阵为 K exp ( - t 4) , K 为常数对称正定阵。在其
解此方程可得: ΓeT = 1. 593 62。对于 2M + 1 点的短 ) ) , 而实际阻 时序列, T = 2M f s。令 Γ= - R e ( f s ln ( Κ 尼 Γr = Γ- Α Γe Σ, 其中 Σ 为修正系数, Α为由算法决定 的常数。 对于基于二阶统计量的子空间方法, Α ≈ 0. 8。 这样 Σ 的边界满足: Γr = 0 时 Σ= 1 而且 Γr →+ ∞ 时 Σ→0。为求得修正系数 Σ, 由式 ( 14) 出发, 假定 Σ 按 下式变化 Σ= Ε=
2 特征阻尼的适度修正
由式 ( 8) 给出的阻尼因子估计是偏大的, 这是因 为信号短时相关运算导致了附加阻尼。 由于自相关计算式
2 M
式
R 11 R 12 { H = [Q 1 Q 2 ] 0 R 22
( 5)
其中最右边矩阵的对角元按绝对值递减。 若序列 { rk ( i) i= 1, 2, …, 2M - 1} 由N 个特征振动组成, 那 { 的前N 列H { (: , 1: N ) 构成主列, 对H { (: , 1: N ) 作 么H
( 4)
… … ω
rk (M ) rk (M + 1)
{ H =
rk ( 2 )
rk (M ) rk (M + 1) … rk ( 2 M - 1) θ ( = (Η M )) k ( 1 ) Η k ( 2 ) …Η k ( { 对 H ankel 矩阵 H 作Q R 分解, 假设分解具有如下形
2 M
Η k ( i) = Ρk ( i) =
∑y
j= 0
k- M + j
M + j+ i
Ν R e ( f s ln ( Κ i = i) )
( i = 1～ N )
( 8)
2 M
∑y
j= 0
k- M + j
M + j+ i
式中 Κ i 为矩阵A 的特征值, f s 为采样频率, Ξi 和 Ν i 分别为角频率和阻尼比。 为反映各特征成分的强弱, 本文从序列{ rk ( i ) i = 1, 2, …, M } 中计算对应于各 特征值的振动分量的能量指标。 构造矩阵
Ξ
表示信号中各振动分量的强弱并在时频面上建立一 种时频表示, 以此反映系统的时变特性。 为降低计算 量, 本文对H ankel 矩阵作Q R 分解, 在提高运算速度 的同时不降低参数估计的精度。
1 时变子空间方法
为从随机响应信号中辨识时变特征参数, 时变 子空间方法需使用响应信号的二阶统计量。 对于非 平稳响应, 假定信号是短时平稳的, 即特征参数变化 在一定时间内可忽略不计 ( 如一个振荡周期) 。 考虑如下的线性时变系统 x n+ 1 = A (n ) x n + B (n ) u n y n = C (n ) x n + D (n ) u n
( 9)
式中
rk ( 0 ) H = rk ( 1 ) rk (M ) rk ( 1 ) rk ( 2 ) rk (M + 1)
… … ω … … …
式中 + i = ( Κ i Κ i … Κ ) , E 是M × N 复矩阵。
T
这样, 各特征分量的展开系数计算如下 ϖT - 1ϖT { F = ( E E ) E H (: , 1) { (: , 1) 中的第 i 阶分量为 在H
( 1)
式 中 A ( n ) , B ( n ) , C ( n ) , D ( n ) 分别为N ×N 系统 矩阵, N ×P 激励矩阵, L ×N 输出矩阵以及 L ×P 直接传输矩阵。 x 为状态向量, y 为输出向量, n = 0,
1, 2, 3, …, 为时间离散点。 u 为外界激励, 假设具有
ቤተ መጻሕፍቲ ባይዱ
2 dΕ 1 ( 15) [ 1 - exp ( - Γr T ) ]
ΓrT
式 ( 15) 的超越方程并不便于直接求解, 为此对 Σ- Γr 关系采用近似处理, 这样便得到近似关系式 ( 16) 。 1 ( 16) Σ= 1 + 0. 625 ΓrT + 0. 1 ( ΓrT ) 3 结合式 ( 16) 与 Γr = Γ- Α Γe Σ 便可求得 Γr 的较好近似。 需要指出, 上述修正方法是近似的, 称为 “适度 修正” 。Α ≈ 0. 8 是为了保证修正后的阻尼比大于零, 因为 Γr = 0 时的修正量= Α Γe。
T
因此, Q R 分解只需对 由Q 1 可估计 M ×N 矩阵进行。 矩阵 C , A [ 13～ 14 ]
C = Q 1 ( 1, : ) A = Q ( 1: M - 1, : ) Q 1 ( 2: M , : )
+ 1
进行估计, rk ( i ) 实际上是其真实值与一个三角窗的 乘积。 正是三角窗的作用使得自相关估计存在附加 阻尼[ 15 ]。 对于长度给定的分析窗, 信号分量的阻尼 越大, 三角窗的影响就越小。 因此, 修正方法必须考 虑具有不同衰减常数的多分量信号的阻尼恢复。 ) , 而三 设信号分量 s ( t) = A exp ( - Γt) co s ( Ξt+ Υ 角窗在 t= 0 和 t= T 处分别为 1 和 0。为修正阻尼估
特征参数时变系统是现实存在的, 如机器人的 柔性臂、 飞机机翼等结构。 对那些因结构参数、 边界 条件等因素变化而引起特征参数显著变化的系统, 时变特性对运行具有不可忽略的影响。 因此, 识别这 类系统的时变特征参数对掌握系统的时变特性或进 行状态监测具有实际意义。 时变系统特征参数的辨识方法包括参数化及非 参数化方法[ 1～ 14 ]。 非参数化方法, 如时频分析和小 波变换, 采用直接变换的方式, 数值运算量小。 参数 化方法, 如自适应A R 建模和子空间建模, 是通过建 立输出和 ( 或) 输入信号的参数化模型来识别系统的 时变特征, 具有分辨率高的特点。 目前的参数化方法 基本是从输入输出信号中辨识频率及阻尼比[ 2, 4, 6 ] , 而仅使用输出信号的参数化方法主要用于谱估计, 目前还不能给出时变阻尼比的估计。 基于输出信号 的辨识方法在定常系统的参数辨识中具有重要的应 用价值, 近些年得到了很大的发展。 然而, 从输出信 号辨识时变特征参数的方法却相当有限, 目前常用 的做法是将用于定常系统的辨识方法在一定的假设 下转变成自适应方法, 从而实现时变特征参数辨 识[ 8 ]。 在这些自适应方法中, 时变子空间方法是一类 重要的参数化方法, 具有良好的数值性能。 本文给出 以子空间辨识为基础的时变特征估计的新方法, 通 过结合阻尼修正和数字滤波技术, 提高频率和阻尼 比的估计精度, 增强时变子空间方法的辨识功能。 本文仅考虑线性时变系统, 而且假定系统输出 信号是短时平稳的。 这样, 时变特征可通过对响应信 号在局部时间段上建模来获得。 使用分量能量指标
E = [ + 1 + 2 … + N ]
1 2
M i
式中 i = 0, 1, 2, …, 2M , 当 k - M + j + i > 2M 时, [ 12 ～ 13 ] y k - M + j + 1 = 0。 结合式 ( 1)～ ( 2) 可以得到
H = #( + 5 2
( 3)
rk (M ) rk (M + 1) rk ( 2 M )
vi = + 2 i- 1 F 2 i1
( 10)
Ρk ( 0) 2= Ρk ( 1) Ρk (M )
Ρk ( 1) Ρk ( 2) Ρk (M + 1)
Ρk (M ) Ρk (M + 1) Ρk ( 2M )
+ + 2 iF 2 i , ( i = 1, 2, …, N 2) ( 11)
ω … ( ( ) ( ) ( ) ) ( = Η …Η k 0 Η k 1 k M
( 7)
式中 Q
+ 1
为Q 1 的广义逆, Q 1 ( 1: M - 1, : ) 表示矩阵
～ M - 1 行。 由矩阵A 可估计瞬时频率和阻尼 Q 1 的1
216
振 动 工 程 学 报
第 17 卷
计, 设三角窗的等效阻尼 Γe 按如下原则等效
exp ( ∫
0
T
Γe t) d t = T 2
式中 F 2 i 分别为展开系数的第 2 i 个元素。而在时刻
n = k 附近, 各振动分量的能量指标定义为
CE I ( k , i) =
vi vi ( i = 1, 2, …, N
γT
2) ( 12)
对于分析窗内的响应信号, 可等效认为它是某定常 系统在冲击激励下的响应 ( 在二阶统计的意义上) , 即对于 i> 0, Ρk ( i) = 0, 这样式 ( 3) 可简化为 { θ H = # ( 式中
rk ( 1 ) rk ( 2 ) rk ( 3 )
从计算式 ( 12) 可以看出, CE I 反映了响应信号 中各特征分量在分析窗内的能量, 因而可用来评价 特征分量的振动强弱。 由于预先选择的模型阶次一 般要高于实际系统的可观测阶次, 这样在估计的特 征值中就存在虚假成分。 然而, 根据CE I 可在时频面 内建立一种时频表示, 它可以反映时变特征的连续变 化过程。 对于虚假的特征分量, 可根据其在时频面内 的无规则分布特征以及能量弱小的特点加以识别。
第 17 卷第 2 期 2004 年 6 月
振 动 工 程 学 报
Jou rna l of V ib ra t ion Eng ineering
Vol . 17 N o. 2 J un. 2004
基于系统输出的时变特征参数辨识
张志谊1 胡 芳2 樊江玲1 华宏星1
Ξ
( 1 上海交通大学振动、 冲击、 噪声研究所 上海, 200030) ( 2 武汉理工大学材料科学与工程学院 武汉, 430070)
白噪声特性。 现以离散时间 n = k 为中心, 从输出响 应 y 中前后各取M 点构成一个短时序列 v ( k ) = {y k - M , y k - M + 1 , …, y k , …, y k + M - 1 , y k + M }T , 同时认为 A (n ) , B ( n ) , C ( n ) , D ( n ) 在短时窗内几乎为常数矩 阵, 即 A ( n ) ≈ A , B ( n ) ≈ B , C ( n ) ≈ C , D ( n ) ≈ D 考虑L = 1 的情形, 并由序列 v ( k ) 定义如下的短 时相关序列
Q R 分解可以得到 { H (: , 1: N ) = Q 1R 11 ( 6)
rk ( i ) =
∑y
j- 0
k- M + j
y k-
M + j+ i
( 13)
仅使用短时序列 v ( k ) = {y k - M , y k - M + 1 , …, y k , …,
y k+ M - 1, y k+ M }
摘 要 在子空间辨识方法的基础上给出时变特征参数辨识的新方法。 将特征阻尼修正与数字滤波相结合, 由输 出信号估计特征频率和阻尼比。对于系统输出, 假定其具有短时平稳特性, 通过对响应信号在局部时间段上建模得 到时变特征。使用分量能量指标指示信号中各振动分量的强弱, 并通过能量指标在时频面上建立一种时频表示, 以 此反映系统的时频特性。 对 H ankel 矩阵使用Q R 分解, 在提高运算速度的同时不降低参数估计的精度。 关键词: 参数辨识; 时变系统; 子空间建模 中图分类号: O 324; TN 911
收稿日期: 2003204204; 修改稿收到日期: 2003211213
第 2 期
2 M
张志谊等: 基于系统输出的时变特征参数辨识
215
rk ( i ) =
∑y
j= 0
k- M + j
y kx ku k-
比
M + j+ i
Ξi =
( 2)
f s ln ( Κ i) f s ln ( Κ i)
中的一个质点处施加白噪声激励, 而在另一个质点 处记录振动加速度, 采样率= 1 000 H z。为保证可控 可观性, 激励及测量点取在两个端部质点处。 对响应 数据使用式 ( 2) , 式 ( 6)～ ( 8) 以及阻尼修正和滤波措 施可得到图 1 ～ 图 4 的估计结果, 估计中使用了 8 阶 的参数模型, 每段40 点数据。 图1 表明频率估计与精 确值比较吻合, 而虚假特征频率的分布具有随机性。 按式 ( 9)～ 12) 得到的CE I 分布说明了虚假特征的能
( 14)
d iag ( [ 1 0. 5 1 ] ) , 阻尼阵 = d iag [ 10 10 10 ], 刚度阵为 K exp ( - t 4) , K 为常数对称正定阵。在其
解此方程可得: ΓeT = 1. 593 62。对于 2M + 1 点的短 ) ) , 而实际阻 时序列, T = 2M f s。令 Γ= - R e ( f s ln ( Κ 尼 Γr = Γ- Α Γe Σ, 其中 Σ 为修正系数, Α为由算法决定 的常数。 对于基于二阶统计量的子空间方法, Α ≈ 0. 8。 这样 Σ 的边界满足: Γr = 0 时 Σ= 1 而且 Γr →+ ∞ 时 Σ→0。为求得修正系数 Σ, 由式 ( 14) 出发, 假定 Σ 按 下式变化 Σ= Ε=
2 特征阻尼的适度修正
由式 ( 8) 给出的阻尼因子估计是偏大的, 这是因 为信号短时相关运算导致了附加阻尼。 由于自相关计算式
2 M
式
R 11 R 12 { H = [Q 1 Q 2 ] 0 R 22
( 5)
其中最右边矩阵的对角元按绝对值递减。 若序列 { rk ( i) i= 1, 2, …, 2M - 1} 由N 个特征振动组成, 那 { 的前N 列H { (: , 1: N ) 构成主列, 对H { (: , 1: N ) 作 么H
( 4)
… … ω
rk (M ) rk (M + 1)
{ H =
rk ( 2 )
rk (M ) rk (M + 1) … rk ( 2 M - 1) θ ( = (Η M )) k ( 1 ) Η k ( 2 ) …Η k ( { 对 H ankel 矩阵 H 作Q R 分解, 假设分解具有如下形
2 M
Η k ( i) = Ρk ( i) =
∑y
j= 0
k- M + j
M + j+ i
Ν R e ( f s ln ( Κ i = i) )
( i = 1～ N )
( 8)
2 M
∑y
j= 0
k- M + j
M + j+ i
式中 Κ i 为矩阵A 的特征值, f s 为采样频率, Ξi 和 Ν i 分别为角频率和阻尼比。 为反映各特征成分的强弱, 本文从序列{ rk ( i ) i = 1, 2, …, M } 中计算对应于各 特征值的振动分量的能量指标。 构造矩阵
Ξ
表示信号中各振动分量的强弱并在时频面上建立一 种时频表示, 以此反映系统的时变特性。 为降低计算 量, 本文对H ankel 矩阵作Q R 分解, 在提高运算速度 的同时不降低参数估计的精度。
1 时变子空间方法
为从随机响应信号中辨识时变特征参数, 时变 子空间方法需使用响应信号的二阶统计量。 对于非 平稳响应, 假定信号是短时平稳的, 即特征参数变化 在一定时间内可忽略不计 ( 如一个振荡周期) 。 考虑如下的线性时变系统 x n+ 1 = A (n ) x n + B (n ) u n y n = C (n ) x n + D (n ) u n
( 9)
式中
rk ( 0 ) H = rk ( 1 ) rk (M ) rk ( 1 ) rk ( 2 ) rk (M + 1)
… … ω … … …
式中 + i = ( Κ i Κ i … Κ ) , E 是M × N 复矩阵。
T
这样, 各特征分量的展开系数计算如下 ϖT - 1ϖT { F = ( E E ) E H (: , 1) { (: , 1) 中的第 i 阶分量为 在H
( 1)
式 中 A ( n ) , B ( n ) , C ( n ) , D ( n ) 分别为N ×N 系统 矩阵, N ×P 激励矩阵, L ×N 输出矩阵以及 L ×P 直接传输矩阵。 x 为状态向量, y 为输出向量, n = 0,
1, 2, 3, …, 为时间离散点。 u 为外界激励, 假设具有
ቤተ መጻሕፍቲ ባይዱ
2 dΕ 1 ( 15) [ 1 - exp ( - Γr T ) ]
ΓrT
式 ( 15) 的超越方程并不便于直接求解, 为此对 Σ- Γr 关系采用近似处理, 这样便得到近似关系式 ( 16) 。 1 ( 16) Σ= 1 + 0. 625 ΓrT + 0. 1 ( ΓrT ) 3 结合式 ( 16) 与 Γr = Γ- Α Γe Σ 便可求得 Γr 的较好近似。 需要指出, 上述修正方法是近似的, 称为 “适度 修正” 。Α ≈ 0. 8 是为了保证修正后的阻尼比大于零, 因为 Γr = 0 时的修正量= Α Γe。
T
因此, Q R 分解只需对 由Q 1 可估计 M ×N 矩阵进行。 矩阵 C , A [ 13～ 14 ]
C = Q 1 ( 1, : ) A = Q ( 1: M - 1, : ) Q 1 ( 2: M , : )
+ 1
进行估计, rk ( i ) 实际上是其真实值与一个三角窗的 乘积。 正是三角窗的作用使得自相关估计存在附加 阻尼[ 15 ]。 对于长度给定的分析窗, 信号分量的阻尼 越大, 三角窗的影响就越小。 因此, 修正方法必须考 虑具有不同衰减常数的多分量信号的阻尼恢复。 ) , 而三 设信号分量 s ( t) = A exp ( - Γt) co s ( Ξt+ Υ 角窗在 t= 0 和 t= T 处分别为 1 和 0。为修正阻尼估
特征参数时变系统是现实存在的, 如机器人的 柔性臂、 飞机机翼等结构。 对那些因结构参数、 边界 条件等因素变化而引起特征参数显著变化的系统, 时变特性对运行具有不可忽略的影响。 因此, 识别这 类系统的时变特征参数对掌握系统的时变特性或进 行状态监测具有实际意义。 时变系统特征参数的辨识方法包括参数化及非 参数化方法[ 1～ 14 ]。 非参数化方法, 如时频分析和小 波变换, 采用直接变换的方式, 数值运算量小。 参数 化方法, 如自适应A R 建模和子空间建模, 是通过建 立输出和 ( 或) 输入信号的参数化模型来识别系统的 时变特征, 具有分辨率高的特点。 目前的参数化方法 基本是从输入输出信号中辨识频率及阻尼比[ 2, 4, 6 ] , 而仅使用输出信号的参数化方法主要用于谱估计, 目前还不能给出时变阻尼比的估计。 基于输出信号 的辨识方法在定常系统的参数辨识中具有重要的应 用价值, 近些年得到了很大的发展。 然而, 从输出信 号辨识时变特征参数的方法却相当有限, 目前常用 的做法是将用于定常系统的辨识方法在一定的假设 下转变成自适应方法, 从而实现时变特征参数辨 识[ 8 ]。 在这些自适应方法中, 时变子空间方法是一类 重要的参数化方法, 具有良好的数值性能。 本文给出 以子空间辨识为基础的时变特征估计的新方法, 通 过结合阻尼修正和数字滤波技术, 提高频率和阻尼 比的估计精度, 增强时变子空间方法的辨识功能。 本文仅考虑线性时变系统, 而且假定系统输出 信号是短时平稳的。 这样, 时变特征可通过对响应信 号在局部时间段上建模来获得。 使用分量能量指标
E = [ + 1 + 2 … + N ]
1 2
M i
式中 i = 0, 1, 2, …, 2M , 当 k - M + j + i > 2M 时, [ 12 ～ 13 ] y k - M + j + 1 = 0。 结合式 ( 1)～ ( 2) 可以得到
H = #( + 5 2
( 3)
rk (M ) rk (M + 1) rk ( 2 M )
vi = + 2 i- 1 F 2 i1
( 10)
Ρk ( 0) 2= Ρk ( 1) Ρk (M )
Ρk ( 1) Ρk ( 2) Ρk (M + 1)
Ρk (M ) Ρk (M + 1) Ρk ( 2M )
+ + 2 iF 2 i , ( i = 1, 2, …, N 2) ( 11)
ω … ( ( ) ( ) ( ) ) ( = Η …Η k 0 Η k 1 k M
( 7)
式中 Q
+ 1
为Q 1 的广义逆, Q 1 ( 1: M - 1, : ) 表示矩阵
～ M - 1 行。 由矩阵A 可估计瞬时频率和阻尼 Q 1 的1
216
振 动 工 程 学 报
第 17 卷
计, 设三角窗的等效阻尼 Γe 按如下原则等效
exp ( ∫
0
T
Γe t) d t = T 2
式中 F 2 i 分别为展开系数的第 2 i 个元素。而在时刻
n = k 附近, 各振动分量的能量指标定义为
CE I ( k , i) =
vi vi ( i = 1, 2, …, N
γT
2) ( 12)
对于分析窗内的响应信号, 可等效认为它是某定常 系统在冲击激励下的响应 ( 在二阶统计的意义上) , 即对于 i> 0, Ρk ( i) = 0, 这样式 ( 3) 可简化为 { θ H = # ( 式中
rk ( 1 ) rk ( 2 ) rk ( 3 )
从计算式 ( 12) 可以看出, CE I 反映了响应信号 中各特征分量在分析窗内的能量, 因而可用来评价 特征分量的振动强弱。 由于预先选择的模型阶次一 般要高于实际系统的可观测阶次, 这样在估计的特 征值中就存在虚假成分。 然而, 根据CE I 可在时频面 内建立一种时频表示, 它可以反映时变特征的连续变 化过程。 对于虚假的特征分量, 可根据其在时频面内 的无规则分布特征以及能量弱小的特点加以识别。
第 17 卷第 2 期 2004 年 6 月
振 动 工 程 学 报
Jou rna l of V ib ra t ion Eng ineering
Vol . 17 N o. 2 J un. 2004
基于系统输出的时变特征参数辨识
张志谊1 胡 芳2 樊江玲1 华宏星1
Ξ
( 1 上海交通大学振动、 冲击、 噪声研究所 上海, 200030) ( 2 武汉理工大学材料科学与工程学院 武汉, 430070)
白噪声特性。 现以离散时间 n = k 为中心, 从输出响 应 y 中前后各取M 点构成一个短时序列 v ( k ) = {y k - M , y k - M + 1 , …, y k , …, y k + M - 1 , y k + M }T , 同时认为 A (n ) , B ( n ) , C ( n ) , D ( n ) 在短时窗内几乎为常数矩 阵, 即 A ( n ) ≈ A , B ( n ) ≈ B , C ( n ) ≈ C , D ( n ) ≈ D 考虑L = 1 的情形, 并由序列 v ( k ) 定义如下的短 时相关序列
Q R 分解可以得到 { H (: , 1: N ) = Q 1R 11 ( 6)
rk ( i ) =
∑y
j- 0
k- M + j
y k-
M + j+ i
( 13)
仅使用短时序列 v ( k ) = {y k - M , y k - M + 1 , …, y k , …,
y k+ M - 1, y k+ M }
摘 要 在子空间辨识方法的基础上给出时变特征参数辨识的新方法。 将特征阻尼修正与数字滤波相结合, 由输 出信号估计特征频率和阻尼比。对于系统输出, 假定其具有短时平稳特性, 通过对响应信号在局部时间段上建模得 到时变特征。使用分量能量指标指示信号中各振动分量的强弱, 并通过能量指标在时频面上建立一种时频表示, 以 此反映系统的时频特性。 对 H ankel 矩阵使用Q R 分解, 在提高运算速度的同时不降低参数估计的精度。 关键词: 参数辨识; 时变系统; 子空间建模 中图分类号: O 324; TN 911
收稿日期: 2003204204; 修改稿收到日期: 2003211213
第 2 期
2 M
张志谊等: 基于系统输出的时变特征参数辨识
215
rk ( i ) =
∑y
j= 0
k- M + j
y kx ku k-
比
M + j+ i
Ξi =
( 2)
f s ln ( Κ i) f s ln ( Κ i)