g-RIPS SENDAI 2023
MITSUBISHI-B

Final Paper

Authors:
Aya Akahori1
Emma Thomas2
Kaede Watanabe3
Cameron Watt4

Mentors:
Dr. Toshiaki Yachimura+
Ryuhei Takahashi*
Dr. Hiroshi Kameda*

1 Okayama University
2 University of Wisconsin-Madison
3 Tohoku University
4 University of Pittsburgh+ Academic Mentor, Tohoku

University* Industrial Mentor, Mitsubishi

September 26, 2023

Page 1 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Contents
1 Physical background
1.1 Relationship between Doppler spectrum and ocean wave spectra . . . . . . .
1.2 Radar Range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
3

2 Mathematical Background
2.1 Tikhonov Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Gradient Descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Dynamic Step Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Riesz Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
4
5
8

3 Methods
3.1 Defining the Cost Function for the Tikhonov Regularization Problem . . . .
3.2 Calculating the Gradient of the Cost Function . . . . . . . . . . . . . . . . .

9994 Applying the Gradient Descent Algorithm
4.1 Numerical Technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Tuning The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
12
14

5 Further Research
5.1 Improving the Gradient Descent Algorithm with Backtracking . . . . . . . .
5.2 Testing Different Wave Spectra . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Environmental Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16
16
17
17

1

Page 2 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group1

Physical background

Ocean HF radar works by emitting radio waves of different frequencies. The radar thenreceives radio waves that are reflected back by the ocean’s surface. The signal of thesereflected waves is then changed to a frequency domain equation using a Fourier transform
so that the Doppler spectrum of the reflected radio waves can be analyzed. The signals thatthe radar receives back from the ocean are determined by Bragg scattering. Bragg resonanceoccurs when radio waves that are reflected by consecutive ocean wave peaks are in phasewith one another and have maximum constructive interference. This requires the distance
between the ocean wave peaks to be half of the radio wavelength. These reflected radio waves
are said to be Bragg scattered.
First-order Bragg scattering occurs when an ocean wave moving directly towards or awayfrom the radar has a wavelength that is half the radio wavelength. We can represent bothocean waves and radio waves using wavevectors. The magnitude of each wavevector is thewavenumber of the wave, and the direction is the direction of the physical wave. If k0 is a
wavevector representing the radio wave and kd is an ocean wave moving directly towards oraway from the radar, then first order Bragg scattering requires thatkd = ±2k0 .

(1.1)

Second-order Bragg scattering involves multiple ocean waves moving in different directions. It arises from both electromagnetic and hydrodynamic effects. In the case of electrodynamics effects, the second-order Bragg scattering is a result of a radio wave bouncing off
first one ocean wave and then another in such a way that Bragg resonance occurs. Let k1
represent the first wave and k2 the second wave. In this situation, Bragg resonance occurs
when
k1 + k2 = −2k0 .
(1.2)
In the case of hydrodynamic effects, the second-order Bragg scattering occurs when the wave
reflects off the intersection of two ocean waves. As long as the two intersecting ocean waves
fulfill (1.2), Bragg resonance will occur and result in second-order Bragg scattering.
1.1

Relationship between Doppler spectrum and ocean wave spectra

Let σ(ω) be the Doppler spectrum of Bragg scattering, where ω is the Doppler frequency inradians per second (ω = 2πf ). If σ(1) is the Doppler spectrum of first-order Bragg scatteringand σ(2) is the Doppler spectrum of second-order Bragg scattering, thenσ(ω) = σ(1) + σ(2)(1.3)

There exist theoretical equations describing the relationship between the ocean waves’spectra and the Doppler spectra of first and second-order Bragg scattering. Our research
focuses on retrieving the wave spectrum, S(k) of an ocean wave, k, from the theoreticalequation relating S(k) to σ(2) given below.

2Page 3 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group6

σ(2) (ω) = 2 πk04

X

X ZZ

m1 =±1 m2 =±1|Γ(m1 k1 , m2 k2 )2 |S(m1 k1 )S(m2 k2 )
R2


p
p 
×δ ω − m1 gk1 − m2 gk2 dp dq.
(1.4)
Here, Γ is the coupling coefficient described below, k0 = |k0 | = 2π
where λ is the distanceλbetween ocean wave peaks, and delta is the Dirac delta function. The wave vectors k1 andk2 are parameterized by p and q as follows:
k0 = (k0 , 0), k1 = (p − k0 , q), k2 = (−p − k0 , −q).

(1.5)

The coupling coefficient, Γ(m1 k1 , m2 k2 ), is calculated byΓ(m1 k1 , m2 k2 ) = ΓE (k1 , k2 ) + ΓH (m1 k1 , m2 k2 )

(1.6)

