何案華 王言章

1 吉林大學地球信息探測儀器教育部重點實驗室，長春市民主大街938號，130026 2 吉林大學儀器科學與電氣工程學院，長春市民主大街938號，130026 3 應急管理部國家自然災害防治研究院，北京市安寧莊路1號，100085

井水位觀測值動態(tài)表達式一般為：

yn=xn+Pn+Tn+Rn+εn

(1)

式中，yn為水位觀測值；xn為水位殘差；Pn為氣壓效應；Tn為潮汐效應(體應變潮汐理論值)；Rn為降雨效應；εn為測量噪聲，假定其為均值為0、方差為σ2的高斯白噪聲[1]。為簡化計算，暫不考慮井水位降雨效應。從已有研究來看，井水位中潮汐效應與氣壓效應并非是簡單的一元線性關系，而是存在延時、滯后且持續(xù)一段時間。考慮實際過程，本文側重討論一元線性關系的計算，并對一元線性擬合后的數(shù)學問題進行簡單討論。

1 井水位潮汐因子與氣壓因子計算

井水位潮汐與氣壓效應計算過程如圖1所示，主要步驟如下：

圖1 井水位潮汐與氣壓效應計算過程Fig.1 The calculating process of tidal and barometric effects of groundwater level

1)源數(shù)據預處理。①剔除明顯錯誤數(shù)據；②缺數(shù)處理，由于MATLAB許多計算過程要求數(shù)據完整，故采用一定數(shù)學方法對缺數(shù)進行插值；③降采樣率，由于觀測數(shù)據采樣率多為1次/min，甚至1次/s，而地球固體潮汐、氣壓等周期都是以小時為單位，為減少計算的復雜性，進行適當降采樣率處理，本文采用1次/10 min。

2)潮汐理論值計算。常規(guī)計算方法為采用專用軟件，如MAPSIS按指定要求進行計算，得到相應數(shù)據文件后再參與計算，但這會增加處理的復雜性。通過梳理重力固體潮汐理論，采用H_gravity_earthtide自定義函數(shù)進行實時計算，即根據臺站基礎信息(經緯度)以及數(shù)據時間進行計算。

3)潮汐因子計算。采用低通濾波提取趨勢變化，然后用源數(shù)據減去低通濾波的趨勢性變化即得到僅保留潮汐成分的水位，利用一元線性回歸計算潮汐因子，進而剔除水位中的潮汐成分。

4)氣壓因子計算。由于氣壓效應存在滯后效應[2]，故先求出最大相關系數(shù)，并取得滯后時間；再根據滯后時間進行一元線性回歸計算，求取氣壓因子并剔除氣壓干擾。

根據上述過程，MATLAB代碼如下：

fc=1/(24*60*60*2);%截止頻率

fs=1/(10*60);%采樣頻率

[fdata,ddata]=H_butter_lowpass( fc,fs,sw );

subplot(211);

plot(dates,sw,'b.',dates,fdata,'r--','linewidth',2);

set(gca,'xlim',[xmin xmax],'xtick',xtk,'xticklabel',xtkl,'fontsize',12,'fontweight','bold','ydir','reverse');

xlabel('日期(2022年)','fontsize',12,'fontweight','bold');

ylabel('水位(m)','fontsize',12,'fontweight','bold');

legend('水位原始觀測數(shù)據','水位低通濾波后數(shù)據','location','best');

subplot(212);

plot(dates,ddata,'linewidth',2);

set(gca,'xlim',[xmin xmax],'xtick',xtk,'xticklabel',xtkl,'fontsize',12,'fontweight','bold','ydir','reverse');

xlabel('日期(2022年)','fontsize',12,'fontweight','bold');

legend('去掉低頻波成分后水位數(shù)據','location','best');

ylabel('水位(m)','fontsize',12,'fontweight','bold');

將水位原始數(shù)據通過取均值方法[3]降采樣率為1次/10 min，低通濾波取截止頻率fc=1/(24·60·60·2)，圖2為低通濾波后產品。

圖2 井水位低通濾波及殘差Fig.2 Low pass filtering and residual error of groundwater level

td =[];%按臺站信息求出重力潮汐理論值

for i=1:length(dates)

td=[td H_gravity_earthtide(dates(i),stationlongitude,stationlatitude)];

end

figure(2);

subplot(221);

plottwo(dates,ddata,td,'reverse','normal','水位(m)','重力潮汐(uGal)','日期(2022年)',

datenum(0,0,5),'mm/dd');

