# Geotechnical Parrot Tales (GPT): Harnessing Large Language Models in
geotechnical engineering

Krishna Kumar

##  1 Introduction

The rise of large language models (LLMs) has generated widespread interest due
to their ability to answer questions, generate text/code, and summarize. One
of the most popular LLM is OpenAI’s ChatGPT, released in November 2022.
ChatGPT is a probabilistic LLM built on OpenAI’s GPT-3.5 (Generative
Pretrained Transformers) with a chat user interface to the GPT model. Despite
their impressive capabilities, ChatGPT can produce plausible-sounding but
false outputs, blurring the line between facts and hallucinations [1].

Prompt engineering is creating effective inputs to help LLMs generate the
desired output. Prompt engineering is crucial to mitigate these risks and
harness the full potential of ChatGPT in geotechnical engineering. Prompt
engineering refers to carefully designing input prompts to elicit accurate and
valuable responses from LLMs [3]. By developing a deep understanding of
ChatGPT’s strengths and limitations, engineers can create prompts that
effectively guide the model toward generating valuable insights while
minimizing the risks of falsehoods and hallucinations.

In this article, we will explore ChatGPT and its applications in geotechnical
engineering. We will discuss the challenges and pitfalls associated with these
models and highlight the importance of prompt engineering in navigating these
challenges to ensure reliable and accurate outcomes.

##  2 Demystifying ChatGPT: The technical bits

The groundbreaking paper “Attention is All You Need” introduced Transformers
[5], the ‘T‘ in the GPT models. Unlike previous natural language models,
transformers could look at the whole context. Imagine we are reading a book.
As you progress through the sentences, you naturally anticipate the next word
based on your understanding of language and the context provided by the story.
Think of the Transformer as a dedicated reader who has read thousands of books
and learned to recognize patterns in how words and sentences are structured.
When encountering a new sentence or paragraph, it uses this knowledge to make
educated guesses about what word should come next. For example, if the
Transformer reads a sentence like ”During site investigation, the engineer
discovered that the soil has a high degree of…”, it would likely predict the
next word as something like ”saturation” or ”compaction,” based on its
understanding of common phrases and concepts in the geotechnical field. The
Transformer model achieves this by utilizing layers of attention mechanisms.
These attention mechanisms help the model focus on the most relevant words or
phrases in the given context, just like how we might focus on specific details
when trying to predict while reading a book.

Figure 1a illustrates how ChatGPT might complete the sentence, ”The
geotechnical engineer discovered that a site with dry collapsible soil in arid
Arizona  .” GPT chooses each subsequent word based on probabilities derived
from its extensive training corpus and the text generated thus far. This
example selects the highest-ranking word at every step, resulting in the
response, ”… was not suitable for the proposed construction.” While this
method works well for short responses, it can take longer responses to appear
robotic, stiff, and unimaginative. To address this, GPT provides various
controls for adjusting word selection. One such control is the temperature
setting. With a temperature of 0, GPT always chooses the word with the highest
probability (fig. 1a). However, by increasing the temperature to 1.0 (ChatGPT
uses 0.7), GPT occasionally selects lower-ranked words, producing text that
often feels more natural and inventive. Figure 1b demonstrates this: a non-
zero temperature setting generates an impossible and inaccurate response ”was
highly saturated with water” from the same input prompt. A zero temperature
means the GPT model will be more consistent (the output will not change
regardless of how often you run the same query), but it may still be wrong.
Most scientific applications would require a temperature of zero to guarantee
consistent results.

ChatGPT is susceptible to hallucinations, meaning that it can occasionally
generate responses that are not based on factual information or are unrelated
to the context. These hallucinations can happen even when the temperature
parameter is zero. A model’s capability refers to its objective or goal, while
alignment focuses on the desired behavior of the model compared to its actual
training. ChatGPT’s goal is to produce human-like text based on input, while
its alignment aims to ensure that generated responses meet the expectations
and values of its users. Models such as GPT-3 are capable but misaligned, as
their training does not always capture higher-level meaning. This creates a
noticeable gap between how these models are trained and how we want to use
them. ChatGPT employs Reinforcement Learning from Human Feedback to better
align with human values and expectations. GPT-3 has 175 billion parameters,
while GPT-4 [2] is speculated to have 1 trillion parameters to encode all
training information. ChatGPT generates plausible-sounding responses using
these parameters alone, which leads to hallucinations, limited
interpretability, and potentially biased output. Involving humans in the loop
can introduce new biases. They may lack the required expertise for specific
tasks (ChatGPT trainers may not include experts in Geotechnical Engineering),
resulting in less accurate feedback.

ChatGPT can be described as a stochastic parrot because it learns from vast
amounts of text data and generates human-like text based on the patterns it
identifies. However, due to its training strategies, it may produce
inconsistent, biased output or lack deeper understanding, much like a parrot
that can mimic human speech but does not fully grasp the meaning behind the
words. While ChatGPT is competent, addressing hallucination and alignment
issues remains crucial to ensure its outputs align with human values and
expectations.

Figure 1: GPT text generation using transformers.

##  3 Applications of ChatGPT to Geotechnical Engineering

Even with the current advancements in GPT-4, the models will hallucinate,
i.e., lie or confidently make things up. Although ChatGPT is widely used to
showcase its generative power, like writing emails and reports, we have not
explored the innovative use cases of LLMs in engineering. The best use of
ChatGPT is not as an intelligent answering machine but as a reasoning engine
or a “calculator of words.” We will explore advanced ChatGPT use cases that
can accelerate geotechnical engineering design and research.

###  3.1 Prompt Engineering

ChatGPT has a huge potential to accelerate the adoption of data-driven methods
and computational modeling in Civil and Geotechnical Engineering. When we
provide ChatGPT with a specific request called a “prompt,” it uses its vast
knowledge to create a suitable response. However, crafting an effective prompt
is crucial, as it guides the GPT in understanding what we want – “prompt
engineering.” A well-designed prompt ensures GPT comprehends the task and
delivers accurate, efficient, and functional code. Context is vital in prompt
engineering for code generation because it helps the ChatGPT to better
understand the specific requirements and nuances of the task at hand.
Providing sufficient context allows the model to generate accurate and
relevant code for the problem we are trying to solve. Without adequate
context, ChatGPT might generate code that appears functional but lacks the
necessary details to address the issue effectively.

Consider the prompt: ‘Write a Python code for USCS soil classification.’ In
this case, ChatGPT hallucinates and responds with an incorrect code without
including Atterberg’s limit, yet indicating it wrote the code based on ASTM
D2487-17 (see fig. 2). The prompt assumes that ChatGPT possesses all the
necessary knowledge to accomplish the task. However, users should primarily
use ChatGPT as a translator, providing the necessary context and knowledge to
accomplish the task. We can provide a context by including “Context: ” in the
chat prompt and adding the relevant information. For the USCS soil
classification, we will treat ChatGPT as a translator from English to Python.

Prompt: Write a Python code for USCS soil classification using the following
conditions and using the A-line equation as 0.73 * (LL - 20).
Context:
Inputs: % passing sieve 4 and 200, D10, D30, D60, Liquid limit (LL), and
Plastic Limit (PL).
Compute auxiliary variables: Plasticity Index $(PI)=LL-PL$, $Cu=D60/D10$, and
$Cc=D30^{2}/(D10*D60)$
$<$Showing only Gravel classification for brevity$>$

  1. 1.

If the percentage passing through the No. 200 sieve (Sieve200) is less than or
equal to 50, the soil is coarse-grained:

    1. (a)

If the percentage of gravel is greater than the percentage of sand, classify
the soil as follows:

      1. i.

If the percentage passing through the No. 200 sieve is less than 5 and the
coefficient of uniformity (Cu) is greater than 4, and the coefficient of
curvature (Cc) is strictly between 1 and 3, classify the soil as ”GW”.

      2. ii.

If the above conditions are not met, classify the soil as ”GP”.

    2. (b)

If the percentage passing through the No. 200 sieve is greater than 12,
classify the soil as follows:

      1. i.

If the plasticity index (PI) is below the A-line, classify the soil as ”GM”.

      2. ii.

If the plasticity index is above the A-line, classify the soil as ”GC”.

      3. iii.

If the plasticity index is on the A-line, classify the soil as ”GM-GC”
(borderline case).

    3. (c)

For borderline cases $(5\%\leq\textit{}Sieve200\leq 12\%)$, consider the
following classifications:

      1. i.

If Cu $>4$ and Cc is strictly between 1 and 3, check the plasticity index: If
PI is below the A-line, classify the soil as ”GW-GM”. If PI is above the
A-line, classify the soil as ”GW-GC”. If PI is on the A-line, classify the
soil as ”GW-GM-GC” (borderline case).

      2. ii.

If the conditions above are unmet, check the plasticity index: If PI is below
the A-line, classify the soil as ”GP-GM”. If PI is above the A-line, classify
the soil as ”GP-GC”. If PI is on the A-line, classify the soil as ”GP-GM-GC”
(borderline case).