where
1 (k1 k0 )(k2 k0 )/k02 − 2(k1 k2 )√2
( k1 k2 − k0 ∆))i(k1 k2 − k1 k2 )(ω 2 + ωB2 )√ΓH (m1 k1 , m2 k2 ) = −k1 + k2 −2
m1 m2 k1 k2 (ω 2 − ωB2 )with ∆ being the average normalized impedance of the seawater.
ΓE (k1 , k2 ) =

1.2

(1.7)
(1.8)

Radar Range

The area of the sea’s surface that the radar reflects off of is broken up into a grid. The
signals received back by the radar are measured for each grid section, where the sections
are separated by direction and range, as shown in Figure 1. The width of the direction andrange separations is determined by the radar’s azimuth and range resolutions, which in thisproject are approximately 15 deg and 1500 m. The total range is about 120 deg, and thetotal azimuth is about 50 km.

Figure 1: Radar Range Direction
3

Page 4 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group2

Mathematical Background
2.1

Tikhonov Regularization

Tikhonov regularization is a technique used to solve ill-posed inverse problems of the form
F (x) = y0

(2.1)

where F : X → Y is a continuous operator between Hilbert spaces X and Y . The strategy
is to consider
L(x) = ∥F (x) − y0 ∥(2.2)

and find the x which minimizes this function. For ill-conditioned F , x0 = arg min L(x),xmay not accurately approximate the true solution [10]. In such cases, we introduce theregularization term λ∥x∥ for some regularization parameter λ > 0
L(x) = ∥F (x) − y0 ∥ + λ∥x∥.

(2.3)

The regularization term controls the fitting of the solution to data. Small values of λ provide
a weak constraint to the error function and, therefore, may result in an overfitted solution.
An overfitted solution matches the data points very accurately, but may not provide a good
approximation of the true solution. Larger values of λ will overregularize the solution resulting in underfitting. For finite-dimensional linear operators A, we can choose a regularization
parameter to satisfy∥Axλ − y0 ∥ = τ ∥e∥(2.4)

xλ = ∥(A∗ A + λ2 I)−1 ∥A∗ y0

(2.5)

where

τ is a tolerance and e is the error due to noise [10]. Therefore, if we have a good estimate
for the expected amount of noise ∥e∥, then it is relatively easy to find a good regularization
parameter to use. However, for infinite dimensional non-linear problems, choosing the rightparameter can be a much more difficult task. As such, we will need to rely on experimental
methods to determine what regularization parameter to use.
2.2

Gradient Descent

Gradient descent is an iterative method used for optimizing differentiable functions. Gradient descent algorithms can be used to locate the minima of finite dimensional cost functionsof the form f : Rn → R. The first step in the gradient descent algorithm is to pick an initial
guess, x0 . Next, calculate the gradient of f at x0 . The vector −∇f (x0 ) then gives us the
direction of the steepest descent of our cost function. After calculating the gradient, we set
x1 = x0 − α∇f (x0 )
4

(2.6)

Page 5 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Groupfor a chosen step size α. For sufficiently small α, we have f (x1 ) < f (x0 ). We then iterate
this process by setting
xn+1 = xn − α∇f (xn ).

(2.7)

For small enough values of α, we will have a monotone decreasing sequence. Then if we
assume our cost function is bounded below, then the monotone converge theorem ensuresthe existence of a minimum y ∗ = min f (xn ). However, in general, the existence of lim xn isn→∞n→∞
not guaranteed.
In order for lim xn to converge to a minimizer, we need to have a starting guess that isn→∞already sufficiently close to a minimizer. This can be difficult if there is little information
about the behavior of the cost function f . One way to ensure convergence of the gradientdescent method to a global minimizer is to prove the strict convexity of the cost function. If
a function is not convex, we could attempt to locate the global minimizer of the cost function
by running the gradient descent algorithm for many different starting points. This may be
useful for determining more general properties of the cost function, such as the locationsof many different minimizers and what the basins of attraction are for different minima.
However, running the gradient descent algorithm for many different starting points may take
too long to compute.
2.3

Dynamic Step Size

When using the gradient descent method, dynamic step sizes can be useful for reducing
the number of iterations and also for finding global minima. Small step sizes will give moreaccurate trajectories toward the minimizer of the cost function. However, if the step sizes
are too small, then the algorithm will take too many iterations to converge. Larger step sizes
may not converge to a minimizer even with a good initial guess because the gradient only
provides local information about the cost function.
When using a dynamic step size, we modify equation 2.7 by
xn+1 = xn − αn ∇f (xn ).

(2.8)