subplot(222);

plot(td,ddata,'*'); hold on;

x=[td' ones(1,length(td))'];

y=ddata';

[b,bint,r,rint,stats]=regress(y,x,0.05);

x=-170:10:100;

y=b(1).*x+b(2);

plot(x,y,'r--','linewidth',2); hold off;

set(gca,'xlim',[-170 100],'xtick',-170:20:100);

xlabel('重力潮汐(uGal)','fontsize',12,'fontweight','bold');

ylabel('水位(m)','fontsize',12,'fontweight','bold');

strtemp=sprintf('y=%.06f*x+%.04f',b(1),b(2));

text(-70,0.1,strtemp,'fontsize',12,'fontweight','bold');

采用自定義的H_gravity_earthtide計算重力潮汐理論值。通過MATLAB自帶的regress完成一元線性回歸，可通過stats查詢回歸結果，在本實例中得到回歸結果R2=0.930 2，說明回歸效果非常理想。求得潮汐因子為-0.000 479 m/μGal，負號是由于水位為靜水位數(shù)據，其數(shù)值大小與重力潮汐呈負相關(圖3)。

圖3 井水位潮汐因子計算Fig.3 The tidal factor calculating of groundwater level

figure(3);

sw_eliminate_td=detrend(sw-(b(1).*td+b(2)));%剔除潮汐干擾

qytemp=qy(1:144*3);%計算水位與氣壓最大相關系數(shù)，確定氣壓效應滯后時間

cor=[];

for i=1:1:20

swtemp=sw_eliminate_td(i:i+144*3-1);

cor=[cor min(min(corrcoef(swtemp,qytemp)))];

end

[val pos]=max(cor);

subplot(221);

plottwo(dates(pos:end),sw_eliminate_td(pos:end),qy(1:end-pos+1),'reverse','reverse','水位(m)','氣壓(hPa)','日期(2022年)',datenum(0,0,5),'mm/dd');

swtemp=sw_eliminate_td(pos:end);

qytemp=qy(1:end-pos+1);

x=[qytemp' ones(1,length(qytemp))'];

y=swtemp';

[b,bint,r,rint,stats]=regress(y,x,0.05);

subplot(222);

plot(qytemp,swtemp,'*'); hold on;

x=1008:1027;

y=b(1)*x+b(2);

plot(x,y,'r--','linewidth',2); hold off;

xlabel('氣壓(hPa)','fontsize',12,'fontweight','bold');

ylabel('水位(m)','fontsize',12,'fontweight','bold');

strtemp=sprintf('y=%.06f*x+%.04f',b(1),b(2));

text(1016,-0.04,strtemp,'fontsize',12,'fontweight','bold');

x=qytemp;

y=swtemp-(b(1)*qytemp+b(2));

subplot(212);

plot(dates(pos:end),y,'linewidth',2);

[xmin,xmax,xtk,xtkl]=GetXInfo(dates(pos:end),datenum(0,0,5),'mm/dd');

set(gca,'xlim',[xmin xmax],'xtick',xtk,'xticklabel',xtkl,'fontsize',12,'fontweight','bold','ydir','reverse');

text(datenum(2022,1,15,19,30,0),-0.01,'downarrow湯加火山','color','r','fontsize',12,'fontweight','bold');

采用查找最大相關系數(shù)方法，先求取氣壓效應的滯后時間，本例中氣壓效應滯后時間pos=6，采樣率為1次/10 min，即氣壓效應滯后時間為60 min。利用相同線性回歸求得氣壓因子為0.004 631 m/hPa。通過去除氣壓效應，可清晰觀察到原來隱藏在觀測數(shù)據中的湯加火山爆發(fā)事件(圖4)。

圖4 井水位氣壓因子計算及最終井水位動態(tài)Fig.4 The barometric factor calculating of groundwater level and the final groundwater level dynamics

2 重力潮汐理論值計算

重力潮汐理論值的計算原理與計算過程可參考文獻[4-5]，但上述文獻未給出詳細的MATLAB實現(xiàn)，本文在此提供。

function Gg=H_gravity_earthtide(date,longtitude,latitude)

longtitude=longtitude * pi/180;

latitude=latitude * pi/180;

Phi=latitude;

Phip=atan(power(1-1/298.256,2)*tan(Phi));

[y,m,d,HH,MM,SS]=datevec(date);

t=HH+MM/60+SS/3600;