Figure 2: ChatGPT response to write a Python code for soil classification.

Similarly, ChatGPT can summarize articles and reports, but it is important to
provide a context. For example, a summarization prompt without a context, such
as, Write a summary of “Crocker, Kumar, and Cox. (2023) Using explainability
to design physics-aware CNNs for solving subsurface inverse problems.” only
gives a hallucinated summary. ChatGPT does not know the paper, as its
knowledge is limited to articles published before September 2021. Even if we
ask ChatGPT to summarize a paper published before September 2021, without
context, ChatGPT will hallucinate. Accurate summarizing is only possible when
the whole paper is provided as context. To summarize a soil report, we could
ask ChatGPT specific information that needs summarization “Summarize the
strength parameters of the clay layer in the soil report: $<$insert
report$>$”.

###  3.2 A context-specific search engine

This section offers an approach to proving the correct context to LLMs. In the
following sections, we will use OpenAI’s GPT-3.5-turbo Application Programming
Interface (API) rather than the chat interface for context-specific Q&As. LLMs
are reasoning engines that analyze textual information to extract answers
based on context. Unless we provide an appropriate context, GPT will suffer
from hallucinations. For example, if we ask GPT, “What is the XML tag to store
plastic limit in DIGGS?” we get the following incorrect response: “<TestType
>Plastic Limit</TestType>.” However, the XML tag test type is a generic tag
GPT invented to store Plastic Limit, although GPT has been trained on this
data and is aware of DIGGS. For GPT to generate answers, we need to provide
appropriate context. In this case, we need to identify the appropriate
sections of the DIGGS schema and provide it as a context. The number of tokens
(words) we can provide as context is limited (32k), and the compute cost is
related to the number of tokens. Hence, it is important to find the
appropriate context rather than providing the entire knowledge base, which may
not be possible in many cases. How can we build a context-specific search
engine that uses local knowledge bases to find the correct context for GPT to
answer the question?

Figure 3: Context-specific semantic search engine using vector database and
latent space embedding.

Semantic search aims to understand the intent and contextual meaning of search
queries, providing more accurate and relevant results. Figure 3 shows an
overview of the semantic search engine architecture. Using OpenAI’s text
embedding model, we convert the textual knowledge base (e.g., DIGGS schema)
into context vectors (see fig. 3). The context vectors of the knowledge base
are stored in a fast-accessible vector database like FAISS. When a user enters
a search term, we convert the search term to a vector using the same embedding
model. We perform a cosine similarity on the vector database to identify the
appropriate context. Cosine similarity involves computing the dot product
between the search vector and every knowledge vector in the database. OpenAI
uses an absolute positional embedding, which means a cosine similarity of 1.0
means the context is an exact match to the search query. We then rank the
contexts based on the dot product values (search term vs. knowledge vector).
Once we have identified the top five contexts matching the search query, we
feed those as a context to GPT API (we use a GPT generation model like `text-
davinci-003`) to create a customized answer. This context-specific GPT will
guarantee deterministic answers and prevent hallucinations. Using this
context-specific GPT, we get the correct answer “The DIGGS XML tag for plastic
limit is <diggs_geo:waterContent>within the parent tag
<diggs_geo:plasticLimitTrial>” We can reuse the search engine to do more
complex tasks, such as asking it to generate a DIGGS XML file (see fig. 4). We
can also finetune the semantic search engine using supervised learning trained
on a list of questions and answers about DIGGS.

Prompt: I have a list of three plastic limit trials: [11.9 %, 11.7 %, and
11.4%]. Autoincrement trial number starting with one. Add a ‘gml:id‘ string
for each trial. Create the XML file.

Response: Each plastic limit trial is represented by a
‘<diggs_geo:PlasticLimitTrial>’ element inside a
‘<diggs_geo:plasticLimitTrial>’ element. The ‘gml:id’ attribute is set to
“tr1”, “tr2”, and “tr3” for the three trials, respectively. The ‘trialNo’
element is auto-incremented starting from 1 for each trial, and the
‘waterContent’ element contains the percentage values for the trials. The
‘isManual’ element is set to true for all trials, but you can change it to
false if needed.

Figure 4: Using GPT to generate DIGGS Schema.

###  3.3 Reasoning with a train of thought

Can we make GPT do geotechnical analyses? [6] showed that instead of
finetuning a separate language model for each new task, it is possible to just
”prompt” the model with a handful of input-output examples to demonstrate the
task. Let us now consider if we can provide examples to GPT to classify fine-
grained soils based on liquid and plastic limits.

Prompt: Using soil classification of fine-grained soils as CL, CH, ML, or MH,
classify the following soil with a Liquid Limit of 30, Plastic Limit of 10.

Response: In this case, PI = LL - PL = 30 - 10 = 20. Since the PI of the given
soil is 20, it falls into the category of High Plasticity Soil. Therefore, the
soil can be classified as either CH (High Plasticity Clay) or MH (High
Plasticity Silt) based on visual and textural characteristics.
GPT incorrectly classifies the soil as CH or MH when the liquid limit is less
than 50%. How can we build a reasoning engine that understands the underlying
logic? We explore how reasoning abilities emerge naturally in large language
models via a simple method called chain-of-thought prompting. Imagine solving
a geotechnical engineering problem, like calculating the number of trucks
needed to construct an embankment. We typically break down the problem into
smaller steps and solve each before arriving at the final answer: ”If the
embankment requires 500 cubic meters of soil and each truck carries 25 cubic
meters, we need 20 trucks… then, considering a 10% compaction loss, we need an
additional two trucks… so the total number of trucks required is 22.” We show
how to provide language models to classify soils with a similar thought
process—a coherent series of intermediate reasoning steps leading to the
answer.

Prompt: Context: Fine-grained soil classes (“CL,” “CH,” “ML, “MH”)
Q: What is the classification of fine-grained soil with a LL = 60 and a PL =
50?
A: LL boundary of 50 classifies if it is “L” low plastic or “H” high plastic.
LL = 60 >50, so “H”. PI = LL - PL = 60 - 50 = 10. Check if PI is above A-line:
If PI >0.73 * (LL - 20) = 10 >0.73 * (60 - 20), if it is above the A-line, it
is “C” for clay, else “M” for silt. It is not so “M” silt. Soil is “MH.”
Q: What is the classification of fine-grained soil with a LL = 30 and PL = 10?

Response A:LL boundary of 50 classifies if it is “L” low plastic or “H” high
plastic. LL = 30 <50, so “L”. PI = liquid limit - plastic limit = 30 - 10 =
20. Check if PI is above A-line: If PI >0.73 * (LL - 20) = 20 >0.73 * (30 -
20). It is above the A-line, so “C” clay. So the soil is “CL.”

GPT now correctly classifies the soils as low-plastic clay, following a
similar reasoning process. The train-of-thought prompting provides a reasoning
ability to GPT and facilitates knowledge discovery.

Let us consider another example, a task of assessing the safety of a gravity
retaining wall. The wall has specific dimensions and is acting upon particular
soil properties. We’re to ensure the wall is safe against sliding,
overturning, and bearing, applying a safety factor of 1.25. In train-of-
thought reasoning, we will break down each step into separate prompts and
provide sufficient context for GPT to solve the problem. For e.g., “the first
step involves evaluating the safety against sliding. This means we have to
compute the ratio of the driving force (active soil pressure) to the resisting
force (frictional resistance at the base of the wall). If the calculated
factor of safety is above 1.25, it implies that the wall is safe against
sliding. In this scenario, we’re assuming passive pressure isn’t considered,
as it may be necessary to remove the earth in front of the wall for placing
services, leaving the active pressure to be resisted only by the friction at
the base of the wall.” Similarly, we can provide detailed context for GPT to
perform calculations to check the safety against overturning by asking it to
calculate the total moment about the base center to find the eccentricity ’e’
and state that If the eccentricity is within the middle third of the base
width, it would imply that the wall is safe against overturning. Finally,
prompt GPT to check for safety in bearing pressure by using Meyerhof’s reduced
base width approach. If the pressure exerted by the wall is less than the
bearing capacity of the soil (considering the safety factor of 1.25), the wall
would be considered safe against bearing failure. In order to guide the model
in understanding these steps, prompt engineering must be carefully designed,
with each step detailed out clearly. This would allow the GPT model to follow
a logical sequence and provide a solution that aligns with the geotechnical
reasoning required in these evaluations.

Although, train-of-thought process offers a reasoning engine with GPT, it is
still limited by the formulas used by GPT, which may be hallucinated, leading
to incorrect solutions. The next section explores how to build a robust
engineering framework with GPT.

###  3.4 GPT-enabled engineering workflows

Let us now consider a complex geotechnical engineering workflow of calculating
the maximum load supported by a foundation. Typically, most engineering
workflows need several calculation stages. Each stage feeds into the next one,
leading to the final result. Language models, like GPT, have limitations. They
cannot solve these complex workflows using only prompt engineering. In such
situations, GPT-like models might ‘hallucinate’. They might give answers that
seem correct. However, these answers could mislead designers who don’t have
the necessary domain expertise. The following examples illustrate this issue.