Here we see that the step size an can now change as the descent algorithm runs. In orderfor the gradient descent method to converge quickly, we want to have large step sizes when∥xn − x∗ ∥ is large and small step sizes when ∥xn − x∗ ∥ is small. For the first case, the larger
step sizes cause the values of xn to quickly move through the domain of f , which allows us
to quickly find a point that is relatively close to a minimum. As we get closer to a minimum,
we then use smaller steps to hone in on a precise value for the minimizer. One problem withusing dynamic steps is determining the optimal step size. You cannot simply calculate thevalue of ∥xn − x∗ ∥ as described because that would require us to know what the minimizeris, to begin with. If you have some knowledge about what the minimum of the cost function
is, as is the case for some least squares problems, then you can adjust the step size according
to ∥f (xn ) − f (x∗ )∥. For many problems, though, we do not know the expected minimum. In
this situation, we need a different method for determining the step size.
There are a few things to consider when attempting to determine what step size to use.
First, we will note that a sufficiently large step size from xn may result in f (xn+1 ) > f (xn )5Page 6 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Figure 2: Gradient descent algorithm with insufficient decrease. Reprinted from [8].
This is a situation we definitely want to avoid since we cannot locate the minimum of f if thevalues of f (xn ) are not monotonically decreasing. Therefore, the simplest restriction that we
can introduce isf (xn+1 ) ≤ f (xn ).

(2.9)

This condition, however, is not sufficient to ensure quick convergence, as can be seen in
figure 2. We can remedy this by introducing another term that requires the next estimate to
decrease by a certain amountf (xn+1 ) ≤ f (xn ) − c1 αn ∥∇f (xn )∥(2.10)

for some constant c1 ∈ (0, 1). This is referred to as the Armijo condition. The Armijo
constant c1 is a tuning parameter that indicates how strict of a decrease we are demandingfor our function. The Armijo condition is depicted in figure 3. Values of c1 closer to 0 are lessstrict. In figure 3, this corresponds to l(α) having a smaller slope which would increase thesize of the acceptable regions. A value of c1 = 0 would correspond to l(α) being a horizontal
line and gives us the condition in equation 2.9. We exclude 0 from the possible values of the
Armijo condition because we want to improve upon the condition in equation 2.9. Valuesof c1 that are closer to 1 require a more strict decrease. A value of c1 = 1 corresponds tol(α) being the tangent line. We exclude this value because there might not be any values off (x) in the direction of the gradient that is below the tangent line, as is the case in figure 3.
For values of c1 ∈ (0, 1), we get a secant line l(α) with a slope greater than the tangent line.
There are guaranteed to be values of f (x) on this line in the direction of the gradient thatlies below the line. One issue with choosing a strict Armijo constant is that the upper boundline l(α) may be so steep that the acceptable region only exists close to the initial guess. Thisis undesirable as it results in a very short step size giving little improvement and requiring
more iterations. We, therefore, need to carefully choose our step size in order to accomplish afaster computation. As with the regularization parameter, we rely on experimental methods
to determine what parameter to choose.
In addition to the Armijo condition, we can add more conditions to restrict the acceptableregion for the step size. We would like for the next step in the iteration to be closer to a6

Page 7 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Figure 3: Example of Armijo Condition. Reprinted from [8].

Figure 4: Example of Wolfe Condition. Reprinted from [8]minimum than the current step. Without a good understanding of the cost function, it isdifficult to know if you are getting closer or further from a minimum. However, we do knowthat any minimum will have a gradient of 0 and, therefore, a directional derivative of 0
in every direction. Thus, we would like the next iteration to have a less steep directional
derivative than the current directional derivative in the direction of ∇f (xn )|⟨∇f (xn+1 ), ∇f (xn )⟩| ≤ |⟨∇f (xn ), ∇f (xn )⟩|.

(2.11)

As with the Armijo rule, this restriction is useful, but we can introduce a parameter c2 todemand further improvement. This gives us our second condition|⟨∇f (xn+1 ), ∇f (xn )⟩| ≤ c2 |⟨∇f (xn ), ∇f (xn )⟩|(2.12)

for some constant c2 ∈ (c1 , 1). This is referred to as the strong Wolfe condition. The constantc2 determines how shallow of a slope we demand. Values of c2 that are close to 1 will moreclosely approximate condition 2.11 and are therefore less restrictive. Values closer to c1 willcause the desired slope to be flatter, which provides a stricter improvement demand.
One important concern when adding these restrictions to the possible step size is whether
there exist step sizes that can actually satisfy both these conditions simultaneously. We will7

Page 8 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

prove that such a step length will always exist. However, we will first add the assumptionthat f (x) be bounded below [8].
Proof. First, we will try to find a range of values such that the Armijo condition is satisfied.
Let ϕ(α) = f (xn − α∇f (xn )) and l(α) = f (xn ) − c1 α⟨∇f (xn ), ∇f (xn )⟩. First note thatϕ(0) = l(0). Furthermore, note that−c1 ∥f (xn )∥ > −∥f (xn )∥l′ (0) > ϕ′ (α).

(2.13)(2.14)

