programme de calcul du backscatter issu des données ADCP

backscatter/target_strength.m

+function out = target_strength(EA,EA0,S,T,sspd,tempx,sspdx,binlength,blank,beamangle,freq)
+% function out = targetstrength(EA,EA0,S,T,sspd,tempx,sspdx,binlength,blank,beamangle)
+%
+% freq : khz
+% EA:  echo amplitudes in counts;         [d x t]
+% EA0: noisefloor from ADCP electronic noise; a scalar
+% S:   salinity on ADCP-grid;             [d x t]
+% T:   in-situ temperature on ADCP-grid;  [d x t]
+% sspd: soundspeed on ADCP-grid;          [d x t]
+% tempx: transducer temerature e.g. from Thermosal or ADCP binary data; a vector [1 x t]
+% sspdx: transducer soundspeed e.g. from Thermosal or ADCP binary data; a vector [1 x t]
+% binlength: ADCP configurated length of bin; a scalar
+% blank:     ADCP configurated blank length; a scalar
+% beamangle: ADCP beam angle to vertical; a scalar
+%
+% units: EA, EA0:     counts (rawdata-output)
+%        T, tempx:    degrees Celsius
+%        sspd, sspdx: m/s
+%        beamangle:   degrees
+%        bin, blank:  m
+%
+% out is a structure with
+% out.ts: target strength in dB (relative)
+% out.csspd: 'cumulative' soundspeed = effective sspd for this bin
+% out.R: slant range along beam
+% out.alpha: absorption coefficient = average of water column until this bin
+%
+% calculations following Plimpton, Freitag, McPhaden
+% 'Processing of subsurface ADCP data in the equatorial Pacific'
+%
+% last changed 11/02/18 TF
+
+pH = 8.1;  % seawater pH
+
+szm = size(EA);
+nbins = szm(1);
+
+dist = ([1:nbins]'-0.5)*binlength+blank;
+si = [dist(1);diff(dist)];
+tges = cumsum(repmat(si,1,szm(2))./sspd);
+csspd = repmat(dist,1,szm(2))./tges;
+out.csspd = csspd;
+
+Rconst = (blank+0.31/2+([1:nbins]'-0.5)*binlength+binlength/4)/cos(beamangle*pi/180);
+R = repmat(Rconst,1,szm(2)).*csspd./repmat(sspdx,nbins,1);
+out.R = R;
+
+%% Calculation of sound absorption (Appendix 2 Plimpton et al 2004)
+
+% Pure water contribution:
+A3 = 4.937e-4-2.59e-5*T+9.11e-7*T.^2-1.5e-8*T.^3;
+A31 = 3.964e-4-1.146e-5*T+1.45e-7*T.^2-6.5e-10*T.^3;
+A3(T>20) = A31(T>20);
+zz = R*cos(beamangle*pi/180);
+P3 = 1-3.83e-5*zz+4.9e-10*zz.^2;
+% MgSO4 contribution:
+A2 = 21.44*S.*(1+0.025*T)./sspd;
+P2 = 1-1.37e-4*zz+6.2e-9*zz.^2;
+f2 = 8.17*(10.^(8-1990./(273.16+T)))./(1+0.0018*(S-35));
+% Boric acid contribution:
+A1 = 8.86*10^(0.78*pH-5)./sspd;
+f1 = 2.8*sqrt(S/35).*10.^(4-1245./(273.16+T));
+
+alpha = (A3.*P3*freq^2+A2.*P2.*f2.*freq^2./(freq^2+f2.^2)+A1.*f1.*freq^2./(freq^2+f1.^2))/1000; % dB/m
+alpham = cumsum(alpha)./repmat([1:nbins]',1,szm(2));
+out.alpha = alpham;
+
+ampcor = EA-EA0;
+ampcor(ampcor<0) = 0;
+out.ts = 10*log10((repmat(tempx,nbins,1)+273.16)./sspd)+10*log10(10.^(127.3./(repmat(tempx,nbins,1)+273.16).*ampcor/10)-1)+2*alpham.*R+20*log10(R);

backscatter/template_get_adcp_TS.m

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% template_get_adcp_TS.m
+% -------------------------------
+% Author : Jrmie HABASQUE - IRD
+% -------------------------------
+% INPUTS:
+% - ADCP data processed
+% - CTD file
+% OUTPUTS:
+% - 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+close all
+clear all
+
+% path
+addpath('..\moored_adcp_proc');
+
+% Directory for outputs
+fpath_output = '..\data_example\';
+
+% load ADCP data
+load('..\data_example\_15258_instr_01_int_filt_sub.mat');
+
+Z = data.Z;
+z = data.z_bins;
+time = data.time;
+freq = raw.config.sysconfig.frequency;
+blen = raw.config.cell;    % bin length
+nbin = raw.config.ncells;  % number of bins
+bin  = 1:nbin;
+
+% prepare S, T and sspd on ADCP-grid [d x t]
+% S:   salinity on ADCP-grid;             [d x t]
+% T:   in-situ temperature on ADCP-grid;  [d x t]
+% sspd: soundspeed on ADCP-grid;          [d x t]
+
+EA = data.ea(1:length(Z),:);
+EA0 = 18;
+tempx = data.temp;
+sspdx = data.sspd;
+[sspd,dummy] = meshgrid(sspdx,1:length(Z));
+binlength = adcp.config.cell;
+blank = adcp.config.blank;
+beamangle = adcp.config.sysconfig.angle;
+% mean salinity and temperature profile at deployment and recovery
+file_ctd = 'Z:\PIRATA-FR25\data-processing\CTD\OS_PIRATA-FR25_CTD.nc';
+lat = ncread(file_ctd, 'LATITUDE') ;
+lon = ncread(file_ctd, 'LONGITUDE') ;
+depth = ncread(file_ctd, 'PRES');
+temp = ncread(file_ctd, 'TEMP');
+sal = ncread(file_ctd, 'PSAL');
+TIME = ncread(file_ctd, 'TIME') ;
+TIME_CTD = TIME * 3600 *24;
+for k=1:length(TIME_CTD)
+    TIME_CTD(k)=datenum([2012 1 1 00 00 TIME_CTD(k)]);
+end
+for k=1:length(time)
+    [yr, mn, dy, hr]= gregorian(time(k));
+    time_adcp(k)=datenum(yr, mn, dy, hr, 00, 00);
+end 
+% look for the closest CTD profile
+ind_profile = find(TIME_CTD>time_adcp(1), 1, 'first' );
+sal_adcp = ones(length(Z),1);
+temp_adcp = ones(length(Z),1);
+% for each ADCP bin, find salinity and temperature values
+for i_bin=1:length(Z)    
+    ind_depth = find(depth(:,ind_profile)>Z(i_bin), 1, 'first' );
+    sal_adcp(i_bin) = sal(ind_depth,ind_profile);
+    temp_adcp(i_bin) = temp(ind_depth,ind_profile);
+end
+S = sal_adcp*ones(size(tempx,1),1)'; 
+T = temp_adcp*ones(size(tempx,1),1)'; 
+
+% Compute TS
+out = target_strength(EA,EA0,S,T,sspd,tempx',sspdx',binlength,blank,beamangle,freq);
+
+save([fpath_output,mooring.name, '_TS'],'out');
+
+%% Interpolate data on a regular vertical grid
+Z = fliplr(blen/2:blen:max(z(:))+blen);
+Zmax = max(Z);
+TS_interp = NaN(length(time),length(Z));
+
+for i=1:length(time)
+    % indice correspondant sur la grille finale Z
+    ind = round((Zmax-z(1,i))/blen)+1;
+    % filling the grid
+    npts = min([length(Z)-ind+1 length(bin)]);
+    TS_interp(i,ind:ind+npts-1) = out.ts(1:npts,i);
+end
+    
+% Save interpolated data
+save([fpath_output,mooring.name, '_TS_interp'],'out');
+ 
+%% Figure
+hf=figure('position', [0, 0, 1400, 1000]);
+[C,h] = contourf(time,Z,TS_interp',50); 
+set(h,'LineColor','none');
+h=colorbar;
+ylabel(h,'TS [dB]');
+set(gca,'ydir', 'reverse');
+ylabel('Depth (m)');
+ylim([0,adcp.instr_depth]);
+%change figure label in HH:MM
+gregtick
+title({[mooring.name, ' - Relative backscatter - RDI ',num2str(freq),' kHz']});
+
+graph_name = [fpath_output, mooring.name '_relative_backscatter'];
+set(hf,'Units','Inches');
+pos = get(hf,'Position');
+set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
+print(hf,graph_name,'-dpdf','-r300');
+
+% remove path
+rmpath('.\moored_adcp_proc');