Prompt: Calculate the maximum load that can be supported on the clay layer at
a depth of 5 m below the ground surface due to the leaning tower of Pisa with
a diameter of 20 m.

Response: The maximum load the clay layer can support at a depth of 5 m below
the ground surface due to the leaning tower of Pisa can be calculated using
plane strain bearing capacity theory. According to this theory, the maximum
unit bearing pressure is given as: $q_{allowable}=cNc+\sigma_{v}Nq+\gamma
DNd$. $c$ = Cohesion of clay layer ($kN/m^{2}$), $\sigma_{v}=$ Overburden
Pressure (kPa), and $\gamma D=$ Density of clay layer ($kN/m^{3}$). Assuming
an average cohesion of $1\text{\,}\mathrm{kN}\text{\,}{\mathrm{m}}^{-2}$, an
average overburden pressure of 100 kPa and an average density of
$18.92\text{\,}\mathrm{kN}\text{\,}{\mathrm{m}}^{-3}$, the bearing capacity
factor can be determined using Mohr-Coulomb Failure criterion, resulting in
the following $Nc=15,Nq=5,Nd=3$. Therefore, $q_{allowable}=15\times
1+100\times 5+18.92\times 3=734.6$ kPa. Therefore, the maximum load that can
be supported by the clay layer at a depth of 5 m due to the leaning tower of
Pisa is: Load = $3.14\times(20)^{2}\times 734.6=6,037,295.4$ kN

Despite multiple trials, GPT-3 incorrectly reports a maximum allowable load as
6,037 MN. Although GPT states that it needs to use the bearing capacity
equation, the equation is incorrect missing factors such 0.5 in the last term
and incorrectly assuming $\gamma D$ as Density of the clay layer. GPT also
randomly predicts values of the bearing capacity factors resulting in an
incorrect allowable bearing capacity of 734 kPa, which also has an arithmetic
error (actual answer for that GPT calculation is 571.76 kPa). Finally, GPT
incorrectly assumes the given diameter of the tower of 20 m to be the radius,
and calculates the maximum load, again incorrect calculation, as 6,037 MN
(correct answer for that calculation is 922 MN). The maximum allowable load is
98 MN using an undrained strength of 35 kPa in the clay layer from the soil
report. Hence, GPT reports an answer that is 61 times more than the maximum
allowable load. Using GPT without proper reasoning and engineering
understanding would result in catastrophic engineering failures.

How can we enable GPT to solve complex engineering workflows? To answer the
question, we explore how humans solve problems. An extraordinary
characteristic of human intelligence is the capability to flawlessly integrate
task-oriented activities with verbal cognition (or inner speech) and
maintaining a working memory. Consider the situation of assembling a piece of
furniture. Amid any two distinct steps, we might utilize verbal reasoning to
monitor our progress (“Now that I’ve tightened all the screws, I should attach
the back panel”), to deal with unusual circumstances or modify the strategy
according to the conditions (“I’ve run out of the right screws, so I’ll have
to use these other ones and buy replacements later”), and to identify when we
need additional information (“What is the proper way to secure the drawer?
I’ll need to look that up online”). We may also take actions (read the
assembly instructions, gather tools, verify parts) to facilitate the reasoning
process and to respond to queries (“What part should I assemble next?”). This
strong interplay between “taking action” and “reasoning” allows humans to
rapidly acquire new skills and perform dependable decision-making or
reasoning, even in novel situations or when faced with uncertain information.

The groundbreaking findings by [7] on LLM ReAct suggest the potential of
merging verbal reasoning with interactive decision-making within autonomous
systems. We propose a novel strategy aimed at uniting the reasoning and
linguistic capabilities of GPT (or inner speech) with a suite of high-level
function libraries (action tools). This strategy addresses the challenges GPT
models face, such as performing intricate geotechnical calculations like
estimating bearing capacity. Instead of allowing GPT to perform these
calculations unsupervised, potentially leading to inaccuracies, we introduce
specific engineering tools. An example of such a tool is the
‘BearingCapacityTool’, which inputs soil parameters and foundation dimensions
to generate the allowable bearing pressure. By constraining the inputs, we let
GPT know that in order to use the tool, it needs to provide two inputs and the
output will be the bearing pressure. This allows GPT to string together a
proper workflow, such as using the ‘SoilReportTool’ to get soil parameters.
However, this ReAct frameworks is constrained, and only supports short-term
memory. This restriction means that knowledge from a previous step can only be
transferred to the next immediate step. We propose a novel extension to the
ReAct framework that incorporates a database to store pertinent information,
acting as a long-term memory for the LLM workflows. All LLM Action tools can
access this database, enabling them to utilize both short-term and long-term
memories in their calculations. This unique integration of short and long-term
memory, along with the LLM’s linguistic prowess, creates a unified AI
framework capable of solving complex engineering workflows. An overview of
this framework, which includes an architecture of Action - Observation -
Thought enabled by both Short and Long Term memory, is illustrated in fig. 5a.

Figure 5b shows how GPT accesses the engineering tools using the Action-
Observation-Thought framework to calculate the maximum load accurately for the
same leaning tower problem. When provided with the task of calculating the
load, GPT offers a reasoning workflow and identifies the relevant tools to use
to solve the problem. GPT takes an action with the ‘SoilReportTool’ to extract
strength parameters, makes an observation, then follows it with more tools
‘BearingCapacityTool’, ‘ShapeFactorTool’ and finally the ‘MaxLoadTool’ to
extract the final answer. The GPT reasoning workflow is shown below.

Prompt: Calculate the maximum load that can be supported on the clay layer at
a depth of 5 m below the ground surface due to the leaning tower of Pisa with
a diameter of 20 m using the soil report: pisa_report.pdf.
Response: We need to use the soil report to extract the undrained strength
parameters and then use that to calculate the bearing capacity and maximum
load. We also need to save the foundation dimension.

Action Tool: SoilReport
Action Input: pisa_report.pdf
Action: Extracting soil parameters from report: pisa_report.pdf
Interpolated strength of Clay from a depth of -5 to -17.5 m asl. Su = 35 kPa.

Observation: Su = 35 kPa

Thought: Now that we have the soil report, we can use BearingCapacity to
calculate the bearing capacity of the foundation.

Action Tool: BearingCapacity
Action Input: Su = 35 kPa and shape factor (Sc = 1.11) (from long-term
database)
Action: Calculate bearing capacity of a circular foundation with soil
parameters: ‘Su’: 35.0.

Observation: Bearing capacity = 199.689 kPa

Thought: Now that we have the bearing capacity, we can use MaxLoad to
calculate the maximum load.