In particular, there exists some α1 > 0 such that l(α1 ) > ϕ(α1 ). However, we note thatf (x) is bounded below, so ϕ(α) is bounded below, but l(α) is not. Therefore, there exists
some α2 > α1 such that l(α2 ) < ϕ(α2 ). Thus, by the mean value theorem, there exists an
α′ ∈ (α1 , α2 ) such that l(α′ ) = ϕ(α′ ). Then we let α′ = min{α : l(α) = ϕ(α)}. For α ∈ (0, α′ ),the Armijo condition is satisfied.
Then for the Wolfe condition, we use the mean value theorem to show that there exists
an α′′ ∈ (0, α′ ) such thatϕ(α′ ) − ϕ(0) = α′ ϕ′ (α′′ )l(α′ ) − f (xn ) = −α′ ⟨∇f (xn − α′′ ∇f (xn )), ∇f (xn )⟩−c1 α′ ⟨∇f (xn ), ∇f (xn )⟩ = −α′ ⟨∇f (xn − α′′ ∇f (xn )), ∇f (xn )⟩c1 |⟨∇f (xn ), ∇f (xn )⟩| = |⟨∇f (xn − α′′ ∇f (xn )), ∇f (xn )⟩|.

(2.15)

Furthermore, we note that c2 > c1 , soc2 |⟨∇f (xn ), ∇f (xn )⟩| > |⟨∇f (xn − α′′ ∇f (xn )), ∇f (xn )⟩|.

(2.16)

Therefore, the strong Wolfe condition is satisfied.
2.4

Riesz Representation

Until this point, we have only discussed the mathematical theory for a finite-dimensionalsetting. However, our problem requires us to apply the gradient descent method to an infinite
dimensional cost function. Thus we need to define what the gradient is for a cost function froman infinite dimensional domain. First, we will introduce the concept of the Fréchet derivative.
The Fréchet derivative generalizes the concept of a derivative to functions between normed
vector spaces. Consider a function f : V → W between normed vector spaces. We say f isFréchet differentiable if there exists a bounded linear operator Df : V → W such that∥f (x + h) − f (x) − Df (h)∥W= 0.
∥h∥V →0∥h∥Vlim(2.17)

We call Df , the Fréchet derivative of f . Now we introduce the Riesz representation theorem.
Riesz Representation Theorem 1. Given a Hilbert space, X, and T ∈ L(X, R), ∃!vT ∈ Xsuch that for all h ∈ X
T (v) = ⟨vT , h⟩(2.18)
and||vT ||X = ||T ||L(X,R) .
8(2.19)

Page 9 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Applying the Reisz representation theorem to the Fréchet derivative of f at x gives us avDf (x) such that
Df (x)[h] = ⟨vDf (x), h⟩.

(2.20)

This parallels the equation for a finite-dimensional directional derivative given by
Du f (x) = ⟨∇f (x), u⟩.

(2.21)

Therefore, for a differentiable function f between normed vector spaces, we define the gradientof f to be the Reisz representation of the derivative∇f (x) = vDf (x).

(2.22)

Using this definition allows us to apply the gradient descent algorithm to any function
defined on a normed vector space. In addition, we can use this definition to apply the Armijo
and Wolfe conditions to our problem.

3
3.1

Methods

Defining the Cost Function for the Tikhonov Regularization Problem

The problem of determining the wave spectrum S(k) from a known Doppler spectrum
σ(2) (ω) can be expressed as a Tikhonov regularization problem. First, we define the nonlinear operator A : L(R2 , R) → L(R, R) as
X X ZZ
6
4
A(S) = 2 πk0|Γ(m1 k1 , m2 k2 )2 |S(m1 k1 )S(m2 k2 )m1 =±1 m2 =±1

R2×δ(ω − m1pp
gk1 − m2 gk2 ) dp dq.

(3.1)

The right-hand side is identical to the right-hand side of equation 1.4, but A(S) is a functionof the wave spectrum rather than a function of the Doppler frequency. From this operator,
we define a cost function, L : L(R2 , R) → R2
L(S) = ||A(S) − σ(2) (ω)||2L2 (R) + λ||S||2L2 (R2 )(3.2)

where λ is our regularization parameter. Our goal is to minimize this cost function using a
gradient descent algorithm to determine the wave spectrum.
3.2

Calculating the Gradient of the Cost Function

We need to determine an expression for the gradient of the cost function, but first, we
will show that A(S) is Fréchet differentiable.

9Page 10 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Proof.
6

A(S + h) =2 πk04

X

X ZZ

m1 =±1 m2 =±1|Γ(m1 k1 , m2 k2 )2 |[S(m1 k1 ) + h(m1 k1 )]R2pp
× [S(m2 k2 ) + h(m2 k2 )]δ(ω − m1 gk1 − m2 gk2 ) dp dq(3.3)
X X ZZ=26 πk04|Γ(m1 k1 , m2 k2 )2 | S(m1 k1 )S(m2 k2 ) + S(m1 k1 )h(m2 k2 )m1 =±1 m2 =±1

