全网最全：基于Matlab的多维小波相干（MWC)、偏小波（PWC)与全局相干性

---

前言

在前面我们提到了两个变量的交叉小波分析，一个自变量和一个因变量，但在实际中，某一变量往往受多个其它变量的影响。我们往往关注在不同尺度与时间/空间位置上，因变量如何受多个变量的影响。虽然一些多变量方法例如多变量谱一致性（multiple spectral coherence）、多变量经验模态分解（multivariate empirical mode decomposition）可以分析多变量之间在不同尺度上的相关性，但是这些方法都假设相关性特征不随时间/空间位置变化。另外，与传统相关性方法一样，当其它变量与自变量相关时，其它变量的存在可能会导致两变量小波相干性结果具有误导性，从而不能揭示内在过程的真正机理。
在信号处理和频域分析中，小波变换是一种有效的工具，能够揭示时间序列数据的频率特征和变化规律。小波相干（Wavelet Coherence，WTC）及其衍生的偏小波相干（Partial Wavelet Coherence，PWC）、多小波相干（Multiple Wavelet Coherence，MWC）以及全局小波相干（Global Coherence）等概念，提供了更深入的频域相关性分析，有助于理解信号之间的关联程度及其时间-频率特性。本文将对这些概念进行系统解析，并介绍在MATLAB中如何实现这些分析。
小波周期图在分析非平稳信号上是十分有用的。它将时间序列信号分解到时频空间(Time-Frequency Space)以提取周期信号，水平轴对应原来时间序列的时间，垂直轴代表变化的周期，颜色代表变化周期的强度。
小波相干WTC在时间-频率空间中找到两个时间序列共同变化的区域。交叉小波变换XWT在时间-频率空间中找到时间序列里周期性强度一致的区域。两者很好地反映两个不同时间序列变化之间的“相关性”。
偏相干分析可以作为一种重要的统计方法来分析不同变量之间的相互关系，通过对数据处理和模型建立的不断优化，可以在各种领域中得到广泛的应用。
多元小波相干用来确定多个自变量对一个因变量的影响，最终变量之间的关系可以通过 MWC相干值和显著性检验的面积百分比(PASC)两个指标去评估。

偏小波相干（PWC）

偏小波相干是对小波相干的一种改进，它在分析两个信号之间的相干性时，考虑了第三个信号（或多个信号）的影响。PWC能够剔除第三个信号对两个信号之间关系的干扰，更加准确地反映出两个信号之间的真实相关性。
在MATLAB中，计算PWC的方法与WTC类似，只是需要在计算小波相干度时，加入第三个信号的信息。基本公式如下：

多小波相干（MWC）

多小波相干是对多个信号之间相干性的分析方法，它考虑了多个信号之间的交互关系，能够揭示复杂系统中信号之间的频域耦合情况。MWC的计算与PWC类似，只是需要将多个信号的信息纳入到小波变换和相干度计算中。

全局小波相干（Global Coherence）

全局小波相干是对整个信号集合的相干性进行综合分析的方法。它不仅考虑了信号之间的相互关系，还考虑了整体信号集合的频域特征和共同振荡模式。全局小波相干性分析能够帮助我们理解复杂系统中信号之间的整体协调性和耦合程度。

部分代码

function varargout=pwc(X,varargin)

%---------parse function arguments------------------------------------------
 
[row,col]=size(X);
[y1,dt1]=formatts(X(:,1));
mm=y1(1,1);
nn=y1(end,1);
 
[y2,dt2]=formatts(X(:,2));
mm2=y2(1,1);
nn2=y2(end,1);
 
for i=3:col;
[x,dtx]=formatts(X(:,i));
 
if (dt1~=dtx |dt2~=dtx)
    error('timestep must be equal between time series');
end
 
mm1=x(1,1);
nn1=x(end,1);
 
 
mm=max([mm,mm1,mm2]);
nn=min([nn,nn1,nn2]);
 
x1(:,(i-2))=x(:,1);
x2(:,(i-2))=x(:,2);
 
end
 
t=(mm:dt1:nn)';
 
 
%common time period
if length(t)<4
    error('The three time series must overlap.');
end 
  
n=length(t);
 
%----------default arguments for the wavelet transform-----------
Args=struct('Pad',1,...      % pad the time series with zeroes (recommended)
            'Dj',1/12, ...    % this will do 12 sub-octaves per octave
            'S0',1*dt1,...    % this says start at a scale of 2 years
            'J1',[],...
            'Mother','Morlet', ...
            'MaxScale',[],...   %a more simple way to specify J1
            'MakeFigure',(nargout==0),...
            'MonteCarloCount',300,...
            'BlackandWhite',0,...
            'AR1','auto',...
            'ArrowDensity',[30 30],...
            'ArrowSize',1,...
            'ArrowHeadSize',1);
 