Action Tool: MaxLoad
Action Input: Bearing capacity ($q_{f}=199.689$ kPa and foundation dimension
($\phi=20$ m).
Transfer load to layer at 5 m depth using 2:1 stress transfer mechanism.

Observation: Max. Load = 98022.02586093749 kN

Thought: I now know the final answer.

Final Answer: The maximum load on the clay layer at 5 m depth for the leaning
tower of Pisa with a diameter of 20 m is 98.022 MN.

Figure 5: LLM Action-Observation-Thought workflow with short-term and long-
term memory.

Using the ReAct framework GPT accurately solves the maximum load using both
short and long term memories. By developing tools to formulate workflows, GPT
will revolutionize and accelerate engineering progress by serving as a unified
input interface for complex tasks such as generating a mesh for finite element
analysis, monitoring structural health, analyzing data or designing
foundations, all requiring specialized skills.

###  3.5 The future of geotechnical engineering with LLMs

LLMs, like GPT, offer numerous benefits to geotechnical engineering. Open-
source LLMs, like LlaMA [4] and Alpaca, require less computational power and
democratize LLM access. However, it is crucial to develop context-specific
models to maximize the potential of LLMs. This process involves accessing
relevant data and constructing specialized interfaces tailored to the unique
requirements of the field. While GPT can serve as an effective reasoning
engine, it is essential to remember that context plays a significant role in
its performance. Additionally, GPT has the potential to become a natural
interface for completing complex tasks, such as data analysis and design. By
integrating GPT into geotechnical engineering workflows, professionals can
streamline their work and make informed decisions more efficiently, developing
sustainable and resilient infrastructure systems of the future.

## Acknowledgments

I thank Prof. Ellen Rathje for her insightful comments and discussions in
writing the article.

## References

  * [1] Emily M Bender, Timnit Gebru, Angelina McMillan-Major and Shmargaret Shmitchell “On the Dangers of Stochastic Parrots: Can Language Models Be Too Big?” In _Proceedings of the 2021 ACM conference on fairness, accountability, and transparency_ , 2021, pp. 610–623
  * [2] Sébastien Bubeck et al. “Sparks of artificial general intelligence: Early experiments with gpt-4” In _arXiv preprint arXiv:2303.12712_ , 2023
  * [3] Fabio Petroni et al. “How context affects language models’ factual predictions” In _arXiv preprint arXiv:2005.04611_ , 2020
  * [4] Hugo Touvron et al. “Llama: Open and efficient foundation language models” In _arXiv preprint arXiv:2302.13971_ , 2023
  * [5] Ashish Vaswani et al. “Attention is all you need” In _Advances in neural information processing systems_ 30, 2017
  * [6] Jason Wei et al. “Chain of thought prompting elicits reasoning in large language models” In _arXiv preprint arXiv:2201.11903_ , 2022
  * [7] Shunyu Yao et al. “React: Synergizing reasoning and acting in language models” In _arXiv preprint arXiv:2210.03629_ , 2022


# Bi-Static Sensing for Near-Field RIS Localization

Reza Ghazalian1, Kamran Keykhosravi2, Hui Chen2, Henk Wymeersch2, Riku Jäntti1
1Aalto University, Finland, 2Chalmers University of Technology, Sweden
E-mail: {reza.ghazalian<EMAIL_ADDRESS>{kamrank, hui.chen,
<EMAIL_ADDRESS>
###### Abstract

We address the localization of a reconfigurable intelligent surface (RIS) for
a single-input single-output multi-carrier system using bi-static sensing
between a fixed transmitter and a fixed receiver. Due to the deployment of
RISs with a large dimension, near-field (NF) scenarios are likely to occur,
especially for indoor applications, and are the focus of this work. We first
derive the Cramér-Rao bounds (CRBs) on the estimation error of the RIS
position and orientation and the time of arrival (TOA) for the path
transmitter-RIS-receiver. We propose a multi-stage low-complexity estimator
for RIS localization purposes. In this proposed estimator, we first perform a
line search to estimate the TOA. Then, we use the far-field approximation of
the NF signal model to implicitly estimate the angle of arrival and the angle
of departure at the RIS center. Finally, the RIS position and orientation
estimate are refined via a quasi-Newton method. Simulation results reveal that
the proposed estimator can attain the CRBs. We also investigate the effects of
several influential factors on the accuracy of the proposed estimator like the
RIS size, transmitted power, system bandwidth, and RIS position and
orientation.

###### Index Terms:

RIS Localization, Near-Field, CRB.

##  I Introduction

Reconfigurable intelligent surfaces (RISs) consist of controllable elements
that can alter the characteristics of an electromagnetic (EM) wave, such as
phase, amplitude, frequency, or polarization [1, 2, 3]. By properly
configuring the element coefficients, RISs provide controllable communication
channels to support communication and localization services for wireless
devices [4]. With the properties of being energy-efficient, lightweight, easy
to deploy, and compatible with the existing wireless infrastructures, RISs are
expected to be the vital enabler for the sixth generation (6G) communication
systems [5, 6].

RIS-aided communication and localization have been intensely studied in recent
years, see, e.g., [7, 8, 9, 10, 11]. In the RIS-assisted communications
context, for instance, a location-aware beamforming design in the far-field
(FF) region was proposed in [7], where authors established the geometrical
channel model based on the user’s 3D location. In the localization context,
the 3D position of a user with the aid of RIS has been estimated in [8], where
the authors have considered the FF model. For 3D user localization in the
near-field (NF) region of an RIS acting as a lens, a low-complexity estimator
has been proposed in [9]. For 2D RIS-aided user positioning, the RIS-assisted
received signal strength (RSS)-based localization algorithms have been
proposed in [10]. In [11], the linear RIS phase profile has been designed for
both communications and positioning systems for the FF region in a 2D
scenario. According to the recent works in RIS-aided localization systems, it
can be realized that the location and orientation of the RIS need to be known.
That is because RIS is one of the reference points in these systems [1]. In
the context of RIS-aided communication systems, knowing RIS location can
facilitate the location-based RIS phase profile optimization [12]. However,
the _location and the orientation of the RIS were assumed to be known_ in all
previous studies, which is non-trivial to obtain in practice. The problem of
RIS localization (i.e., determining the position of an RIS) can be seen as a
bi-static sensing problem, and was previously tackled in [13] under FF
conditions. However, the RIS orientation could not be determined, given the
time of arrival (TOA) measurements. When the RIS is large, devices in most
RIS-aided systems will be in the NF, which makes the RIS localization problem
more challenging than in [13].

In this paper, we aim at estimating the location and orientation of an RIS in
a NF scenario with a single-antenna receiver (RX), where a multi-carrier (MC)
single-antenna transmitter (TX) sends signals to the RX via the linear RIS.
Our main contributions are as follows: (i) we formulate the RIS localization
problem in the NF region for a single-input and single-output (SISO)-MC system
with a blocked line of sight (LOS) path, (ii) we derive a Fisher information
matrix (FIM) for RIS location and orientation estimate, (iii) we propose a
multi-stage low-complexity RIS localization algorithm and evaluate it by
comparing with the theoretical error bound; Moreover, we use FF approximation
of the NF observed signals to earn information about the angle-of-arrival
(AOA) and the angle-of-departure (AOD) at the RIS.

Figure 1: System setup, comprising single-antenna TX and RX, as well as an RIS
with unknown location and orientation. $\tilde{\theta}_{\text{tx}}$ and
$\tilde{\theta}_{\text{rx}}$ represent the AOA and AOD at the RIS center when
$\alpha=0$, respectively.

#### Notation

Vectors and matrices are indicated by lowercase and uppercase bold letters,
respectively. The element in the $i$th row and $j$th column of matrix
$\mathbf{A}$ is specified by $\left[\mathbf{A}\right]_{i,j}$. The sub-index
$i:j$ determines all the elements between $i$ and $j$. The complex conjugate,
Hermitian, and transpose operators are represented by $\left(.\right)^{*}$,
$\left(.\right)^{\mathsf{H}}$, and $\left(.\right)^{\top}$, respectively.
$\|.\|$ calculates the norm of vectors or Frobenius norm of matrices. By
$\odot$, we indicate element-wise product. $\jmath=\sqrt{-1}$ and
$\mathbf{1}_{K}$ is a column vector comprising all ones with length $K$. The
function $\text{atan2}(y,x)$ is the four-quadrant inverse tangent function.

##  II System Model

In this section, we describe the considered wireless system for 2D111Our study
is limited to 2D, given the inherent difficulty of RIS localization. The 3D
problem is left for future work. RIS localization and elaborate the signal
model.

###  II-A System Setup

We consider a wireless system consisting of a single-antenna TX, a single-
antenna RX, and an $M$-elements linear RIS in a two-dimensional Cartesian
coordinate system (global reference) as shown in Fig. 1. The RX and TX are
located at NF region222NF region, which is usually defined as the range
smaller than the Fraunhofer distance $d_{F}=2D^{2}/\lambda$ [14]. Here, $D$
and $\lambda$ are the dimension of the RIS and the wavelength at the carrier
frequency, respectively. of the RIS with the LOS path blocked.333This is not a
limiting assumption, but only invoked for notational convenience. When the LOS
path is present, it can be separated from the RIS path as in [13]. When the TX
and RX are synchronized, the LOS path does not convey any information. We also
assume that TX and RX are synchronized (i.e., known clock bias) with known
position ($\mathbf{p}_{\text{tx}}$ and $\mathbf{p}_{\text{rx}}$) (e.g., RX as
calibration agent).

With a rotation angle $\alpha$, the global position of $m$th RIS element can
be expressed as a function of the RIS center $\mathbf{p}_{\text{ris}}$:

$\mathbf{p}_{m}=\mathbf{p}_{\text{ris}}+\mathbf{R}_{\alpha}\left(m-\frac{M-1}{2}\right)\Delta,$
(1)

where $\Delta$ is the inter-element space (e.g., ${\lambda}/{2}$) and
$\mathbf{R}_{\alpha}$ is the 2D rotation matrix as

$\mathbf{R}_{\alpha}=\begin{bmatrix}\cos{(\,\alpha)\,}&-\sin{(\,\alpha})\,\\\
\sin{(\,\alpha)\,}&\cos{(\,\alpha})\,\\\ \end{bmatrix}.$ (2)

The aim of this work is to estimate the orientation ($\alpha$) and the center
position of the RIS ($\mathbf{p}_{\text{ris}}$).

###  II-B Signal Model

We consider a MC scenario, where TX transmits $T$ orthogonal frequency
division multiplexing (OFDM) symbols over time with $N_{c}$ sub-carries for
RIS localization purposes. Without loss of generality, we assume that all the
transmitted symbols are equal to one. It is also assumed that RIS changes the
phase of the incident wave through a random profile over time (denoted by
$t$). The RIS phase profile at time $t$ is shown by the vector
$\boldsymbol{\gamma}(t)\in\mathbb{C}^{M\times 1}$, where
$|[\boldsymbol{\gamma}(t)]_{m}|=1$ for all $m$. Based on these assumptions, we
can write the received signal over time and sub-carrier frequencies as the
following $N_{c}\times T$ matrix

$\displaystyle\mathbf{Y}=\rho
e^{\jmath\phi}\sqrt{P_{t}}\mathbf{d}\left(\tau\right)\mathbf{b}^{\top}\boldsymbol{\Gamma}+\mathbf{W},$
(3)

where $\rho e^{\jmath\phi}$ is an unknown complex gain, modeled with
$\phi\sim\mathcal{U}[0,2\pi)$ and [15]

$\displaystyle\rho=\sqrt{\frac{G_{t}G_{r}A\lambda^{2}}{64\pi^{3}}}\frac{1}{\|\mathbf{p}_{\text{tx}}-\mathbf{p}_{\text{ris}}\|\|\mathbf{p}_{\text{rx}}-\mathbf{p}_{\text{ris}}\|},$
(4)

where $A=\frac{\Delta^{2}}{4}$ denotes the RIS element’s surface area and
$G_{t}$ and $G_{r}$ are the antenna gain at the TX and RX, respectively.
$P_{t}$ denotes the transmitted power. The matrix $\mathbf{\Gamma}$ contains
all of the RIS phase profiles over time, that is
$\mathbf{\Gamma}=\begin{bmatrix}\boldsymbol{\gamma}(0),\dots,\boldsymbol{\gamma}(T)\end{bmatrix}$,
$\mathbf{W}\in\mathbb{C}^{N_{c}\times T}$ is the noise matrix containing zero-
mean circularly-symmetric independent and identically distributed Gaussian
elements with variance $\sigma^{2}$. We consider the effect of the delay on
the signal with the vector $\mathbf{d}(\tau)\in\mathbb{C}^{N_{c}\times 1}$,
which is defined as

$\mathbf{d}(\tau)=[\,1,e^{-\jmath 2\pi\tau\Delta f},\dots,e^{-\jmath
2\pi\tau(N_{c}-1)\Delta f}]\,^{\top},$ (5)

where $\Delta f$ is sub-carrier spacing and
$\tau=\|\mathbf{p}_{\text{tx}}-\mathbf{p}_{\text{ris}}\|/c+\|\mathbf{p}_{\text{rx}}-\mathbf{p}_{\text{ris}}\|/c$
for speed of light $c$. Finally,
$\mathbf{b}=\mathbf{a}_{t}\odot\mathbf{a}_{r}$, where the vectors
$\mathbf{a}_{t}$ and $\mathbf{a}_{t}$ are the NF array response vectors from
the RIS to the TX and RX, respectively. The $m$th elements of these vectors
are defined as [16]

$[\mathbf{a}_{t}]_{m}=e^{-\jmath\frac{2\pi}{\lambda}\left(R^{t}_{m}-R^{t}_{o}\right)},$
(6a)
$[\mathbf{a}_{r}]_{m}=e^{-\jmath\frac{2\pi}{\lambda}\left(R^{r}_{m}-R^{r}_{o}\right)},$
(6b)

where $R^{t}_{m}=\|\mathbf{p}_{m}-\mathbf{p}_{\text{tx}}\|$,
$R^{r}_{m}=\|\mathbf{p}_{m}-\mathbf{p}_{\text{rx}}\|$,
$R^{t}_{o}=\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{tx}}\|$,
$R^{r}_{o}=\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}}\|$. Here, we
assume that the wavelength remains relatively constant over the transmission
bandwidth (BW), so that beam squint effects can be ignored.