R2pp
× S(m2 k2 )h(m1 k1 ) + h(m1 k1 )h(m2 k2 )]δ(ω − m1 gk1 − m2 gk2 ) dp dq (3.4)

X X ZZ
6
4
2=A(S) + 2 πk0|Γ(m1 k1 , m2 k2 ) | S(m1 k1 )h(m2 k2 )m1 =±1 m2 =±1

R2


p
p
+ S(m2 k2 )h(m1 k1 ) δ(ω − m1 gk1 − m2 gk2 ) dp dq + o(h).

(3.5)

Therefore, the Fréchet derivative of A is given by
X X ZZ′
6
4
2
A (S)[h] =2 πk0|Γ(m1 k1 , m2 k2 ) | S(m1 k1 )h(m2 k2 )m1 =±1 m2 =±1

R2


p
p
+ S(m2 k2 )h(m1 k1 ) δ(ω − m1 gk1 − m2 gk2 ) dp dq.

(3.6)

We can then use this to express the Fréchet derivative of our cost function
L′ (S)[h] = 2⟨A(S) − σ2 (ω), A′ (S)[h]⟩L2 (R) + 2λ⟨S, h⟩L2 (R2 ) .

(3.7)

Now we need to determine the Riesz representation of the derivative of the cost function.
The Riesz representation of L′ (S) is a function vS ∈ L2 (R2 ) such that
ZZ′
L (S)[h] =vS (p, q)h(p, q)dpdq.
(3.8)
R2

One way for us to calculate vS is to rewrite the integral for L′ (S)[h] in a way that allows
us to factor out the h. In order to do this, we need to make sure that each instance of h inequation 3.6 has the same argument. This can be accomplished by splitting up A′ (S) into 8terms, each of which takes the form
X X ZZ
p
p
6
4
2
2 πk0
|Γ(m1 k1 , m2 k2 ) | S(mi ki )h(mj kj ) δ(ω − m1 gk1 − m2 gk2 ) dp dq,
m1 =±1 m2 =±1

R2

(3.9)
i, j ∈ {1, 2}, j ̸= i. In order to factor out h from each term, we want to perform substitutions
so that the arguments of h in each term are identical. For each term, we make the substitution(mj kj ) = (x, y). When j=1, this substitution gives us x = m1 (−p − k0 ) and y = −m1 q, thensolve for p and q to get p = mx1 + k0 and q = − my2 . When j=2, this substitution gives us
10Page 11 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Groupx = −m2 (p − k0 ) and y = −m1 q, then solve for p and q to get p = − mx2 − k0 and q = − my2 .
We then substitute these values as the p and q parameters of kj . As an example, we will lookat the term where i = 1, j = 2, and mi = mj = 1. In this case, we make the substitutionx = −p − k0 and y = −q, then solve for p and q to get p = −x − k0 and q = −y. Substituting
these values in as the parameters for k1 gives us: 2X X ZZ642 πk0
Γ (−x − 2k0 , −y), (x, y)S((−x − 2k0 , −y)h(x, y)(3.10)
m1 =±1 m2 =±1

R2



√ 2
√
2
2 14
2 14
× δ ω − g[(x + 2k0 ) + y ] − g[x + y ] dxdy.

(3.11)The remaining terms can be calculated with similar substitutions. Grouping results,
factoring out the h(x, y), and adding the normalization term gives us2m1
vS (x, y) =4(A(S) − σ(2) (ω))Γ − (x + 2m2 k0 ), y), (x, y)m2−∞
m1 =±1 m2 =±2


# 41
"

2

1
1
m
√
√
1
2
2
2
x + 2k0 + y− m2 g[x + y ] 4  S − (x + 2m2 k0 , y) dω×δ ω − m1 g
m2
m2
Z ∞

X

X

+ 2λS.

(3.12)

We can now define∇L(S) = vS (x, y).

(3.13)

Thus, we can use this to define our gradient descent recursion as
Sn+1 = Sn − αn vSn .

(3.14)The expression for the gradient of L also allows us to rewrite the Armijo condition as
L(Sn+1 ) ≤ L(Sn ) − c1 αn ∥vSn ∥L2 (R2 ) .

(3.15)