Args=parseArgs(varargin,Args,{'BlackandWhite'});
 
if isempty(Args.J1)
    if isempty(Args.MaxScale)
        Args.MaxScale=(n*.17)*2*dt1; %auto maxscale;
    end
    Args.J1=round(log2(Args.MaxScale/Args.S0)/Args.Dj);
end
 
ad=mean(Args.ArrowDensity);
Args.ArrowSize=Args.ArrowSize*30*.03/ad;
%Args.ArrowHeadSize=Args.ArrowHeadSize*Args.ArrowSize*220;
Args.ArrowHeadSize=Args.ArrowHeadSize*120/ad;
 
if ~strcmpi(Args.Mother,'morlet')
    warning('MWC:InappropriateSmoothingOperator','Smoothing operator is designed for morlet wavelet.');
end
 
if strcmpi(Args.AR1,'auto')
      for i=1:col
arc(i)= ar1nv(X(:,i));
       end
Args.AR1=arc
    if any(isnan(Args.AR1))
        error('Automatic AR1 estimation failed. Specify it manually (use arcov or arburg).');
    end
end
 
%-----------:::::::::::::--------- ANALYZE ----------::::::::::::------------
 
%Calculate and smooth wavelet spectrum Y1, Y2, and X
  
[Y1,period,scale,coiy1] = wavelet(y1(:,2),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
sinv=1./(scale');
smY1=smoothwavelet(sinv(:,ones(1,n)).*(abs(Y1).^2),dt1,period,Args.Dj,scale);
dte=dt1*.01; 
idx=find((y1(:,1)>=(t(1)-dte))&(y1(:,1)<=(t(end)+dte)));
Y1=Y1(:,idx);
smY1=smY1(:,idx);
coiy1=coiy1(idx);
coi1=coiy1;
 
 
[Y2,period,scale,coiy2] = wavelet(y2(:,2),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
sinv=1./(scale');
smY2=smoothwavelet(sinv(:,ones(1,n)).*(abs(Y2).^2),dt1,period,Args.Dj,scale);
dte=dt2*.01; 
idx=find((y2(:,1)>=(t(1)-dte))&(y2(:,1)<=(t(end)+dte)));
Y2=Y2(:,idx);
smY2=smY2(:,idx);
coiy2=coiy2(idx);
coi2=coiy2;
 
 
 
for  i=3:col
 [XS,period,scale,coix] = wavelet(x2(:,(i-2)),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
 
idx=find((x1(:,(i-2))>=(t(1))-dte)&(x1(:,(i-2))<=(t(end)+dte)));
XS=XS(:,idx);
coix=coix(idx);
 
XS1(:,:,(i-2))=XS;
coi0=min(coi1,coi2);
coi=min(coi0,coix);
end
 
% -------- Calculate Cross Wavelet Spectra----------------------------
 
% -- between dependent variable (or independent variables) and excluding
% factors------
 
for i=1:(col-2)
Wy1x=Y1.*conj(XS1(:,:,i));
sWy1x=smoothwavelet(sinv(:,ones(1,n)).*Wy1x,dt1,period,Args.Dj,scale);
sWy1x1(:,:,i)=sWy1x;
 
Wy2x=Y2.*conj(XS1(:,:,i));
sWy2x=smoothwavelet(sinv(:,ones(1,n)).*Wy2x,dt1,period,Args.Dj,scale);
sWy2x1(:,:,i)=sWy2x;
end
 
% ---- between dependent variable and independent variables------
Wy1y2=Y1.*conj(Y2);
sWy1y2=smoothwavelet(sinv(:,ones(1,n)).*Wy1y2,dt1,period,Args.Dj,scale);
 
% ----between excluding variables and excluding variables------
for i=1:(col-2);
for j=1:(col-2);
Wxx=XS1(:,:,i).*conj(XS1(:,:,j));
sWxx=smoothwavelet(sinv(:,ones(1,n)).*Wxx,dt1,period,Args.Dj,scale);
sWxx1(:,:,i,j)=sWxx;
end
end
 
% --------------- Partial wavelet coherence ---------------------------------
% calculate the partial wavelet coherence 
m=length(scale);
Rsq=zeros(m,n);
Wy1y2x=zeros(m,n);
nofe=col-2; %number of excluding factors
 
if nofe==1;    
% ------Partial wavelet coherence for one excluding variable ------------------
%%%%%%%%% |1-R^2 y,x穁 (s,t)|^2   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:m;
    for jj=1:n;
        sWxx1_i(ii,jj)=inv(sWxx1(ii,jj)); % inversed sWxx1
    end
end
 
sWy1x1_sWxx1=bsxfun(@times,sWy1x1,sWxx1_i); % sWy1x1*sWxx1_i
sWy2x1_c=conj(sWy2x1); %conjugate of smoothed cross-wavelet power spectra
% between independent variable and excluding variables
sWy1x1_sWxx1_sWy2x1=bsxfun(@times,sWy1x1_sWxx1,sWy2x1_c); 
% sWy1x1*sWxx1_i*sWy2x1_c
R2y1y2x=sWy1x1_sWxx1_sWy2x1./sWy1y2; %(Ry,x.Z)^2
abs_one_minus_R2y1y2x=(abs(1-R2y1y2x)).^2; %squared absolute of 1-(Ry,x.Z)^2
%%%%%%%%%%%%%%%% R^2 y,x (s,t)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R2y1y2=(sWy1y2.*conj(sWy1y2))./(smY1.*smY2); %(Ry,x)^2
%%%%%%%%%%%%%%%%%%% 1-R^2 y,Z(s,t)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sWy1x1_c=conj(sWy1x1);
sWy1x1_sWxx1_sWy1x1=bsxfun(@times,sWy1x1_sWxx1,sWy1x1_c);
R2y1x=sWy1x1_sWxx1_sWy1x1./smY1; %(Ry,Z)^2
one_minus_R2y1x=(1-R2y1x);
 
%%%%%%%%%%%%%%%%%%%%%  1-R^2 x,Z(s,t) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sWy2x1_sWxx1=bsxfun(@times,sWy2x1,sWxx1_i);
sWy2x1_sWxx1_sWy2x1=bsxfun(@times,sWy2x1_sWxx1,sWy2x1_c);
R2y2x=sWy2x1_sWxx1_sWy2x1./smY2; %(Rx,Z)^2
one_minus_R2y2x=(1-R2y2x);
%%%%%%%%%%%%%%%%%%  R^2 y,x穁 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rsq_num=roundn(bsxfun(@times,abs_one_minus_R2y1y2x,R2y1y2),-16); 
%numerator part of the equation for squared pwc, if the value is less than 
%10^(-16), 0 is assigned
Rsq_den=real(bsxfun(@times,one_minus_R2y1x,one_minus_R2y2x));  
%denominator part of the equation for squared pwc
Rsq= bsxfun(@rdivide,Rsq_num,Rsq_den); %squared partial wavelet coherence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

出图效果

以下是几篇论文中的出图效果。

土壤湿度的小波相干性结果说明三因子组合被认为是解释每个研究点土壤湿度变化的最具代表性的。上排子图为表土层(0-10 cm)的分析，第二行对应第二层土层(10-40 cm)，第三行对应第三层土层(40-100 cm)，下行为第四层土层(100-200 cm)。淡色区域是影响锥(COI)，其中边缘效应可能会扭曲结果。粗轮廓对红色噪声的显著性水平为5%，颜色表示一致性的强度。
下图为代码出来的效果图。小波交叉谱:U型黑色细实线为小波边界效应影响锥(COI)，在该曲线以外的功率谱由于受到边界效应的影响而不子考虑; 黑色粗黑线包围的范围通过了0.05显著性水平下的红噪声标准谱的检验，箭头代表相位关系， 其中←表示两个时序变量的变化位相相反，即两者呈负相关关系，→ 表示两者变化位相致，即两者呈正相关关系，↑表示指数A落后因子B变化1/4周期，↓表示指数A超前因子B变化1/4周期; 色柱表示能量密度的相对变化; XWT：交叉小波谱，WTC小波相干谱，以此类推。

欢迎扫描关注我的公众号，博士期间科研分享日常，获取更多科研代码和前沿论文资讯等相关内容
若需要以上小波分析的代码，请关注公众号回复“小波分析”关键词

参考文献

1.Zhao D, Xiong D, Zhang B, et al. Long-term response of runoff and sediment load to spatiotemporally varied rainfall in the Lhasa River basin, Tibetan Plateau[J]. Journal of Hydrology, 2023, 618: 129154.
2.He, K., Shi, H., Chen, C., Cheng, Y., & Liu, J. (2021). The study on the time lag of water level in the Three Gorges Reservoir under the regulation processes. Hydrology Research, 52, 734-748
3.Hu, W., and B.C. Si (2021), Technical Note: Improved partial wavelet coherency for understanding scale-specific and localized bivariate relationships in geosciences, Hydrol. Earth Syst. Sci.,25,321-331.
4.Gu, X., Jamshidi, S., Sun, H., & Niyogi, D. (2021). Identifying multivariate controls of soil moisture variations using multiple wavelet coherence in the U.S. Midwest. Journal of Hydrology, 602, 126755
5.Zhao, R., Biswas, A., Zhou, Y., Zhou, Y., Shi, Z., & Li, H. (2018). Identifying localized and scale-specific multivariate controls of soil organic matter variations using multiple wavelet coherence. Science of the total environment, 643, 548-558