利用多分辨率小波分解与重构算法将原始信号分解成不同频率的信号成分。多分辨率小波分解与重构采用Mallat提出的算法(Daubechies,1988),具体简述如下:

设{V_i}_j∈Z是空间L²(R)的一列互相嵌套的闭子空间,φ和ψ分别是相应的尺度函数和小波函数,对于数据序列c⁰构造一个函数f(t)=∑_nc⁰_nφ_0n,对某个j∈Z,j>0,函数f(t)∈V_j,有

f(t)=S₁ f(t)+D₁ f(t)=S₂ f(t)+∑²_k=1D_k f(t)=…=S_j f(t)+∑^j_k=1D_K f(t).(1)

其中,S_j f(t)=∑^∞_n=-∞c^j_nφ_j,n(t)为f(t)在尺度2^j分辨率下的连续逼近; D_j f(t)=∑^∞_n=-∞d^j_nψ_j,n(t)为f(t)在尺度2^j分辨率下的连续细节。也可以理解成:S_j f(t)为f(t)的分辨率尺度不超过2^-j的成分,D_j f(t)为f(t)的分辨率尺度介于2^-j和2^-(j-1)之间的成分。

当 {〈φ_j+1k,φ_jn〉=h(n-2k)

〈ψ_j+1k,ψ_jn〉=g(n-2k)时,(2)

则{c^j_n(t)=∑_nh(n-2k)·c^j-1_n(t),

d^j_n(t)=∑_ng(n-2k)·c^j-1_n(t).(3)

利用(1)式还可以对已分解的信号进行重构。采用重构算法对小波分解后的信号进行重构可以增加信号点数,使重构信号和原始信号的点数一样,从而得到不同分辨率尺度的变化信息。

本文采用db5小波,分解尺度j取8。db5小波是对给定支集宽度,具有极值相位和最高消失矩的紧支集小波,其相关的尺度滤波器是最小相位滤波器。db5小波的尺度函数和小波函数如图1所示。

图1 db5小波的尺度函数(a)和小波函数(b)
Fig.1 Scaling function(a)and the wavelet function(b)of db5 wavelet

对小波分解后的信号,按照正则变换小波谱的定义(Meneveau,1991)计算小波能量谱:

E(m)=∑_n|D_m,n|².(4)

其中D_m,n为小波系数。目的是以此识别出所计算时间序列的主尺度或主模(敬少群等,2002)。

对识别出的主尺度信号,进行傅里叶频谱分析,获取主尺度的主频率f_主。

按(3)式对原始观测序列进行三角拟合

y_主=Asin(2π·f_主·t+φ),(5)

并从原始序列中消除主频率的影响

y_残差=y_原-y_主.(6)

安乡井位于洞庭盆地北部的湖南省安乡市,井点坐标(29°24'N,112°10'E),井深767.32 m。该井揭露的地层较为简单,观测层是井下300～680 m间的中生代白垩系和新生代第三系沅江组一段石英(粉)砂岩多层层隙承压水。该井水位年动态平稳,水位潮汐现象不明显,自记曲线光滑,具有一定的映震能力(刘德宽等,1994)。

本文选取1993年1月至1999年12月安乡井水位日均值序列进行小波多尺度分析。分解结果见图2,由图可以看出,安乡水位日均值时间序列存在多种频率成分。当分解尺度分别为7和8时,分解后的细节部分有半年和年周期变化,且原始信号的能量主要集中在这部分。由于1996年进行洗井,干扰了2^-5～2^-7分辨率间连续细节的变化规律,半年后才恢复。

由式(4)计算的小波谱可知,2^-8分辨率下的连续细节为水位日均值序列的主要信息,傅里叶分析显示其主要成分为约365天的周期变化。

同时选取1993年1月至1999年12月安乡井气压日均值序列,进行小波多尺度分析。并对分解后同一尺度信息的气压、水位进行一元回归分析,计算相关系数。结果表明尺度8的水位与气压的相关系数最高,且该尺度的气压系数为7 mm/hPa。

图2 安乡井水位日均值序列的多尺度特征
(a)水位日均值;(b)尺度2-j分辨率下的连续逼近;(c)～(j)尺度2-j
分辨率下的连续细节
Fig.2 Multi-scale characteristic of the daily mean value of water level in Anxiang Well (a)Daily mean value of water level;(b)Continuous approximation on the resolution of scale 2-j; (c)～(j)Continuous details on the resolution of scale 2-j

分别对水位、气压日均值序列按(5)式进行拟和去除年变。水位与气压的年变振幅比为8 mm/hPa,与所计算的尺度8的气压系数相当。消除年变化后的结果见图3。

图3 安乡井去除主频信息后的水位、气压与日降雨总量分布图
Fig.3 Water level,pressure and the sum of daily rainfall after removing dominant frequency in Anxiang Well

对消除年变化后的气压、水位日均值与日累计降雨总量进行比较,笔者注意到,所研究时间段共发生过两次日降雨强度大于100 mm的强降雨,且在强降雨发生之后,在消除年变的水位日均值系列图上,均迭加有明显的方波。而且,这两个方波形态相似,起点均为倒“V”形,终点为垂直下降形,只是变化幅度随当日降雨强度增大而增加,持续天数随累计强降雨的天数增加而延长。由此表明,安乡井水位对当日降雨强度大于100 mm的强降雨引起的雨荷效应反映明显(图3)。

浅牛6井位于河北省永清县龙虎庄乡,井点坐标(39.220°N,116.421°E)。该井1996年以前为自流热水井,1996年断流后改为静水位观测井。井孔位于牛驼镇凸起北端的南侧,井深1 274 m,观测层顶板埋深1 065.3 m,观测层为震旦亚界迷雾山组白云岩岩溶裂隙承压含水层。该井水位年动态平稳,水位潮汐现象不明显,具有较强的映震能力(张宝民,肖爱华,1990; 刘喜兰等,1990)。

为了避开观测方式变更和缺数对分析结果产生的影响,取1997年8月至2007年4月的数据进

图4 浅牛6井水位日均值序列的多尺度特征
Fig.4 Multi-scale characteristic of the daily mean value in Qianniu No.6 Well

对水位日均值序列按(5)式拟和去除主频信息,结果见图5。由图中可以看出,消除主频后的时间序列存在明显的半年左右的周期变化。通过仔细对比图中各周期的形态特征,可以发现在1997年12月初至1998年2月中旬、2002年底至2003年3月初、2003年11月底至2004年3月中旬3个时间段内,在各周期的谷底迭加有其它的周期成分(图中阴影部分),表明这期间的井孔—含水层系统迭加有附加应力,从而导致在去除主频后的井孔水位不下降反而上升的情况发生。

图5 浅牛6井去除主频信息后的时序图
Fig.5 Sequence diagram of daily mean value of water level after removing dominant frequency in Qianniu No.6 Well