Similarly, we can rewrite the Wolfe condition as|⟨vSn+1 , vSn ⟩L2 (R2 | ≤ c2 |⟨vS , vS ⟩|L2 (R2 ) .

4(3.16)

Applying the Gradient Descent Algorithm
4.1

Numerical Technique

We suppose that we are given a second-order Doppler spectrum σ(2) (ω) with the aim offinding the wave spectra S(k) that produced it. First, we interpolate the data of σ(2) (ω)
to get a function to R. The radar provides data for ω ∈ [−4.6116π, 4.6116π]. We thenchoose an initial wave spectrum S0 to begin the gradient descent algorithm. A partition of
the pq-plane is chosen to numerically calculate the integral of A(S0 ). The current version of
11

Page 12 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

our code breaks the relevant segment of the pq-plane into a 30 × 30 grid. We note that the
double integral in equation 3.1 includes a delta function. This delta function defines a curveof integration, but the argument of the delta function is too complicated to find an explicit
parameterization. Therefore, we approximate the delta function as(1
if |x| < ϵ2ϵ
δ(x) ≈.
(4.1)
0 otherwise
Once an approximation for A(S0 ) is calculated, the Reisz representation is numerically calculated using the approximation for δ(x) again. Lastly, the original function S0 is adjusted
according to equation 3.14.
4.2

Tuning The Algorithm

There are a couple of challenges that we need to address when using the Gradient descentalgorithm. The first is that the behavior of L(S) is difficult to predict. Therefore, we areunable to determine if the solution we get from the gradient descent algorithm is the actual
solution or if there is another wave spectrum that produces the same minimum. One way
to fix this would be to figure out how to determine a good starting guess for the gradientdescent algorithm. However, this requires more information about the characteristics that
could affect the wave spectra.
Additionally, we know that L(S) is bounded below 0, but we have no information about
how small the minimum value of L(S) is expected to be. Therefore, it is unknown if we areclose to a minimum or not at any step in the iteration. This makes it difficult to determine
when to stop the algorithm. One way to determine a stopping point would be to check forthe rate of improvement and stop once the improvement slows to a chosen amount.
Our main method for determining what parameters to use for the gradient descent algorithm will be to train the algorithm on Doppler spectra created from known wave spectra.
When constructing the wave Doppler spectra, we want to use functions that closely model
real wave spectra. We use the model from Kataoka [6]. First, we describe the wave spectra
as a function of frequency and direction:
S(k) = Sf (f )G(θ|f )(4.2)

where Sf is the frequency spectra and G is the wave spectra distribution. The frequency
spectra are given bySf (f ) = 0.257Hs2 Ts (Ts f )−5 exp −1.03(Ts f )−4(4.3)
where Hs is the height of the significant wave and Ts is the period of the significant wave.
The wave spectra distribution is given byθ − θ02s̃(f )
(4.4)
G(θ|f ) = G0 cos
2where G0 is a normalization constant and s̃ is the wave index. The wave index is given by 5smax f
if f ≤ fpf p −2.5(4.5)
s̃(f ) =f
s
if f > fp
max fp12

Page 13 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

where smax is a parameter and fp is the peak wave frequency. Since G(θ|f ) is a distribution,
we require that
Z πG(θ|f )dθ = 1
(4.6)
−π

Therefore,
Z π
G0 =

2s̃(f )cos−π

θ − θ0
2

−1


dθ

.

(4.7)

An example of a wave spectrum constructed from this model is shown in figure 5. The figure
indicates a dominant wave traveling with a frequency of f = T1s = 0.33 in the directionof θ0 = 0.5. The mass of the wave spectra is concentrated near T1s , θ0 and continuously
decreases to 0 at further points.

Figure 5: Model wave spectrum for Hs = 1, Ts = 3, smax = 10, fp = 0.317, θ0 = 0.5 in thefrequency-direction plane.
The wave spectrum is constructed in the f θ-plane with f ∈ [10−2 , 1] and angle values ofθ ∈ [−π, π]. The spectra is then converted to the pq-plane with the relations pθ = arctan
(4.8)
q
1√ 2
f=
g(p + q 2 )1/4 .
(4.9)
2πThe model wave spectrum is then used to construct a second-order Doppler Wave spectrum as described in section 4.1. The second-order Doppler spectrum resulting from thewave spectrum in figure 5 is depicted in figure 7. For this computation, ϵ = 0.01 was used
for the√ delta function approximation. One important thing to note is that for a frequency
ω = 2gk0 , there is only first-order Bragg scattering, so we expect an absolute power of 0 atthese values.
13

Page 14 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Figure 6: Model wave spectrum for Hs =1, Ts = 3, smax = 10, fp = 0.317, θ0 = 0.5
in the pq-plane.

4.3

Figure 7: Second order Doppler spectrum
constructed from the model wave spectrum in figure 5.

Numerical simulation

We solved the direct problem by using the wave spectrum with Hs = 1, Ts = 3, smax = 10,
fp = 0.317, θ0 = 0, and constructing σ2 . These are the contour lines of the sample wave
spectrum (figure 8) and constructed σ2 (figure 9).

Figure 8: The model of wave spectrumwith parameters Hs = 1, Ts = 3, smax =10, fp = 0.317, θ0 = 0.

Figure 9: Second-order Doppler spectrum
constructed from the model wave spectrum in the figure.

By using constructed σ2 , we tried to reconstruct the initial value. The algorithm is as
below.
1. Set initial S0 , n = 0.
2. Using the Riesz representation formula, find the steepest gradient, ∇L(S).
3. From the formula, Sn+1 = Sn − α∇L(S), generate next S value.
14

Page 15 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

4. Iterate the process from 2 to 3 until convergence.
By setting the initial condition to get σ2 as S0 , the result of the numerical simulation is in
Figure 10. The parameters are step size = 1, iteration number = 30, regularization term =10−4 , and calculation time is about 54 sec. This is the almost same as the initial condition.

Figure 10: The same figure 8

Figure 11: Reconstruction.

Next, we add perturbation to the initial data. Let ϵ be the amplitude of perturbation.
When ϵ = 0.001, the initial condition is figure 11.

Figure 13:0.001.

Figure 12: ϵ = 0.001

Reconstruction when ϵ =We tried to reconstruct the initial condition (figure 8) from this and changed the iteration
number = 200. The calculation time is about 267 sec, and we change the value less than
0.002 to 0 (figure 12).
When ϵ = 0.003, the initial condition is figure 14, and we reconstructed from this withthe same parameters of ϵ = 0.003. The calculation time is 263 sec, and the result is figure 15
with changing the value from less than 0 to 0.
15

Page 16 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group

Figure 15:0.003.

Figure 14: ϵ = 0.003

Reconstruction when ϵ =From these results, we can almost reconstruct the initial condition with regularization
terms and changing minimum values. However, it is difficult to determine what parameters
to use to reconstruct the initial conditions. In the future, we need to find out a good way to
set the parameters and a stronger system to reconstruct the wave spectrum without changing
minimum values.

5
5.1

Further Research

Improving the Gradient Descent Algorithm with Backtracking

The main improvement that we would like to add to the gradient descent algorithm is
the Armijo and Wolfe conditions mentioned in section 2.3. These conditions would hopefully increase the rate of convergence of the algorithm. In addition to improving convergence
speed, the Armijo and Wolfe conditions would also improve the ability of the gradient descent algorithm to avoid local minima and locate a global minimum. In order to implement
these conditions, we would use a backtracking line search algorithm. The step size for a
backtracking line search at the nth step in the gradient descent is defined by
αn = M β mn

(5.1)

where M > 0 and β ∈ (0, 1) are constants, and mn is the first non-negative integer suchthat both the Armijo and Wolfe conditions are satisfied. The backtracking algorithm iterates
through test step sizes
α = M βm(5.2)

over m until the conditions are met, and then the nth step size for the gradient descent isdefined accordingly. One issue with this algorithm is that for poorly chosen M and β, theArmijo and Wolfe conditions may never be satisfied for the tested step sizes. In section 2.3,
we proved the existence of a step size that satisfies both conditions. However, this does not
16

Page 17 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Groupguarantee that the backtracking line search will find such a value. The backtracking algorithm
tests for the conditions at discrete values defined by equation 5.2. When moving from one
test step to the next, some information about the cost function is ignored. Referring back tofigure 4, we can see that the acceptable regions are positive intervals on the real axis that do
not contain 0 as an endpoint. This is an issue because it is possible for a backtracking searchto skip over these acceptable regions, and if the test step ends up between 0 and the infimumof the acceptable regions, then the algorithm will be stuck in an infinite loop. In addition toskipping over acceptable regions, poorly chosen backtracking parameters may also result inthe test steps skipping over a global minimum, which could lead to an inaccurate result. In
order to gain more information about the line search, we need to use a larger backtracking
parameter β. This will give a better chance of finding a value in an acceptable region, but it
may result in unnecessary calculations if the test step is far from an acceptable region.
The backtracking algorithm requires us to choose two parameters: the max possible step
M and the backtracking parameter β. The magnitude ∥S∥L2 (R2 ) depends on the significant
wave height and frequency. Therefore, if reasonable bounds can be obtained for either ofthese values, the Kataoka model [6] could be integrated for these extreme values to obtain abound S for ∥S∥L2 (R2 ) . Then for two arbitrary wave spectra S1 and S2 ,∥S1 − S2 ∥L2 (R2 ) ≤ ∥S1 ∥L2 (R2 ) + ∥S2 ∥L2 (R2 ) ≤ 2S.

(5.3)

If bounds for the wave height and frequency cannot be obtained, the max step size could bedetermined experimentally. The backtracking parameter is much more difficult to determine
because of the unpredictable behavior of the cost function 3.2. One strategy for deciding
the backtracking constant would be to start with large values of β, which give the mostinformation about the cost function. Since a large value of β would result in the backtracking
algorithm scanning the cost function for detail, these values would be expected to locate thecorrect minima. Then the value of β could be decreased to a βl where the gradient descent nolonger converges to the desired minima. This would provide a lower bound for β. Running
this procedure for many different gradient descent trials will produce a set B of βl . Then wecan choose β to be greater than or equal to sup(B).
5.2

Testing Different Wave Spectra

Another way to improve the gradient descent is to test the gradient descent algorithm
by attempting to recover more complicated wave spectra. The main wave spectra we would
consider would be a summation of wave spectra from the Kataoka model [6]. Using a sumof two or more wave spectra with different parameters would allow us to represent surface
waves composed of multiple wave trains. Testing the gradient descent algorithm on more
complicated wave spectra could allow us to fine-tune the parameters to work in a broader
range of scenarios. This may make the gradient descent code more accurate when testing on
a Doppler spectrum of unknown wave spectra.
5.3

Environmental Effects

Lastly, to optimize our algorithm, we would like to use environmental factors to try to
determine a good initial guess for the gradient descent. We hope to use the Kataoka model [6]17Page 18 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Groupfor starting spectra. However, this requires us to define parameters for the significant wave
height Hs , significant wave period, Ts , and wave direction θ0 . There are many environmental
factors that influence surface waves, including ocean currents and ocean depth.
First, we will discuss the effects of currents. The main effect of ocean currents is that wave
trains tend to align with ocean currents [13]. Therefore, if we have information about thedominant currents in the radar range, then we can use this to define a direction for the initialwave spectrum. It is important to note, however, that the current direction and dominant
wave train direction do not necessarily coincide. Therefore, it might be useful to look for analternative method to determine the dominant wave train direction. When water currents
oppose the wave currents, the waves are slowed. Since the wave energy is determined by its
speed and height, conservation of energy causes the wave height to increase [13]. Therefore,
if we have knowledge about the average wave height in the radar range, then data about theocean currents could be used to determine what wave height to use for our initial guess.

Figure 16: Bathymetric maps of Muroto coast and Ocean Current map
The ocean depth also affects the wavefronts. In particular, shallow areas tend to cause
waves to move slower. Again, this leads to an increase in wave height. Additionally, if the
wavefronts are not aligned with the ocean depth contours, then the slowing in shallower
waters causes the wave direction to change due to refraction. This could also be a usefulfactor in determining the wave direction for an initial guess. However, it would still require
some initial knowledge of the wave direction at some ocean depth to be useful. Figure 16
shows the depth of the ocean around the Muroto coast. The yellow line represents a depth of
200m. The red line represents a depth of 1000m. At this contour line, there is a sudden dropin ocean depth. Waves move faster in deeper water, but as they enter shallow water, theyslow, and the conservation of energy increases wave height. We hope that these effects canbe explored in more detail to determine how this affects the parameters in our initial guess.

References[1] Peter Mathé Abhishake Rastogi Gilles Blanchard. Convergence Analysis of Tikhonov
Regularization for Non-Linear Statistical Inverse Problems. 2020.

18

Page 19 of 19

g-RIPS Sendai 2023, MITSUBISHI-B Group[2] Japan Meteorological Agency. Daily 50m Current. 2023. url: https://www.data.
jma.go.jp/kaiyou/data/db/kaikyo/daily/current_HQ.html?areano=8&stryy=
2023&strmm=08&strdd=03.
[3] D. E. Barrick. “Extraction of Wave Parameters from Measured HF Radar Sea-Echo
Doppler Spectra”. In: Radio Sci, 415–424 (1977).
[4] Christian Cason. Regularization of Inverse Problems. Feb. 2021.
[5] J.M.J. Journée and W.W. Massie. Offshore Hydromechanics. Delft University of Technology, 2001.
[6] Tomoya Kataoka. “A Bayesian Inversion for a Directional Spectrum of Ocean Waves
in Shallow Water Using HF Radar”. In: (2016).
[7] Thuy T. Le and Loc H. Nguyen. “The gradient descent method for the convexification tosolve boundary value problems of quasi-linear PDEs and a coefficient inverse problem”.
In: Journal of Scientific Computing (2022).
[8] Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, 2006.
[9] J. Paduan and H. Graber. “Introduction to High-Frequency Radar: Reality and Myth”.
In: Oceanography 10 (1997).
[10] Y. Park et al. “Parameter Determination for Tikhonov Regularization Problems in
General Form”. In: Journal of Computational and Applied Mathematics (2018).
[11] Aditya Tiwari. “Understanding Overfitting and Underfitting in Machine Learning”. In:(2019).
[12] S. Twomey. Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements. 1977.
[13] Trevor Whittaker and Jonathan Hancock. Environmental Loads, Motions, and Mooring
Systems. 2020.

19