一种基于区域最优化拟合求取地震计传递函数的方法

1.本发明涉及地震探测技术领域，特别是涉及一种基于区域最优化拟合求取地震计传递函数的方法。


背景技术：

2.地震计作为地震探测的核心元件，其传递函数的求取直接影响地震探测数据的精度。目前广泛使用的传递函数标定方法是阶跃响应标定法和正弦波标定法。
3.阶跃响应法主要基于响应曲线拟合，该方法一方面使用的收敛算法运算量大；另一方面由于无法剔除环境噪声的影响，计算精度较低。正弦波标定法标定精度较高，但为保证其精度，必须加大频率点的密度，且频率越低周期越长。长时间测试可能导致两方面问题，一是出现环境噪声的变化，引起测试精度的降低；二是对于甚宽频带和超宽频带地震计的测试，标定工作量极其巨大。


技术实现要素：

4.为了克服现有技术的不足，本发明的目的是提供一种基于区域最优化拟合求取地震计传递函数的方法以解决地震计传递函数标定精度低的问题。
5.为实现上述目的，本发明提供了如下方案：
6.一种基于区域最优化拟合求取地震计传递函数的方法，包括：
7.获取目标地震计待求解的传递函数；
8.基于阶跃响应初步标定待求解的传递函数得到初步求解的传递函数；
9.对所述初步求解的传递函数进行傅里叶反变换，得到目标地震计的理论脉冲响应函数；
10.在所述目标地震计输入端加载脉冲信号得到目标地震计的实际脉冲响应波形；
11.根据所述理论脉冲响应函数和所述实际脉冲响应波形构建区域最优化目标函数；
12.对所述区域最优化目标函数进行求解得到标定完成的传递函数。
13.优选地，所述待求解的传递函数为：
[0014][0015]
其中，h(s)为待求解的传递函数，ξ为阻尼系数，ωn为固有频率，a为常系数，a＝ki0/m，k为地震计灵敏度，i0为施加给地震计的直流电流，m为地震计动圈质量。
[0016]
优选地，基于阶跃响应初步标定待求解的传递函数得到初步求解的传递函数，包括：
[0017]
对所述待求解的传递函数进行拉氏逆变换得到地震计阶跃响应曲线；
[0018]
利用所述地震计阶跃响应曲线上的极值点和拐点构建传递函数初步求解公式；
[0019]
根据所述传递函数初步求解公式得到初步求解的传递函数。
[0020]
优选地，所述地震计阶跃响应曲线为：
[0021][0022]
其中，e(t)为地震计阶跃响应曲线。
[0023]
优选地，所述传递函数初步求解公式为：
[0024][0025]
其中，t1为地震计阶跃响应曲线的拐点，t0为地震计阶跃响应曲线的极值点。
[0026]
优选地，所述根据所述理论脉冲响应函数和所述实际脉冲响应波形构建区域最优化目标函数，包括：
[0027]
采用公式：
[0028][0029]
构建区域最优化目标函数；其中，δ为区域最优化目标函数，bi(t)为实际脉冲响应波形，bi(t)为理论脉冲响应函数，n为bi(t)的长度。
[0030]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0031]
本发明提供了一种基于区域最优化拟合求取地震计传递函数的方法，与现有技术相比，本发明通过利用理论脉冲响应函数和实际脉冲响应波形构建区域最优化目标函数，并对其进行求解，不仅可以极大的降低传递函数标定工作的运算量，而且还可在此基础上提高传递函数的标定精度。
附图说明
[0032]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0033]
图1为本发明提供的实施例中的一种基于区域最优化拟合求取地震计传递函数的方法流程图；
[0034]
图2为本发明提供的实施例中的一种基于区域最优化拟合求取地震计传递函数的方法原理图。
具体实施方式
[0035]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0036]
在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本技术的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同
的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
[0037]
本技术的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤、过程、方法等没有限定于已列出的步骤，而是可选地还包括没有列出的步骤，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤元。
[0038]
本发明的目的是提供一种基于区域最优化拟合求取地震计传递函数的方法以解决地震计传递函数标定精度低的问题。
[0039]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0040]
请一并参阅图1-2，一种基于区域最优化拟合求取地震计传递函数的方法，包括：
[0041]
步骤1：获取目标地震计待求解的传递函数；
[0042]
在实际应用中，地震计在低频段的传递函数可用二阶系统来描述：
[0043][0044]
其中，ξ为阻尼系数，ωn为固有频率，a为常系数，
[0045]
a＝ki0/m
[0046]
k为地震计灵敏度，i0为施加给地震计的直流电流，m为地震计动圈质量。
[0047]
步骤2：基于阶跃响应初步标定待求解的传递函数得到初步求解的传递函数；
[0048]
进一步的，步骤2包括：
[0049]
对所述待求解的传递函数进行拉氏逆变换得到地震计阶跃响应曲线；
[0050]
利用所述地震计阶跃响应曲线上的极值点和拐点构建传递函数初步求解公式；
[0051]
根据所述传递函数初步求解公式得到初步求解的传递函数。
[0052]
下面结合具体的实施例对本发明上述的初步求解做进一步的说明：
[0053]
对(1)式取拉氏逆变换，根据地震计运动方程，可得：
[0054][0055]
式(2)即为地震计阶跃响应的曲线。
[0056]
对e(t)求导，令e'(t)＝0，可得响应曲线的极值点，设e'(t0)＝0，有
[0057][0058]
输出电压最大处
[0059][0060]
等价于
[0061][0062]
式(5)中，t0为响应曲线的极值点，是可测量的，只有ξ是未知量，求出ξ即可得到ωn。
[0063]
对e'(t)求导，并令e”(t)＝0，可得
[0064][0065]
由式(6)可知，t1＝2t0.t1是响应曲线的拐点。
[0066]
将t0、t1代入式(2)，相除，有：
[0067][0068]
将式(4)代入式(7)有：
[0069][0070]
求解式(8)，即可求得ξ和ωn。
[0071]
将ξ、ωn和t0值代入式(2)，即可求得灵敏度
[0072][0073]
由此，初步求出地震计传递函数。
[0074]
步骤3：对所述初步求解的传递函数进行傅里叶反变换，得到目标地震计的理论脉冲响应函数；
[0075]
步骤4：在所述目标地震计输入端加载脉冲信号得到目标地震计的实际脉冲响应
波形；
[0076]
步骤5：根据所述理论脉冲响应函数和所述实际脉冲响应波形构建区域最优化目标函数；其中，本发明的区域最优化目标函数为：
[0077][0078]
构建区域最优化目标函数；其中，δ为区域最优化目标函数，bi(t)为实际脉冲响应波形，bi(t)为理论脉冲响应函数，n为bi(t)的长度。
[0079]
步骤6：对所述区域最优化目标函数进行求解得到标定完成的传递函数。
[0080]
下面本发明结合具体的实施例对目标函数的构建过程及其求解过程做说明：
[0081]
对上述步骤求得的传递函数h(s)做傅里叶反变换，得到地震计的理论脉冲响应函数b(t)。
[0082]
在待测地震计输入端加载脉冲信号，记录地震计的脉冲响应波形b(t)，b(t)长度为n，将b(t)各点横坐标代入b(t)。
[0083]
拟求解的区域最优化目标函数为：
[0084][0085]
将上述步骤中求得的理论传递函数的各项系数设置为自变量向量(x)，在x的附近求使得目标函数δ最小的参数。
[0086]
本发明采用区域最优化算法即nelder-mead单纯形算法来对区域最优化目标函数进行求解。此算法对n维向量x使用n+1各点组成的单纯形。首先向x0添加各分量x0(i)的5％，以围绕初始估计值x0生成一个单纯形。然后使用两个上述n个向量作为单纯形的除x0之外的元素，按照以下过程反复修改单纯形：
[0087]
用x(i)表示当前单纯形中的点列表i＝1,
……
,n+1。
[0088]
按最小函数值δ(x(1))到最大函数值δ(x(n+1))的顺序对单纯形中的点进行排序。在迭代的每个步骤中，此算法都会放弃当前的最差点x(n+1)，并接受单纯形中的另一个点。
[0089]
生成反射点
[0090]
r＝2m-x(n+1)，
[0091]
其中
[0092][0093]
并计算δ(r).
[0094]
如果δ(x(1))≤δ(r)《δ(x(n))，则接受r并终止迭代。
[0095]
如果δ(r)《δ(x(1))，则计算延伸点s
[0096]
s＝m+2(m-x(n+1))
[0097]
并计算δ(s).
[0098]
如果δ(s)《δ(r)，接受s并终止迭代。
[0099]
否则，接受r并终止迭代。
[0100]
如果δ(r)≥δ(x(n))，则在m和x(n+1)或r之间执行。
[0101]
如果δ(r)《δ(x(n+1))，则计算
[0102]
c＝m+(r-m)/2
[0103]
并计算δ(c).如果δ(c)《δ(r)，则接受c并终止迭代。
[0104]
否则，继续执行步骤7.
[0105]
如果δ(r)≥δ(x(n+1))，则计算
[0106]
cc＝m+(x(n+1)-m)/2
[0107]
并计算δ(cc)。如果δ(cc)《δ(x(n+1))，则接受cc并终止迭代。
[0108]
否则，继续执行步骤7。
[0109]
计算n点
[0110]
v(i)＝x(1)+(x(i)-x(1))/2
[0111]
并计算δ(v(i))，i＝2,
…
,n+1。下一迭代中的单纯形为x(1)，v(2)，
…
，v(n+1)。
[0112]
下面本发明结合具体的应用场景对传递函数的标定过程做进一步的说明：
[0113]
第一步，阶跃响应初步标定地震计的传递函数：
[0114]
本实施例的地震数据采集器选用edas-24gn，其误差小于1/10000。在待测地震计的输入端加载+1v阶跃电压，采样率设置为800sps，记录地震计输入端的电压信号和地震计的阶跃响应波形。
[0115]
根据阶跃响应波形，读取峰值时间t0，由式(4)和式(6)可得t1，读取响应曲线e(t)的e(t0)和e(t1)值，将e(t0)和e(t1)代入式(8)，求解关于ξ的非线性一元函数。本实施例采用对分法求解ξ，对分法的上下界定为[1.001，10]。将ξ和t0代入式(5)，求出固有频率ωn。
[0116]
将ξ、ωn和t0值代入式(2)，求得灵敏度k。
[0117]
由此，初步求出地震计传递函数h(s)的各项系数a1,b1,b2,b3的理论值。
[0118][0119]
第二步，使用区域最优化算法求地震计的精确传递函数
[0120]
在待测地震计的输入端加载单位脉冲信号，采样率设置为800sps，记录地震计的脉冲响应波形b(t)。
[0121]
使用区域最优化算法，根据上述步骤中求得的理论传递函数和实测的脉冲响应函数求取精确的传递函数。
[0122]
首先，在matlab中调用impulse函数，可直接由式(11)中的各项系数求得理论脉冲响应曲线b(t)。
[0123]
然后，设最优化算法的自变量x＝[a1,b1,b2,b3]，将实测值b(t)导入matlab。自变量初始值x0设置为第一步中求得的理论值，最优化算法的目标函数为理论响应曲线和实测响应曲线的误差平方和δ，见式(10)。
[0124]
最后，对目标函数δ进行区域最优化拟合，本实施例调用fminsearch函数进行最优化计算，计算结果作为式(10)的各项系数，即为地震计的精确传递函数。
[0125]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0126]
本发明通过利用理论脉冲响应函数和实际脉冲响应波形构建区域最优化目标函数，并对其进行求解，不仅可以极大的降低传递函数标定工作的运算量，而且还可在此基础
上提高传递函数的标定精度。
[0127]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的方法而言，由于其与实施例公开的装置相对应，所以描述的比较简单，相关之处参见装置部分说明即可。
[0128]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。