2
|
1 |
#======================================================================
|
|
2 |
# A C O U S T I C _ B A C K S C A T T E R . P L
|
|
3 |
# doc: Wed Oct 20 13:02:27 2010
|
|
4 |
# dlm: Thu Dec 30 22:22:02 2010
|
|
5 |
# (c) 2010 A.M. Thurnherr
|
|
6 |
# uE-Info: 131 36 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# HISTORY:
|
|
10 |
# Oct 20, 2010: - created
|
|
11 |
# Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
|
|
12 |
# when run with uplooker data
|
|
13 |
# Dec 30, 2010: - adapted for use with [LADCP_w]
|
|
14 |
|
|
15 |
#----------------------------------------------------------------------
|
|
16 |
# Volume Scattering Coefficient, following Deines (IEEE 1999)
|
|
17 |
# NOTES:
|
|
18 |
# - instrument specific! (300kHz Workhorse)
|
|
19 |
# - no sound-speed correction applied
|
|
20 |
# - R in bin center, instead of last quarter
|
|
21 |
# - transmit power assumes 33V batteries
|
|
22 |
#----------------------------------------------------------------------
|
|
23 |
|
|
24 |
# NB:
|
|
25 |
# - correction seems to work for a subset of bins (~bins 3-9 for
|
|
26 |
# 2010 P403 station 46)
|
|
27 |
# - this may imply that noise level depends on bin
|
|
28 |
# - far bins are important for seabed detection, i.e. cannot simply
|
|
29 |
# be discarded at this stage
|
|
30 |
|
|
31 |
sub Sv($$$$$)
|
|
32 |
{
|
|
33 |
my($temp,$PL,$Er,$R,$EA) = @_;
|
|
34 |
my($C) = -143; # RDI Workhorse monitor
|
|
35 |
my($Ldbm) = 10 * log10($PL);
|
|
36 |
my($PdbW) = 14.0;
|
|
37 |
my($alpha) = 0.069;
|
|
38 |
my($Kc) = 0.45;
|
|
39 |
|
|
40 |
return $C + 10*log10(($temp+273)*$R**2) - $Ldbm - $PdbW
|
|
41 |
+ 2*$alpha*$R + $Kc*($EA-$Er);
|
|
42 |
}
|
|
43 |
|
|
44 |
#----------------------------------------------------------------------
|
|
45 |
# Calculate per-bin backscatter profiles
|
|
46 |
# input: first and last valid LADCP ensemble
|
|
47 |
# output: sSv[$depth][$bin] sum of volume scattering coefficients
|
|
48 |
# nSv[$depth][$bin] number of samples in bin
|
|
49 |
#----------------------------------------------------------------------
|
|
50 |
|
|
51 |
my(@sSv,@nSv);
|
|
52 |
|
|
53 |
sub calc_backscatter_profs($$)
|
|
54 |
{
|
|
55 |
my($LADCP_start,$LADCP_end) = @_;
|
|
56 |
|
|
57 |
my(@Er) = (1e99,1e99,1e99,1e99); # echo intensity reference level
|
|
58 |
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
|
|
59 |
$Er[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0]
|
|
60 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0] < $Er[0]);
|
|
61 |
$Er[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1]
|
|
62 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1] < $Er[1]);
|
|
63 |
$Er[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2]
|
|
64 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2] < $Er[2]);
|
|
65 |
$Er[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3]
|
|
66 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3] < $Er[3]);
|
|
67 |
}
|
|
68 |
debugmsg("per-beam noise levels = @Er\n");
|
|
69 |
|
|
70 |
my($cosBeamAngle) = cos(rad($LADCP{BEAM_ANGLE}));
|
|
71 |
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
|
|
72 |
my(@bd) = calc_binDepths($ens);
|
|
73 |
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
|
|
74 |
my($depth) = int($bd[$bin]);
|
|
75 |
# next if ($depth < 0); # enable to use this code for uplookers
|
|
76 |
my($range_to_bin) = ($bd[$bin] - $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})
|
|
77 |
/ cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))
|
|
78 |
/ $cosBeamAngle;
|
|
79 |
if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP})) {
|
|
80 |
$sSv[$depth][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
|
|
81 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
82 |
$Er[0],$range_to_bin,
|
|
83 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
|
|
84 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
|
|
85 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
86 |
$Er[1],$range_to_bin,
|
|
87 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
|
|
88 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
|
|
89 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
90 |
$Er[2],$range_to_bin,
|
|
91 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
|
|
92 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
|
|
93 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
94 |
$Er[3],$range_to_bin,
|
|
95 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
|
|
96 |
} else {
|
|
97 |
$sSv[$depth][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
|
|
98 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
99 |
$Er[0],$range_to_bin,
|
|
100 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
|
|
101 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
|
|
102 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
103 |
$Er[1],$range_to_bin,
|
|
104 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
|
|
105 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
|
|
106 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
107 |
$Er[2],$range_to_bin,
|
|
108 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
|
|
109 |
Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
|
|
110 |
$LADCP{TRANSMITTED_PULSE_LENGTH},
|
|
111 |
$Er[3],$range_to_bin,
|
|
112 |
$LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
|
|
113 |
}
|
|
114 |
$nSv[$depth][$bin]++;
|
|
115 |
}
|
|
116 |
}
|
|
117 |
}
|
|
118 |
|
|
119 |
#----------------------------------------------------------------------
|
|
120 |
# determine location of seabed from backscatter profiles
|
|
121 |
# input: depth below seabed can possibly be (e.g. max CTD depth)
|
|
122 |
# output: median/mad of estimated water depth
|
|
123 |
#----------------------------------------------------------------------
|
|
124 |
|
|
125 |
sub find_backscatter_seabed($)
|
|
126 |
{
|
|
127 |
my($search_below) = int($_[0]); # grid index to begin search
|
|
128 |
my(@wdepth); # list of water_depth indices
|
|
129 |
|
|
130 |
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # find backscatter min/max below $search_below in each bin
|
|
131 |
my($minSv,$maxSv,$depthmaxSv,$lastvalid) = (1e99,-1e99,-1,-1);
|
|
132 |
for (my($depth)=$search_below; $depth<@nSv; $depth++) {
|
|
133 |
next unless ($nSv[$depth][$bin] > 0);
|
|
134 |
my($Sv) = $sSv[$depth][$bin] / $nSv[$depth][$bin];
|
|
135 |
$lastvalid = $depth;
|
|
136 |
$minSv = $Sv if ($Sv < $minSv);
|
|
137 |
$maxSv = $Sv, $depthmaxSv = $depth if ($Sv > $maxSv);
|
|
138 |
}
|
|
139 |
if ($maxSv-$minSv>10 && $depthmax!=$lastvalid) { # ignore boundary maxima & scatter
|
|
140 |
push(@wdepth,$depthmaxSv);
|
|
141 |
}
|
|
142 |
}
|
|
143 |
|
|
144 |
my($wd) = median(@wdepth);
|
|
145 |
return ($wd,mad2($wd,@wdepth));
|
|
146 |
|
|
147 |
}
|
|
148 |
|
|
149 |
1;
|