##  III Fisher Information Analysis

Introducing the noiseless part of observation given in (3) as
$\boldsymbol{\mu}=\rho
e^{\jmath\phi}\sqrt{P_{t}}\mathbf{d}(\,\tau)\,\mathbf{b}^{\top}\mathbf{\Gamma}$
, and the $6\times 1$ vector of unknowns
$\boldsymbol{\eta}=[\rho\quad\phi\quad\tau\quad\mathbf{p}_{\text{ris}}^{\top}\quad\alpha]^{\top}$,
and
$\boldsymbol{\zeta}=\left[\rho\quad\phi\quad\mathbf{p}_{\text{ris}}^{\top}\quad\alpha\right]^{\top}$,
the FIM $\mathbf{J}(\boldsymbol{\eta})\in\mathbb{R}^{6\times 6}$ is defined as
[17, Sec. 3.9]

$\displaystyle\mathbf{J}(\boldsymbol{\eta})=\frac{2}{\sigma^{2}}\sum_{t=1}^{T}\sum_{n_{c}=1}^{N_{c}}\Re\Bigg{\\{}\frac{\partial[\boldsymbol{\mu}]_{t,n_{c}}}{\partial\boldsymbol{\eta}}\left(\frac{\partial[\boldsymbol{\mu}]_{t,n_{c}}}{\partial\boldsymbol{\eta}}\,\right)^{\mathsf{H}}\,\Bigg{\\}}.$
(7)

Using (7), we can calculate the Cramér-Rao bound (CRB) of the $\tau$, i.e.,
the time error bound (TEB) as

$\displaystyle\sqrt{\mathbb{E}\left[\left(\tau-\hat{\tau}\right)^{2}\right]}\geq\text{TEB}\triangleq\sqrt{\left[\mathbf{J}(\boldsymbol{\eta})^{-1}\right]_{3,3}}\quad,$
(8)

where $\hat{\tau}$ is the estimate of the true $\tau$. Since there is a
relationship between $\tau$ and $\mathbf{p}_{\text{ris}}$, the FIM of position
and orientation $\mathbf{J}\left(\boldsymbol{\zeta}\right)$ can be derived as
$\mathbf{J}\left(\boldsymbol{\zeta}\right)=\mathbf{T}^{\top}\mathbf{J}\left(\boldsymbol{\eta}\right)\mathbf{T}$,
where $\mathbf{T}\in\mathbb{R}^{6\times 5}$ is a Jacobian matrix [17, Eq.
(3.30)]. Its element on the $\ell$th row and the $q$th column is defined

$\displaystyle\left[\mathbf{T}\right]_{\ell,q}=\frac{\partial\left[\boldsymbol{\eta}\right]_{\ell}}{\partial\left[\boldsymbol{\zeta}\right]_{q}}\quad.$
(9)

Considering (7) and (9), the position error bound (PEB) and orientation (EOB)
can be written as

$\displaystyle\sqrt{\mathbb{E}\left[\|\mathbf{p}_{\text{ris}}-\hat{\mathbf{p}}_{\text{ris}}\|^{2}\right]}\geq\text{PEB}\triangleq\sqrt{\mathrm{tr}(\left[\mathbf{J}(\boldsymbol{\zeta})^{-1}\right]_{3:4,3:4})}\quad,$
(10a)
$\displaystyle\sqrt{\mathbb{E}\left[\left(\alpha-\hat{\alpha}\right)^{2}\right]}\geq\text{OEB}\triangleq\sqrt{\left[\mathbf{J}(\boldsymbol{\zeta})^{-1}\right]_{5,5}}\quad,$
(10b)

where $\hat{\mathbf{p}}_{\text{ris}}$ and $\hat{\alpha}$ are the estimates of
the true RIS position and orientation, respectively. The matrices $\mathbf{J}$
and $\mathbf{T}$ are calculated in the Appendix.

##  IV RIS Location and Orientation Estimator

###  IV-A Maximum Likelihood Estimator

Based on the observation given in (3), the Maximum Likelihood (ML) estimator
is defined as

$\displaystyle\left[\hat{g}_{r},\hat{\alpha},\hat{\mathbf{p}}_{\text{ris}}\right]$
$\displaystyle=\arg\max_{\alpha,\mathbf{p}_{\text{ris}}}f\left(\mathbf{Y}|g_{r},\alpha,\mathbf{p}_{\text{ris}}\right)$
(11)
$\displaystyle=\arg\min_{g_{r},\alpha,\mathbf{p}_{\text{ris}}}\|\mathbf{Y}-g_{r}\sqrt{P_{t}}\mathbf{d}\left(\tau\right)\mathbf{b}^{\top}\boldsymbol{\Gamma}\|^{2},$

where $g_{r}=\rho e^{\jmath\phi}$. To solve (11), we can use any gradient
descent method (e.g., Newton method). However, the function described in (11)
is non-convex and has many local optima around the global optimum. Thus
finding a suitable initial estimate for gradient descent methods is
challenging. To tackle this issue, we propose a low-complexity estimator,
which can provide such an initial estimate.

###  IV-B Low Complexity Estimator

In this section, we develop a low-complexity estimator for RIS orientation and
position estimation. To this end, first, we estimate $\tau$. Then we use the
FF approximation of the NF to gain information about the AOA and the AOD at
the RIS. Based on these estimations, we can find the RIS position and
orientation through a line search. Finally, we adopt this estimate as an
initial guess for the quasi-Newton algorithm to solve the problem (11). We
indeed use the quasi-Newton algorithm to refine the initial estimate.

####  IV-B1 Estimation of $\tau$

