MATLAB在地下流體數據分析中的基礎應用(Ⅱ)

何案華 王言章

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

井水位觀測值動態表達式一般為：

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

(1)

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

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

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

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

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

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

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

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

根據上述過程，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('水位原始觀測數據','水位低通濾波后數據','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('去掉低頻波成分后水位數據','location','best');

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

將水位原始數據通過取均值方法[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，負號是由于水位為靜水位數據，其數值大小與重力潮汐呈負相關(圖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);%計算水位與氣壓最大相關系數，確定氣壓效應滯后時間

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');

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

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

2 重力潮汐理論值計算

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

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 結 語

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

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