#======================================================================
# L A D C P P R O C . B A C K S C A T T E R
# doc: Wed Oct 20 13:02:27 2010
# dlm: Fri Mar 6 15:52:50 2015
# (c) 2010 A.M. Thurnherr
# uE-Info: 32 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Oct 20, 2010: - created
# Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
# when run with uplooker data
# Jun 15, 2011: - added calculation of backscatter profiles from
# subset of bins
# Jul 7, 2011: - use $BEAM? vars to clarify code
# - save SV values to use in BT code
# Jul 14, 2011: - implemented Sv_V04
# Jul 15, 2011: - modified Sv_V04 to take noise levels into account (Sv_T11)
# Apr 10, 2012: - BUG: ensonified volume in Sv_T11 was wrong
# May 16, 2012: - BUG: find_backscatter_seabed() used wrong start depth for
# search
# - BUG: same used bins entirely below seabed (only possible
# for shallow casts)
# Mar 21, 2014: - adapted to new [LADCPproc.utils]
# Mar 27, 2014: - adapted to depthOfBinAlongBeam()
# May 25, 2014: - made search_below integer
# Jul 27, 2014: - moved depthOfGI() to [LADCPproc.utils]
# Aug 5, 2014: - BUG: find_backscatter_seabed() discarded everything if
# LADCP bin 1 had the backscatter max at the edge of
# the search domain, which could easily happen when
# the bottom stop was a long way from the seabed
my($BEAM1) = 0;
my($BEAM2) = 1;
my($BEAM3) = 2;
my($BEAM4) = 3;
#----------------------------------------------------------------------
# Volume Scattering Coefficient
#
# 1) There are currently 4 different methods for correcting echo
# amplitudes as a function of range. The method to calculate Sv
# of Deines (1999) is the only published one, although Visbeck's
# correction was presumably used by Visbeck and Thurnherr (2009).
# The third method (Thurnherr 2011) improves on Visbeck's code.
# The fourth method consists in no correction.
# 2) The following observations stand out when plotting echo amplitudes
# against bin number in clear water (e.g. below 1400m in first
# downcast of IWISE yoyo cast 160; [Documentation/Sv_comparison.eps]):
# - the corrections works approximately for the valid bins
# (2-13 in case of IWISE) but it amounts to a slight
# overcorrection, i.e. the resulting corrected backscatter
# profiles increase with increasing distance from the
# transducer
# - in terms of remaining vertical structure, the different
# methods are nearly identical
# - however, Visbeck (2004) shows fairly large instrument
# dependent biases because (instrument-specific) noise
# is not considered
# - the vertical structure can be improved by decreasing the
# attenuation parameter to something like 0.05; I have
# verified this for the Deines and Thurnherr expressions
# with a single IWISE station; however, there remain
# small (but noticeable) instrument-specific biases
# 3) Therefore, the default code is Deines (1999); if this is
# not good enough, it is probably possible to derive an
# empirical correction from in-situ data.
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
# NOTES:
# - instrument specific! (300kHz Workhorse)
# - no sound-speed correction applied
# - R in bin center, instead of last quarter
# - transmit power assumes 33V batteries
# Results:
# - correction seems to work for a subset of bins (~bins 3-9 for
# 2010 P403 station 46)
# - this may imply that noise level depends on bin
# - far bins are important for seabed detection, i.e. cannot simply
# be discarded at this stage
#----------------------------------------------------------------------
sub log10 {
my $n = shift;
return log($n)/log(10);
}
sub Sv_D99($$$$$)
{
my($temp,$PL,$Er,$R,$EA) = @_;
my($C) = -143; # RDI Workhorse monitor
my($Ldbm) = 10 * log10($PL);
my($PdbW) = 14.0;
my($alpha) = 0.069;
# my($alpha) = 0.048;
my($Kc) = 0.45;
return $C + 10*log10(($temp+273)*$R**2) - $Ldbm - $PdbW
+ 2*$alpha*$R + $Kc*($EA-$Er);
}
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Visbeck (code 2004)
#
## function [ts,bcs]=targ($EA,$R,$alpha,$PL,$EAS,$ap)
## Target strength of EA for volume scatterer
## $EA = echoamp in dB
## $R = distance in m
## $alpha = attenuation dB/m
## $PL = pulse/bin legth in m
## $EAS = source level
## $ap = aperature in degree
## M. Visbeck 2004
#
# Results:
# - overall, correction with distance works similarly well to Deines (1999)
# - constant bias, which could be taken care of by changing EAS
# - however, much more serious UL/DL differences than Deines (1999)
#----------------------------------------------------------------------
sub Sv_V04($$$$$)
{
my($temp,$PL,$Er,$R,$EA) = @_; # only uses $R and $EA
# my($alpha) = 0.039; ## attenuation dB/m for 150 kHz
my($alpha) = 0.062; ## attenuation dB/m for 300 kHz
my($EAS) = 100; ## source level in dB
my($ap) = rad(2); ## beam aperature in DEGREE convert to radian
my($r1) = tan($ap)*($R-$PL/2); ## radius of top and bottom of each bin
my($r2) = tan($ap)*($R+$PL/2);
my($V) = $PI*$PL/3 * ($r1**2+$r2**2+$r1*$r2); ## ensonified volume
my($TL) = 20*log10($R) + $alpha*$R; ## transmission loss
my($TS) = 0.45*$EA - $EAS + 2*$TL - 10*log10($V); ## target strength
my($BCS) = exp($TS/10); ## backscatter cross section
return $TS;
}
#----------------------------------------------------------------------
# Volume Scattering Coefficient, Thurnherr [2011]
#
# in an attempt to improved on Visbeck's code, I made several
# modifications:
# - add 5cm to beam radius (transducer radius)
# - take beam noise into account (as in Deines [1999])
# - use different equation for ensonified volume
# results:
# - inclusion of beam noise greatly reduces UL/DL differences in IWISE data
# - constant offset from Deines values, which could be taken care of by changing EAS
# - remaining changes don't have great effects
#----------------------------------------------------------------------
sub Sv_T11($$$$$)
{
my($temp,$PL,$Er,$R,$EA) = @_;
# my($alpha) = 0.067; # Deines [1999]
my($alpha) = 0.062; # Visbeck [2004]
# my($alpha) = 0.048; # attenuation dB/m for 300 kHz
my($EAS) = 100; # source level in dB
my($ap) = rad(2); # beam aperature in DEGREE convert to radian
my($r0) = 0.05; # beam radius at source (transducer radius)
my($d1) = $R - $PL/2; # distance to top and bottom of bin
my($d2) = $R + $PL/2;
my($r1) = $r0 + $d1*tan($ap); # radius of top and bottom of each bin
my($r2) = $r0 + $d2*tan($ap);
my($V) = $PI/3 * ($d2*$r2**2 - $d1*$r1**2); # ensonified volume
my($TL) = 20*log10($R) + $alpha*$R; # transmission loss
my($TS) = 0.45*($EA-$Er) - $EAS + 2*$TL - 10*log10($V); # target strength
return $TS;
}
#----------------------------------------------------------------------
# no correction for attenuation losses
#
# however, noise level is still taken into consideration
#----------------------------------------------------------------------
sub Sv_nocorr($$$$$)
{
my($temp,$PL,$Er,$R,$EA) = @_;
return 0.45*($EA-$Er);
}
#======================================================================
sub mk_backscatter_profs($$)
{
my($LADCP_start,$LADCP_end) = @_;
my(@Er) = (1e99,1e99,1e99,1e99); # echo intensity reference level
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
$Er[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1]
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1] < $Er[$BEAM1]);
$Er[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2]
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2] < $Er[$BEAM2]);
$Er[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3]
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3] < $Er[$BEAM3]);
$Er[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4]
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4] < $Er[$BEAM4]);
}
print(STDERR "\n\t\@per-beam noise levels = @Er") if ($opt_d);
my($Svfunc); # Sv method
if ($opt_u =~ /^[dD]/) {
$Svfunc = \&Sv_D99;
} elsif ($opt_u =~ /^[vV]/) {
$Svfunc = \&Sv_V04;
} elsif ($opt_u =~ /^[tT]/) {
$Svfunc = \&Sv_T11;
} else {
$Svfunc = \&Sv_nocorr;
}
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
my($range_to_bin) = &rangeToBin($ens,$bin);
for (my($beam)=0; $beam<4; $beam++) {
my($gi) = int(&depthOfBinAlongBeam($ens,$bin,$beam) / $GRID_DZ);
next if ($gi < 0);
$LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam] = &$Svfunc($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
$LADCP{TRANSMITTED_PULSE_LENGTH},
$Er[$beam],$range_to_bin,
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$beam]);
$sSv[$gi][$bin] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam];
$nSv[$gi][$bin]++;
if ($bin>=$Svbin_start && $bin<=$Svbin_end) {
$sSv_prof[$gi] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam];
$nSv_prof[$gi]++;
} # if $bin
} # for $beam
} # for $bin
} # for $end
} # sub
sub find_backscatter_seabed($)
{
my($water_depth) = @_;
my(@wdepth_gi); # water_depth indices
my($search_below) = int(max(0,$water_depth-$BT_begin_search_above));
my($mdgi) = int($search_below/$GRID_DZ); # grid index to begin search
printf(STDERR "\n\t\tlooking for seabed below %d m (gi = [%d..%d])",$search_below,$mdgi,scalar(@nSv))
if ($opt_d);
print(STDERR "\n\t\tseabed-max grid indices:") if ($opt_d);
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { # find backscatter min/max below $search_below in each bin
my($min,$max,$gimax,$firstvalid,$lastvalid) = (1e99,-1e99,-1,-1,-1);
for (my($gi)=$mdgi; $gi<@nSv; $gi++) {
next unless ($nSv[$gi][$bin] > 0);
$firstvalid = $gi if ($firstvalid < 0);
my($avg) = $sSv[$gi][$bin] / $nSv[$gi][$bin];
$lastvalid = $gi;
$min = $avg if ($avg < $min);
$max = $avg, $gimax = $gi if ($avg > $max);
}
if ($max-$min>10 && $gimax!=$firstvalid && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
printf(STDERR " %d",$gimax) if ($opt_d);
push(@wdepth_gi,$gimax);
}
}
return (depthOfGI(avg(@wdepth_gi)),stddev(@wdepth_gi)*$GRID_DZ);
}
1;