For $\tau$ estimation, we take the IFFT with length $N_{F}$ of the columns of
the matrix $\mathbf{Y}$, yielding $\mathbf{Z}=\mathbf{F}\mathbf{Y}$. We then
obtain a coarse estimate of the TOA as

$\displaystyle\tilde{k}=\arg\max_{k}\|\mathbf{g}^{\top}_{k}\mathbf{Z}\|,$ (12)

where $\mathbf{g}_{k}$ is vector comprising $N_{F}$ zeros, except with one in
the $k$th entry. Accordingly, one can obtain a refined estimate $\tau$ as
follows [13]:

$\hat{\tau}=\frac{\tilde{k}}{N_{F}\Delta f}-\tilde{\delta},$ (13)

where

$\displaystyle\tilde{\delta}=\arg\max_{\delta\in\left[0,1/\left(N_{F}\Delta
f\right)\right]}\|\mathbf{g}^{\top}\mathbf{Z}\
\odot\left(\mathbf{d}\left(\delta\right)\mathbf{1}_{T}^{\top}\right)\|.$ (14)

####  IV-B2 Estimating AOA and AOD at the center of RIS

We first remove the effect of $\tau$ from $\mathbf{Y}$, i.e.,

$\displaystyle\mathbf{Y}_{r}$
$\displaystyle=\mathbf{Y}\odot\left(\mathbf{d}\left(-\tau\right)\mathbf{1}_{T}^{\top}\right)$
(15)
$\displaystyle=g_{r}\sqrt{P_{t}}\mathbf{1}_{N_{c}}\mathbf{b}^{\top}\boldsymbol{\Gamma}+\mathbf{W}_{d},$

where
$\mathbf{W}_{d}\triangleq\mathbf{W}\odot\left(\mathbf{d}\left(-\tau\right)\mathbf{1}_{T}^{\top}\right)$.
Next, we sum the signal $\mathbf{Y}_{r}$ across the sub-carriers to obtain

$\displaystyle\mathbf{y}_{r}$
$\displaystyle=\mathbf{Y}_{r}^{\top}\mathbf{1}_{N_{c}}\Rightarrow\mathbf{y}_{r}=N_{c}g_{r}\sqrt{P_{t}}\boldsymbol{\Gamma}^{\top}\mathbf{b}+\mathbf{w}_{t},$
(16)

where $\mathbf{w}_{t}\triangleq\mathbf{W}^{\top}\mathbf{1}_{N_{c}}$. In the NF
regime, AOA and AOD can be defined for each array element. However, estimating
all these angles for all elements is not possible, as they are unobservable.
Hence, we only estimate the AOA and AOD at the RIS center using a FF
approximation. Based on (6), we can write FF approximation of $\mathbf{b}$ as

$\displaystyle[\mathbf{b}]_{m}=e^{-\jmath\frac{2\pi}{\lambda}m\Delta\left(\sin{\left(\theta_{\text{tx}}\right)}+\sin{\left(\theta_{\text{rx}}\right)}\right)},$
(17)

where $\theta_{\text{tx}}\triangleq\tilde{\theta}_{\text{tx}}+\alpha$ and
$\theta_{\text{rx}}\triangleq\tilde{\theta}_{\text{rx}}+\alpha$ are AOA and
AOD at the center of the RIS, respectively (see Fig. 1). Given (17), we can
estimate

$\displaystyle\omega\triangleq\sin{\left(\tilde{\theta}_{\text{tx}}+\alpha\right)}+\sin{\left(\tilde{\theta}_{\text{rx}}+\alpha\right)}.$
(18)

Thus, considering (16)–(18), we can rewrite (16) as

$\displaystyle\mathbf{y}_{r}=N_{c}g_{r}\sqrt{P_{t}}\mathbf{\Gamma}^{\top}\mathbf{b\left(\omega\right)}+\mathbf{w}_{t}.$
(19)

Here, for each value of $\omega$, the value of $g_{r}$ can be estimated by
$\hat{g}_{r}(\omega)=\mathbf{b}^{\mathsf{H}}\left(\omega\right)\boldsymbol{\Gamma}^{*}\mathbf{y}_{r}/(N_{c}\sqrt{P_{t}}\|\mathbf{b}^{\mathsf{H}}\left(\omega\right)\boldsymbol{\Gamma}^{*}\|^{2})$.
Hence, the estimation for $\omega$ can be obtained as

$\displaystyle\hat{\omega}=\arg\min_{\omega}\|\mathbf{y}_{r}-N_{c}\sqrt{P_{t}}\hat{g}_{r}\left(\omega\right)\mathbf{\Gamma}^{\top}\mathbf{b}\left(\omega\right)\|^{2},$
(20)

which can be solved via a simple 1D line search.

####  IV-B3 RIS Position and Orientation Estimation

(a)

(b)

(c)

Figure 2: The evaluation of the RMSE of the $\mathbf{p}_{\text{ris}}$,
$\mathbf{\alpha}$, and $\mathbf{\tau}$ versus $P_{t}$. The simulation
parameters are given in Table I: (a) RMSE of the RIS position (meter) and PEB
versus $P_{t}$; (b) RMSE of the RIS orientation (rad) and OEB versus $P_{t}$ ;
(c) RMSE of the $\tau$ (nsec) and TEB versus $P_{t}$.

Based on the $\hat{\tau}$, the RIS is constrained to lie on the ellipse
defined by
$\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{tx}}\|+\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}}\|=c\hat{\tau}$.
Thus a possible solution of RIS position can be found by parameterizing this
ellipse a function of $\nu\in[0,2\pi)$ as

$\displaystyle\mathbf{p}_{\text{ris}}\left(\nu\right)=\mathbf{R}_{\beta}\mathbf{A}_{\nu}\mathbf{e}+\mathbf{c}_{e}.$
(21)

Here, $\mathbf{c}_{e}=(\mathbf{p}_{\text{rx}}+\mathbf{p}_{\text{tx}})/2$ is
the center of the ellipse, $\mathbf{R}_{\beta}$ is a rotation matrix with
angle $\beta$, where $\beta$ is a angle between positive $x$-axis and the line
between TX and RX.
$\mathbf{A}_{\nu}\triangleq\text{diag}(\left[\cos{\nu},\sin{\nu}\right]^{\top})$
and
$\mathbf{e}=[c\hat{\tau}/2,\sqrt{a^{2}-\|\mathbf{p}_{\text{tx}}-\mathbf{c}_{e}\|^{2}}]^{\top}$.
Furthermore, we can formulate the $\tilde{\theta}_{\text{tx}}$ and
$\tilde{\theta}_{\text{rx}}$ a function of $\nu$ as

$\displaystyle\tilde{\theta}_{\text{rx}}\left(\nu\right)$
$\displaystyle=\text{atan2}\left(y_{\text{n}_{\text{r}}},x_{\text{n}_{\text{r}}}\right)-\frac{\pi}{2},$
(22) $\displaystyle\tilde{\theta}_{\text{tx}}\left(\nu\right)$
$\displaystyle=\text{atan2}\left(y_{\text{n}_{\text{t}}},x_{\text{n}_{\text{t}}}\right)-\frac{\pi}{2},$
(23)

where
$[x_{\text{n}_{\text{r}}},y_{\text{n}_{\text{r}}}]^{\top}=(\mathbf{p}_{\text{rx}}-\mathbf{p}_{\text{ris}}(\nu))$
and
$[x_{\text{n}_{\text{t}}},y_{\text{n}_{\text{t}}}]^{\top}=(\mathbf{p}_{\text{tx}}-\mathbf{p}_{\text{ris}}(\nu))$.

Based on the definition of $\omega$, $\hat{\omega}$, (22) and (23), $\alpha$
is also obtained as a function of $\nu$

$\displaystyle\alpha\left(\nu\right)=\arcsin\left(\frac{\hat{\omega}}{2\cos{\theta_{\text{tr}}^{-}}}\right)-\theta_{\text{tr}}^{+}\quad,$
(24)

where
$\theta_{\text{tr}}^{+}\triangleq(\tilde{\theta}_{\text{tx}}+\tilde{\theta}_{\text{rx}})/2$
and
$\theta_{\text{tr}}^{-}\triangleq(\tilde{\theta}_{\text{tx}}-\tilde{\theta}_{\text{rx}})/2$.
Hence, we can estimate $\nu$ as

$\displaystyle\hat{\nu}$
$\displaystyle=\arg\min_{\nu}\|\mathbf{Y}-\hat{g}_{r}\left(\nu\right)\sqrt{P_{t}}\mathbf{d}\left(\hat{\tau}\right)\mathbf{b\left(\mathbf{p}_{\text{ris}}\left(\nu\right),\alpha\left(\nu\right)\right)}^{\top}\boldsymbol{\Gamma}\|^{2},$
(25)