T0=juliandate(y,m,d,0,0,0);

T=(T0-2451545.0+(t-8)/24)/36525;

S=218.31643+481267.88128*T-0.00161*power(T,2)+0.000005 * power(T,3);

H=280.46607+36000.76980*T+0.00030*power(T,2);

P=83.35345+4069.01388*T-0.01031*power(T,2)-0.000001*power(T,3);

N=125.04452-1934.13626*T+0.00207*power(T,2)+0.000002*power(T,3);

Ps=282.93835+1.71946*T+0.00046*power(T,2)+0.000003*power(T,3);

sigma=23.43929-0.01300*T-0.000000167*power(T,2)+0.0000005*power(T,3);

S=S * pi/180;

H=H * pi/180;

P=P * pi/180;

N=N * pi/180;

Ps=Ps * pi/180;

sigma=sigma * pi/180;

Cm_Rm=1+0.01*cos(S-2*H+P)+0.0545*cos(S-P)+0.0030*cos(2*S-2*P)+0.0009*cos(3*S-2*H-P)...

+0.0006*cos(2*S-3*H+Ps)+0.0082*cos(2*S-2*H);

Lambdam=S-0.00324*sin(H-Ps)-0.00103*sin(2*H-2*P)+0.001*sin(S-3*H+P+Ps)+0.02224*sin(S-2*H+P)...

+0.00072*sin(S-H-P+Ps)-0.00061*sin(S-H)+0.10976*sin(S-P)-0.00053*sin(S+H-P-Ps)...

+0.0008*sin(2*S-3*H+Ps)+0.01149*sin(2*S-2*H)+0.00373*sin(2*S-2*P)-0.002*sin(2*S-2*N)+0.00093*sin(3*S-2*H-P);

Beltam=-0.00485*sin(P-N) -0.00081*sin(2*H-P-N)+0.00303*sin(S-2*H+N)+0.0895*sin(S-N)+0.00097*sin(2*S-2*H+P-N)...

+0.0049*sin(2*S-P-N)+0.00057*sin(3*S-2*H-N);

Cs_Rs=1+0.0168*cos(H-Ps)+0.0003*cos(2*H-2*Ps);

Lambdas=H+0.0335*sin(H-Ps)+0.0004*sin(2*H-2*Ps);

Theta =((t-8)*15*pi/180)+H+longtitude-pi;

sinDeltam=sin(sigma)*sin(Lambdam)*cos(Beltam)+cos(sigma)*sin(Beltam);

cosDeltamcosTaum=cos(Beltam)*cos(Lambdam)*cos(Theta)...

+sin(Theta)*(cos(sigma)*cos(Beltam)*sin(Lambdam)-sin(sigma)*sin(Beltam));

sinDeltas=sin(sigma)*sin(Lambdas);

cosDeltascosTaus=cos(Lambdas)*cos(Theta)+sin(Theta)*cos(sigma)*sin(Lambdas);

cosZm=sin(Phip)*sinDeltam+cos(Phip)*cosDeltamcosTaum;

cosZs=sin(Phip)*sinDeltas+cos(Phip)*cosDeltascosTaus;

F_Phi=0.998327+0.00167*cos(2*Phi);

Gg=-165.17*F_Phi* power(Cm_Rm,3) *(power(cosZm,2)-1.0/3.0)...

- 1.3708 * F_Phi * F_Phi * power(Cm_Rm,4) *(5*power(cosZm,3)-3*cosZm)...

-76.08*F_Phi*power(Cs_Rs,3)*(power(cosZs,2)-1.0/3.0);

end

3 結 語

通過上述過程，可自動計算井水位潮汐因子、氣壓因子。通過剔除井水位潮汐干擾與氣壓干擾，可清晰地觀測到隱藏在井水位原始觀測數(shù)據中由湯加火山爆發(fā)引發(fā)的氣壓擾動事件。

從圖4最終產品可知，井水位動態(tài)中仍然可清晰地看到潮汐成分，這與井水位潮汐和氣壓響應原理有關，兩者之間并非簡單的線性關系，存在時間滯后，且響應不是一個時間點，而是一個作用過程。另外，氣壓動態(tài)自身包含潮汐作用力，因此將潮汐與氣壓分兩步進行剔除會引入氣壓分量中的潮汐作用力，這在上述計算過程中無法避免。要解決該問題，需采用多元線性、偏最小二乘法、卷積等其他回歸算法，在此不過多討論，詳細過程可參考文獻[2]。