where
$\hat{g}_{r}(\nu)=\mathbf{h}^{\mathsf{H}}\left(\nu\right)\mathbf{y}/(\sqrt{P_{t}}\|\mathbf{h}\left(\nu\right)\|^{2})$,
$\mathbf{y}=\mathbf{Y}^{\top}\mathbf{1}_{N_{c}}$ and
$\mathbf{h}(\nu)\triangleq(\mathbf{d}(\hat{\tau})\mathbf{b(\mathbf{p}_{\text{ris}}(\nu),\alpha(\nu))}^{\top}\boldsymbol{\Gamma})^{\top}\mathbf{1}_{N_{c}}$,
$\mathbf{b(\mathbf{p}_{\text{ris}}(\nu),\alpha(\nu))}$ is defined based on
(6a) and (6b). Note that (25) can be solved through a line search. Then, the
position and orientation of the RIS can be obtained via (21)–(25). This
solution is applied as the initial guess for quasi-Newton method to solve
(11).

In summary, the RIS location can be recovered by several 1D line searches,
which leads to very low complexity.

##  V Simulation Results

In this section, we evaluate the proposed estimator and compare it to the
corresponding bound. To do so, we compare the root mean square error (RMSE) of
the estimated parameters (e.g, RIS position) with the derived CRB.

###  V-A Simulation Parameters

To evaluate the performance of the proposed estimator, we average over $500$
noise realizations. We assume that the global coordinate system is set to be
aligned with RX coordinate system. The phase profile at the RIS elements is
randomly drawn from the uniform distribution $[0,2\pi)$. The number of RIS
elements is set to $M=64$. The rest of the simulation parameters are given in
Table I.

TABLE I: Simulation Parameters Parameter | Symbol | Value
---|---|---
Wavelength | $\lambda$ | $1~{}\text{cm}$
RIS element spacing | $\Delta$ | $0.5~{}\text{cm}$
Light speed | $c$ | $3\times 10^{8}~{}\text{m}/\text{sec}$
Number of sub-carriers | $N_{c}$ | $500$
TX antenna gain | $G_{t}$ | $2$
RX antenna gain | $G_{r}$ | $2$
Number of transmissions | $T$ | $50$
Sub-carrier spacing | $\Delta f$ | $120~{}\text{kHz}$
Noise PSD | $N_{0}$ | $-174~{}\text{dBm/Hz}$
RX’s noise figure | $n_{f}$ | $8~{}\text{dB}$
Noise variance | $\sigma^{2}=n_{f}N_{0}N_{c}\Delta f$ | $-88~{}\text{dBm}$
IFFT Size | $N_{F}$ | $4096$
RX position | $\mathbf{p}_{\text{rx}}$ | $\left[2,2\right]^{\top}$
TX position | $\mathbf{p}_{\text{tx}}$ | $\left[0,2\right]^{\top}$
RIS position | $\mathbf{p}_{\text{ris}}$ | $\left[0,0\right]^{\top}$
RIS orientation | $\alpha$ | $\pi/6~{}\text{rad}$

###  V-B Results and Discussion

In Fig. 2 , we study the effect of the transmitted power on the performance of
the proposed estimator. As expected, the RMSE of the
$\mathbf{p}_{\text{ris}}$, $\alpha$, and $\tau$ are decreasing functions of
the transmitted power $P_{t}$ and attains the CRB when $P_{t}\geq 0$ dBm.

Next, we investigate the effect of the system BW on the RMSE of the RIS
localization, shown in Fig. 3. To this end, we increase the number of sub-
carriers while the frequency spacing and transmit power are fixed. From Fig.
3(c), it is clear that as the BW increases, the RMSE of $\tau$ decreases.
Regrading the proposed method, $\tau$ estimation plays a key role on the other
stages of the proposed estimator. The more accurate $\tau$ estimate, the more
accurate $\omega$ and RIS positioning. This can be seen in Fig. 3(a) and Fig.
3(b): with better $\tau$ estimation, the ellipse equation (21) becomes more
accurate. Since the RIS location is on this ellipse, the RIS positioning
accuracy will be enhanced. On the other hand, the accuracy of the $\omega$
estimate directly depends on the $\tau$ estimate (see (15)–(20)). The $\omega$
estimate provides information about the RIS orientation. Therefore, with finer
$\tau$ and $\omega$ estimation, the estimator achieves higher accuracy in RIS
location and orientation estimate. This is the effect of the BW increase. For
bandwidths larger than 100 MHz, we observe that the performance is limited by
the the estimation of $\omega$.

(a)

(b)

(c)

Figure 3: The comparison of the RMSE of the $\mathbf{p}_{\text{ris}}$ and
$\mathbf{\alpha}$ with CRB with respect to the system BW.
$P_{t}=10~{}\text{dBm}$, and the rest of simulation parameters are given in
Table I: (a) RMSE of the RIS position (meter) and PEB versus BW; (b) RMSE of
the RIS orientation (rad) and OEB versus BW; (c) RMSE of $\tau$ (nsec) and TEB
versus BW.

Now, we analyze the impact of the RIS size on the RIS localization as shown in
Fig. 4. As it is shown in Fig. 4(a) and Fig. 4(b), the RMSE of the RIS
localization reduces as the number of the RIS elements grows. One can also see
that the proposed method can localize RIS with sub-centimeter accuracy when
the number of RIS elements reaches $128$. This observation is due to two
compounding effects: a larger RIS size leads to higher SNR, and the larger RIS
size leads to a more pronounced NF effect.

(a)

(b)

Figure 4: The comparison of the RMSE of the $\mathbf{p}_{\text{ris}}$ and
$\mathbf{\alpha}$ with CRB with respect to the RIS size. $P_{t}=10\text{dBm}$,
and the rest of the simulation parameters are given in Table I: (a) RMSE of
the RIS position (meter) and PEB versus the RIS size; (b) RMSE of the RIS
orientation (rad) and OEB versus the RIS size.

Finally, we assess the RIS localization coverage and performance through
contour plots of PEB and OEB. Fig. 5 illustrates the contour plots of PEB and
OEB when the $x$ and $y$ coordinates of the RIS are varied while its
orientation is fixed. To evaluate the effect of the RIS orientation, we depict
contour plots for two cases: (i) $\alpha=0$ (rad); (ii) $\alpha=\pi/6$ (rad).
One can see, in general, that high accurate location and orientation can be
obtained when the RIS is close to the RX or TX, due to high SNR. However,
several blind areas exist (e.g., $[-1,-0.5]^{\top}$ and $[-6,-1]^{\top}$ in
Fig. 5(a)). For $\alpha=0$ (rad), both location and orientation estimate
degrade as the RIS approaches the bisector of the line segment between TX and
RX. Furthermore, we can also see two other areas where RIS localization has a
high error. When the RIS rotates (Fig. 5(c) and Fig. 5(d)) the mentioned lines
rotates accordingly.

(a)

(b)

(c)

(d)

Figure 5: The contour plots of PEB and OEB versus the RIS location for one
random RIS phase profile and $P_{t}=10\,\text{dBm}$. TX and RX positions are
marked by red diamond and black square, respectively: (a) $10\log
10(\text{PEB})$ (meter), $\alpha=0$ (rad);(b) $10\log 10(\text{OEB})$ (rad),
$\alpha=0$ (rad);(c) $10\log 10(\text{PEB})$ (meter), $\alpha=\pi/6$ (rad);
(d) $10\log 10(\text{OEB})$ (rad), $\alpha=\pi/6$.

##  VI Conclusion

In this work, we have proposed a low-complexity bi-static sensing estimator to
determine the location and orientation of RIS in NF with the help of a single-
antenna transmitter and receiver. We have also derived the CRBs on the RIS
position and orientation estimate error. To show the efficiency of the
estimator, the RMSE of the estimations has been compared with the CRBs.
Simulation results have verified the estimator’s accuracy such that its RMSE
attains the CRB. To gain insight into the RIS localization problem, we
investigated the influential factors like transmitted power, the number of RIS
elements, BW, and RIS position and orientation concerning the TX/RX. Results
show that increasing the value of these factors enhances the RIS localization
accuracy. In future work, we would like to study 3D RIS localization.

## Acknowledgment

This work has been funded in part by Academy of Finland Profi-5 (n:o 326346)
and ULTRA (n:o 328215) projects, and the EU H2020 RISE-6G project under grant
101017011.

## Appendix

###  VI-A Derivation of $\mathbf{J}$

In this appendix, we drive the elements of the matrix $\mathbf{J}$. To this
end, we derive
$\partial[\boldsymbol{\mu}]_{t,n_{c}}/\partial\boldsymbol{\eta}$ according to
(7) as
$\partial[\boldsymbol{\mu}]_{t,n_{c}}/\partial\rho=e^{\jmath\phi}\sqrt{P_{t}}e^{-\jmath
2\pi\tau n_{c}\Delta f}\mathbf{b}^{\top}\boldsymbol{\gamma}(t)$,
$\partial[\boldsymbol{\mu}]_{t,n_{c}}/\partial\phi=\jmath\rho
e^{\jmath\phi}\sqrt{P_{t}}e^{-\jmath 2\pi\tau n_{c}\Delta
f}\mathbf{b}^{\top}\boldsymbol{\gamma}(t)$,
$\partial[\boldsymbol{\mu}]_{t,n_{c}}/\partial\boldsymbol{\tau}=-\jmath 2\pi
n_{c}\Delta f\rho e^{\jmath\phi}\sqrt{P_{t}}e^{-\jmath 2\pi\tau n_{c}\Delta
f}\mathbf{b}^{\top}\boldsymbol{\gamma}(t)$, and

$\displaystyle\frac{\partial[\boldsymbol{\mu}]_{t,n_{c}}}{\partial\boldsymbol{p}_{\text{ris}}}=\rho
e^{\jmath\phi}\sqrt{P_{t}}e^{-\jmath 2\pi\tau n_{c}\Delta
f}\mathbf{E}^{\top}\boldsymbol{\gamma}(t)$ (26)
$\displaystyle\frac{\partial[\boldsymbol{\mu}]_{t,n_{c}}}{\partial\alpha}=\rho
e^{\jmath\phi}\sqrt{P_{t}}e^{-\jmath 2\pi\tau n_{c}\Delta
f}\mathbf{v}^{\top}\boldsymbol{\gamma}(t)$ (27)

where
$\mathbf{E}\triangleq\partial\mathbf{b}/\partial\mathbf{p}_{\text{ris}}$,$\mathbf{v}\triangleq\partial\mathbf{b}/\partial\alpha$,
and we have

$\displaystyle\left[\mathbf{E}\right]_{:,m}=-\jmath\frac{2\pi}{\lambda}\left(\boldsymbol{\kappa}_{m}^{t}+\boldsymbol{\kappa}_{m}^{r}-\boldsymbol{\kappa}_{o}^{t}-\boldsymbol{\kappa}_{o}^{r}\right)\left[\mathbf{b}\right]_{m}$
(28)

where
$\boldsymbol{\kappa}_{m}^{t}\triangleq(((m-(M-1)/2)\Delta-\mathbf{p}_{\text{tx}})/R_{m}^{t}$,
$\boldsymbol{\kappa}_{m}^{r}\triangleq(((m-(M-1)/2)\Delta-\mathbf{p}_{\text{rx}})/R_{m}^{r}$,
$\boldsymbol{\kappa}_{o}^{t}\triangleq(\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}})/R_{o}^{t}$,
and
$\boldsymbol{\kappa}_{o}^{r}\triangleq(\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}})/R_{o}^{r}$,

$\displaystyle\left[\mathbf{v}\right]_{m}=-\jmath\frac{2\pi}{\lambda}\left(\boldsymbol{\psi}_{m}^{t}+\boldsymbol{\psi}_{m}^{r}\right)\left[\mathbf{b}\right]_{m},$
(29)

where
$\boldsymbol{\psi}_{m}^{t}\triangleq((\mathbf{p}_{m}-\mathbf{p}_{\text{tx}})/R_{m}^{t})\mathbf{R}^{\prime}_{\alpha}(((m-(M-1)/2)\Delta)$,
$\boldsymbol{\psi}_{m}^{r}\triangleq((\mathbf{p}_{m}-\mathbf{p}_{\text{rx}})/R_{m}^{r})\mathbf{R}^{\prime}_{\alpha}(((m-(M-1)/2)\Delta)$,
and
$\mathbf{R}^{\prime}_{\alpha}\triangleq\partial\mathbf{R}_{\alpha}/\partial\alpha$.

###  VI-B Derivation of the Jacobian Matrix $\mathbf{T}$

Based on (9) and the relation between $\tau$, $\phi$, $\rho$ and
$\mathbf{p}_{\text{ris}}$, the Jacobian matrix elements can be written as
$\left[\mathbf{T}\right]_{1:2,1:2}=\partial\left[\rho,\phi\right]^{\top}/\partial\left[\rho,\phi\right]^{\top}=\mathbf{I}_{2\times
2}$, and

$\displaystyle\left[\mathbf{T}\right]_{3,3:4}=\frac{\partial\tau}{\partial\mathbf{p}_{\text{ris}}}=\frac{1}{c}\left(\mathbf{d}_{tx}^{\top}+\mathbf{d}_{rx}^{\top}\right)$
(30)

where
$\mathbf{d}_{tx}=(\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{tx}})/\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{tx}}\|$
and
$\mathbf{d}_{rx}=(\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}})/\|\mathbf{p}_{\text{ris}}-\mathbf{p}_{\text{rx}}\|$.
Besides, we have
$\left[\mathbf{T}\right]_{4:5,3:4}=\partial\mathbf{p}_{\text{ris}}/\partial\mathbf{p}_{\text{ris}}=\mathbf{I}_{2\times
2}$ and $\left[\mathbf{T}\right]_{6,5}=\partial\alpha/\partial\alpha=1$. The
other elements of the matrix $\mathbf{T}$ are zero.

## References

  * [1] H. Wymeersch, J. He, B. Denis, A. Clemente, and M. Juntti, “Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions,” _IEEE Vehicular Technology Magazine_ , vol. 15, no. 4, pp. 52–61, 2020.
  * [2] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE access_ , vol. 7, pp. 116 753–116 773, 2019.
  * [3] J. Yuan, E. De Carvalho, R. J. Williams, E. Björnson, and P. Popovski, “Frequency-mixing intelligent reflecting surfaces for nonlinear wireless propagation,” _IEEE Wireless Communications Letters_ , 2021.
  * [4] E. C. Strinati, G. C. Alexandropoulos, H. Wymeersch, B. Denis, V. Sciancalepore, R. D’Errico, A. Clemente, D.-T. Phan-Huy, E. De Carvalho, and P. Popovski, “Reconfigurable, intelligent, and sustainable wireless environments for 6G smart connectivity,” _IEEE Communications Magazine_ , vol. 59, no. 10, pp. 99–105, 2021.
  * [5] E. Basar, “Reconfigurable intelligent surface-based index modulation: A new beyond mimo paradigm for 6G,” _IEEE Transactions on Communications_ , vol. 68, no. 5, pp. 3187–3196, 2020.
  * [6] H. Chen, H. Sarieddeen, T. Ballal, H. Wymeersch, M.-S. Alouini, and T. Y. Al-Naffouri, “A tutorial on terahertz-band localization for 6G communication systems,” _Accepted for publication in IEEE Commun. Surveys Tuts. arXiv preprint arXiv:2110.08581_ , 2022.
  * [7] Z. Xing, R. Wang, X. Yuan, and J. Wu, “Location-aware beamforming design for reconfigurable intelligent surface aided communication system,” in _2021 IEEE/CIC International Conference on Communications in China (ICCC)_ , 2021.
  * [8] K. Keykhosravi, M. F. Keskin, G. Seco-Granados, and H. Wymeersch, “SISO RIS-enabled joint 3D downlink localization and synchronization,” in _IEEE International Conference on Communications_ , 2021.
  * [9] Z. Abu-Shaban, K. Keykhosravi, M. F. Keskin, G. C. Alexandropoulos, G. Seco-Granados, and H. Wymeersch, “Near-field localization with a reconfigurable intelligent surface acting as lens,” in _ICC 2021-IEEE International Conference on Communications_ , 2021.
  * [10] S. Huang, B. Wang, Y. Zhao, and M. Luan, “Near-field RSS-based localization algorithms using reconfigurable intelligent surface,” _IEEE Sensors Journal_ , 2022.
  * [11] J. He, H. Wymeersch, T. Sanguanpuak, O. Silvén, and M. Juntti, “Adaptive beamforming design for mmwave RIS-aided joint localization and communication,” in _2020 IEEE Wireless Communications and Networking Conference Workshops (WCNCW)_ , 2020.
  * [12] X. Hu, C. Zhong, Y. Zhang, X. Chen, and Z. Zhang, “Location information aided multiple intelligent reflecting surface systems,” _IEEE Transactions on Communications_ , vol. 68, no. 12, pp. 7948–7962, 2020.
  * [13] K. Keykhosravi, M. F. Keskin, S. Dwivedi, G. Seco-Granados, and H. Wymeersch, “Semi-passive 3D positioning of multiple RIS-enabled users,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 10, pp. 11 073–11 077, 2021\.
  * [14] J. Sherman, “Properties of focused apertures in the fresnel region,” _IRE Transactions on Antennas and Propagation_ , vol. 10, no. 4, pp. 399–408, 1962.
  * [15] W. Tang, M. Z. Chen, X. Chen, J. Y. Dai, Y. Han, M. Di Renzo, Y. Zeng, S. Jin, Q. Cheng, and T. J. Cui, “Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 421–439, 2020.
  * [16] B. Friedlander, “Localization of signals in the near-field of an antenna array,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 15, pp. 3885–3893, 2019.
  * [17] S. M. Kay, _Fundamentals of statistical signal processing: estimation theory_. Prentice-Hall, Inc., 1993.