Leandro R Lustosa

Ilya Kolmanovsky

Carlos E S Cesnik

Fabio Vetrano

Clarence L " Kelly

" Johnson

Aided Inertial Estimation of Wing Shape

https://doi.org/10.2514/1.G005368 Advanced large-wing-span aircraft result in more structural flexibility and the potential for instability or poor handling qualities. These shortcomings call for stability augmentation systems that entail active structural control. Consequently, the in-flight estimation of wing shape is beneficial for the control of very flexible aircraft. This paper proposes a new methodology for estimating flexible structural states based on extended Kalman filtering by exploiting ideas employed in aided inertial navigation systems. High-bandwidth-rate gyro angular velocities at different wing stations are integrated to provide a short-term standalone inertial shape estimation solution, and additional low-bandwidth aiding sensors are then employed to bound diverging estimation errors. The proposed filter implementation does not require a flight dynamics model of the aircraft, facilitates the often tedious Kalman filtering tuning process, and allows for accurate estimation under large and nonlinear wing deflections. To illustrate the approach, the technique is verified by means of simulations using sighting devices as aiding sensors, and an observability study is conducted. In contrast to previous work in the literature based on stereo vision, a sensor configuration that provides fully observable state estimation is found using only one camera and multiple rate gyros for Kalman filtering update and prediction phases, respectively. = jth rate gyro output (z axis) r ∈ R 3 = displacement vector r f ∈ R 3 = displacement vector in body frame s ∈ R = wing arc length s 0 ∈ R = arbitrary wing arc length s j ∈ R = arc length of jth rate gyro t ∈ R = time instant w g j ∈ R 3 = jth rate gyro noise x EKF ∈ R 3N3M = extended Kalman filter state xs ∈ R 3 = local wing x axis xb ∈ R 3 = body x axis ŷs ∈ R 3 = local wing y axis ŷb ∈ R 3 = body y axis ẑs ∈ R 3 = local wing z axis ẑb ∈ R 3 = body z axis Γ ∈ R = wing dihedral angle δθ ∈ R N = twist modal amplitudes error δθ i ∈ R = ith twist modal amplitude error δϕ ∈ R N = anhedral modal amplitudes error δϕ i ∈ R = ith anhedral modal amplitude error δψ ∈ R N = sweep modal amplitudes error

δψ i ∈ R = ith sweep modal amplitude error ε j ∈ R 3 = jth rate gyro drift Θ ∈ R M×N = modal to Euler pitch rates transformation θ ∈ R = local twist angle θ ∈ R N = twist modal amplitudes θ i ∈ R = ith twist modal amplitude θi ∈ R = ith twist mode shape Λ ∈ R = wing sweep angle Φ ∈ R M×N = modal to Euler roll rates transformation ϕ ∈ R = local anhedral angle ϕ ∈ R N = anhedral modal amplitudes φ ∈ R 3N
= angular to modal rates transformation φi ∈ R = ith anhedral mode shape

ϕ i ∈ R = ith anhedral modal amplitude Ψ ∈ R M×N = modal to Euler yaw rates transformation ψ ∈ R = local sweep angle ψ ∈ R N = sweep modal amplitudes ψ i ∈ R = ith sweep modal amplitude ψi ∈ R = ith sweep mode shape ω b ∈ R 3 = rigid-body ground truth velocity ωb ∈ R 3 = rigid-body rate gyro output ω j ∈ R 3
= jth station ground truth velocity ωj ∈ R 3 = jth rate gyro output

I. Introduction

T HE resulting increased structural flexibility in high-aspect-ratio wing designs progressively couples aeroelastics with flight dynamics and generally gives rise to instability and/or poor handling qualities. Stable and safe piloting calls for stability augmentation systems (SASs) that simultaneously account for elastic and traditional flight dynamics degrees of freedom. Consequently, in-flight estimation of wing shape is beneficial for the control of very flexible aircraft (VFA). This paper presents a novel technique for wing shape estimation in the presence of large wing deflections based on approaches from the field of inertial navigation.

A state estimator is herein called a flight dynamics model-based estimator if its design requires a mathematical dynamic model of the aircraft relating actuator values (e.g., ailerons deflection, flap configuration) to system state trajectory (e.g., velocity, attitude, structural modal amplitudes). Such estimators are vehicle dependent, and their application to modified aircraft designs entails lengthy reidentification/tuning campaigns, e.g., wind tunnel testing, numerical simulations, and flight testing. Additionally, even when a model is available, it is often not precise or observable enough to support state observer-based controller design. To address this challenge, this paper proposes a model-free approach to wing shape estimation using inertial sensors. The technique avoids strain gauges and fiber Bragg gratings (FBG) [START_REF] Ma | Fiber Bragg Gratings Sensors for Aircraft Wing Shape Measurement: Recent Applications and Technical Analysis[END_REF] sensing due to their structural straindisplacement modeling requirements, although their use as complementary aiding sensors is still feasible.

An alternative flight dynamics model-free approach is to employ sighting devices and stereo-vision tracking of visual references on the wing. Whereas by using existing rivets or joint lines (natural candidates for markers) Tagai et al. [START_REF] Tagai | Wing Deformation Measurement as Bases for In-Flight Aerodynamics[END_REF] failed to achieve satisfactory estimation accuracy on a Beechcraft Type 65 Queen Air, experiments with artificial markers yielded subpixel precision. Alternatively, Meyer et al. [START_REF] Meyer | Optical In-Flight Wing Deformation Measurements with the Image Pattern Correlation Technique[END_REF] employed the image pattern correlation technique (IPCT) to exploit surfaces covered with random patterns (i.e., texture) to estimate wing, slat, and flap deformations on an Airbus A320. Rudnik and Schwetzler [START_REF] Rudnik | High Lift INflight VAlidation (HINVA)-Overview About the 1st Flight Test Campaign[END_REF] exploited the same experimental apparatus to study low-speed high lift and significantly enhance the accuracy and reliability of the prediction of maximum lift for commercial aircraft.

Similarly, Kirmse [START_REF] Kirmse | Recalibration of a Stereoscopic Camera System for In-Flight Wing Deformation Measurements[END_REF] employed IPCT to estimate the wing shape of a VUT100 Cobra airplane. During the flight, aircraft vibration caused significant camera misalignment that ground vibration tests failed to reproduce. To correct for misalignment during postflight analysis, Kirmse [START_REF] Kirmse | Recalibration of a Stereoscopic Camera System for In-Flight Wing Deformation Measurements[END_REF] proposed a technique for recalibration of extrinsic parameters (e.g., camera position and orientation) based on visual data. Camera misalignment also posed a challenge in Ref. [START_REF] Sanches | Designing a Vision-Based System to Measure In-Flight Wing Deformation[END_REF] in the study of a high-altitude long-endurance vehicle. In the latter, intrinsic parameters remained practically constant during flight.

In computer vision, outdoor lighting poses additional problems. In fact, Kurita et al. [START_REF] Kurita | In-Flight Wing Deformation Measurement[END_REF] devised the flight plan for experiments to optimize lighting conditions to study static load deformations on a Cessna Citation Sovereign. Interestingly, because the cameras were installed inside the cabin, the intrinsic optical parameters were a function of cabin pressure due to window displacement (and therefore compound lens reconfiguration) caused by fuselage deformation. Pang et al. [START_REF] Pang | In-Flight Wing Deformation Measurement System for Small Unmanned Aerial Vehicles[END_REF] additionally warn of difficulties due to optimal shutter speed values being dependent on lighting conditions, and light reflection from the wing, which was significantly reduced by application of a light dusting of nonreflective flat white paint for accurate optical tracking of high-intensity light-emitting diodes (LEDs). A similar approach using LEDs and videogrammetry was previously developed for the measurement of wing twist and deflection as a function of time of an F/A-18 research aircraft at NASA's Dryden Flight Research Center [START_REF] Burner | In-Flight Aeroelastic Measurement Technique Development[END_REF]. However, the LED housings, especially outboard toward the trailing edge of the wing, were barely distinguishable from the background and called for manual selection of target image locations. Graves and Burner [START_REF] Graves | Development of an Intelligent Videogrammetric Wind Tunnel Measurement System[END_REF] used retro-reflective tape targets to ensure high-contrast imagery. This allowed for automated target detection employing blob analysis. Liu et al. [START_REF] Liu | Photogrammetric Techniques for Aerospace Applications[END_REF] provide a comprehensive and systematic summary of photogrammetric techniques for aerospace applications.

All computer vision algorithms mentioned above are run offline after flight testing. A different approach in Ref. [START_REF] Deangelis | In-Flight Deflection Measurement of the HiMAT Aeroelastically Tailored Wing[END_REF] pursued in-flight estimation from active LED markers detection by analog electro-optic receivers and allowed for continuous data sampling and processing. However, accurate calibration of analog systems is challenging and prone to noise; thus, engineering practice usually favors digital systems. Nevertheless, after an extensive and complex calibration procedure, this system successfully provided real-time continuous-time data for a 0.44-scale model of a 17,000-lb fighter airplane.

In-flight wing shape estimators for use in SASs must satisfy hard real-time computation constraints and high-bandwidth requirements, thus prohibiting employment of standalone computer vision systems. A practical alternative often pursued in navigation systems is to resort to inertial sensors for high-bandwidth data while exploiting lowbandwidth aiding sensors (e.g., global navigation satellite system receivers, cameras) for bounding divergent estimation errors. In Refs. [START_REF] Kotikalpudi | Real-Time Shape Estimation for a Small Flexible Flying-Wing Aircraft[END_REF][START_REF] Pang | Strain State Estimation of Very Flexible Unmanned Aerial Vehicle[END_REF], extended Kalman filters (EKFs) were implemented, which combined data from flight dynamics models and inertial sensors to obtain accurate results for bending; unfortunately, in Ref. [START_REF] Kotikalpudi | Real-Time Shape Estimation for a Small Flexible Flying-Wing Aircraft[END_REF] the estimator failed to characterize twist satisfactorily. Additionally, the techniques in Refs. [START_REF] Kotikalpudi | Real-Time Shape Estimation for a Small Flexible Flying-Wing Aircraft[END_REF][START_REF] Pang | Strain State Estimation of Very Flexible Unmanned Aerial Vehicle[END_REF] rely on an EKF formulation that precludes simple tuning due to lack of information on the statistics of the flight dynamics modeling errors. Finally, the approach in Ref. [START_REF] Kotikalpudi | Real-Time Shape Estimation for a Small Flexible Flying-Wing Aircraft[END_REF] assumed small wing deflections, thus hindering its application to VFA.

To the best of the authors' knowledge, this paper is the first attempt to use the aided inertial navigation framework to estimate wing shape without resorting to flight dynamics models, and its main contribution is a novel extended Kalman filter-based technique for wing shape estimation based on distributed rate gyros and optional auxiliary sensors. Because of the use of Kalman filtering, there is great flexibility in the choice of aiding sensors type, precision, sampling rates, and number. Herein a case study with sighting devices illustrates the technique and its properties. The resulting real-time estimator 1) does not require a flight dynamics model of the aircraft, 2) has a straightforward tuning procedure, 3) has an adequate bandwidth for VFA control purposes, and 4) allows for large and nonlinear wing deflections.

The remainder of the paper is organized as follows: Sec. II proposes a standalone rate gyro algorithm for wing shape estimation, and Sec. III illustrates the proposed formulation of an extended Kalman filter to incorporate additional sensor data and enhance estimation precision. Section IV illustrates the approach using cameras as aiding devices. Section V investigates the sensor placement impact on filter observability, and Sec. VI summarizes the main findings of the paper.

II. Standalone Rate Gyro Wing Shape Estimation

A. Algorithm Derivation

This section introduces a technique for estimating wing shape using an array of rate gyros. Consider an aircraft equipped with M rate gyros distributed along its right wing (see Fig. 1). Consider additionally, the aircraft inertial navigation system (INS) and, accordingly, its associated global aircraft frame f xb ; ŷb ; ẑb g. The global aircraft frame will be considered as the rigid-body reference to which flexible frames will be defined. The structural deflection is modeled by defining a deformation reference line (DRL) fixed to the wing structure (see Fig. 1). Furthermore, the DRL is defined such that it contains all M rate gyros. The DRL shape is parameterized by Euler angle functions ϕs; t, θs; t, and ψs; t, where s ∈ R denotes the arc length along the reference line (s 0 at the wing root). In this work, the angles ϕs; t, θs; t, and ψs; t are called local anhedral, local twist, and local sweep angles, respectively, and they model the orientation of the local wing frame f xs; ŷs; ẑsg with respect to the body-fixed global aircraft frame f xb ; ŷb ; ẑb g. This formulation is similar to the Frenet-Serret trajectory-based frames approach [START_REF] Kreyszig | Differential Geometry[END_REF] and neglects extension effects. Therefore, to increase method accuracy, the DRL should be located as close as possible to the wing neutral axis, i.e., the axis along which there are no longitudinal stresses or strains.

Figure 2 further illustrates the proposed deformation description and axes conventions. If all three local angles are zero at s s 0 , i.e., ϕs 0 ; t θs 0 ; t ψs 0 ; t 0, then the wing section axes at s s 0 align with the aircraft body axes, as Fig. 2b suggests. Therefore, zero local angles do not necessarily correspond to undeformed wing shape. Additionally, spatially constant angle functions model sweep and dihedral angles (Λ and Γ, respectively), as Fig. 2a suggests. This paper defines the DRL as always tangent to the local ŷs. Figures 2b and 2c illustrate orientation and right-handed sign conventions where the ẑ (ψ)-ŷ (θ)-x (ϕ) rotation order is assumed from body to the local wing section. This rotation order introduces a singularity in θ 90°, which is unlikely to happen even in VFA. On the other hand, an alternative rotation order such as ẑ (ψ)-x (ϕ)-ŷ (θ) introduces a singularity in ϕ 90°, which may occur in VFA.

The jth rate gyro is installed at s s j , and it is considered aligned with f xs j ; ŷs j ; ẑs j g without loss of generality. The corresponding angle functions are described as a superposition of N basis modes, i.e., ϕs j ; t X N i1 φi s j ϕ i t for j 1; : : : ; M (1)

θs j ; t X N i1
θi s j θ i t for j 1; : : : ; M (2)

ψs j ; t X N i1
ψi s j ψ i t for j 1; : : : ; M

where the symbols ϕ, φi , and ϕ i denote deformation angle, mode shapes, and modal amplitudes, respectively. Furthermore, a column vector of modal amplitudes is denoted by a bold symbol, e.g., ϕ ϕ 1 ; : : : ; ϕ N T . The basis modes f φi ; θi ; ψi : i 1; : : : ; Ng; are chosen to minimize N by accounting for most recurrent shapes occurring during flight, for instance, by proper orthogonal decomposition of finite element method simulations to varying loads or eigenmodes of a representative mass condition [START_REF] Medeiros | Computational Aeroelasticity Using Modal-Based Structural Nonlinear Analysis[END_REF].

The jth rate gyro yields angular velocity measurements such that 0

B B @ Pj Qj Rj 1 C C A 0 B @ P j Q j R j 1 C A ε j w g j ( 4
)
where ε j and w g j are, respectively, jth rate gyro's drift and noise. Additionally, ω j P j ; Q j ; R j T and ωj Pj ; Qj ; Rj T denote ground truth and measured angular velocities described in the jth rate gyro local frame f xs; ŷs; ẑsg, respectively. Similarly,

ω b P b ; Q b ; R b T and
P j t -D b j tP b t Q j t -D b j tQ b t R j t -D b j tR b t 1 C C C A (5) a) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) b) ( , )
( , ) = 0 where D b j t is the direction cosine matrix from f xb t; ŷb t; ẑb tg to f xs j ; t; ŷs j ; t; ẑs j ; tg, given by

( , ) = 0 ( , ) ( , ) ( , ) c) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , )
D b j t 2 4
cos θ cos ψ cos θ sin ψ sin θ cos ϕ sin ψ sin ϕ sin θ cos ψ cos ϕ cos ψ sin ϕ sin θ sin ψ sin ϕ cos θ sin ϕ sin ψ cos ϕ sin θ cos ψsin ϕ cos ψ cos ϕ sin θ sin ψ cos ϕ cos θ

To avoid clutter, the compound function notation H ϕ s j ; t, H θ s j ; t, and H ψ s j ; t will be used for the left-hand sides of Eqs. [START_REF] Kurita | In-Flight Wing Deformation Measurement[END_REF][START_REF] Pang | In-Flight Wing Deformation Measurement System for Small Unmanned Aerial Vehicles[END_REF][START_REF] Burner | In-Flight Aeroelastic Measurement Technique Development[END_REF]. Furthermore, note that

∂ ∂t ϕs j ; t X N i1
φi s ϕ i t for j 1; : : : ; M

and thus, using Eq. ( 5) and the definition of ω j , one obtains

H T ϕ ω j -D b j ω b X N i1 φi s j _ ϕ i t
for j 1; : : : ; M (13)

H T θ ω j -D b j ω b X N i1 θi s j _ θ i t for j 1; : : : ; M (14) H T ψ ω j -D b j ω b X N i1 ψi s j _ ψ i t for j 1; : : : ; M (15)
which can be recast as the following systems of linear algebraic equations:

0 B B @ H T ϕ s 1 ; tω 1 -D b j ω b . . . H T ϕ s M ; tω M -D b j ω b 1 C C A 2 6 6 4
φ1 s 1 • • • φN s 1 . . . . . . . . . φ1 s M • • • φN s M 3 7 7 5 |{z} Φ 0 B B @ _ ϕ 1 . . . _ ϕ N 1 C C A (16) 0 B B @ H T θ s 1 ; tω 1 -D b j ω b . . . H T θ s M ; tω M -D b j ω b 1 C C A 2 6 6 4
θ1 s 1 • • • θN s 1 . . . . . . . . . θ1 s M • • • θN s M 3 7 7 5 | {z } Θ 0 B B @ _ θ 1 . . . _ θ N 1 C C A (17) 0 B B @ H T ψ s 1 ; tω 1 -D b j ω b . . . H T ψ s M ; tω M -D b j ω b 1 C C A 2 6 6 4
ψ1 s 1 • • • ψN s 1 . . . . . . . . . ψ1 s M • • • ψN s M 3 7 7 5 | {z } Ψ 0 B B @ _ ψ 1 . . . _ ψ N 1 C C A (18)
Note that Φ, Θ, and Ψ depend only on the choice of basis functions and sensor placement. Assuming M > N and invertibility of Φ T Φ, Θ T Θ, and Ψ T Ψ, the least-squares solutions of Eqs. [START_REF] Medeiros | Computational Aeroelasticity Using Modal-Based Structural Nonlinear Analysis[END_REF][START_REF] Crassidis | Optimal Estimation of Dynamic Systems[END_REF][START_REF] Krener | Measures of Unobservability[END_REF] are given by

0 B B @ _ ϕ 1 . . . _ ϕ N 1 C C A Φ T Φ -1 Φ T 0 B B @ H T ϕ s 1 ; tω 1 -D b j ω b . . . H T ϕ s M ; tω M -D b j ω b 1 C C A (19) 0 B B @ _ θ 1 . . . _ θ N 1 C C A Θ T Θ -1 Θ T 0 B B @ H T θ s 1 ; tω 1 -D b j ω b . . . H T θ s M ; tω M -D b j ω b 1 C C A (20) 0 B B @ _ ψ 1 . . . _ ψ N 1 C C A Ψ T Ψ -1 Ψ T 0 B B @ H T ψ s ; tω 1 -D b j ω b . . . H T ψ s M ; tω M -D b j ω b 1 C C A (21)
Integrating Eqs. (19-21) in time yields the standalone rate gyro wing shape estimation (RG-WSE) algorithm. All estimated outputs from the RG-WSE algorithm are represented in this work with a hat on its corresponding symbol. For instance, φ1 t is the RG-WSE estimate of ϕ 1 t. Additionally, in state-space representation, Eqs. (19-21) can be recast as _

x φx; u, where x ϕ; θ; ψ and u ω 1 ; : : : ; ω M ; ω b . Figure 3 illustrates the algorithm's input and output variables. Inputs are angular velocity measurements at M multiple wing locations and at the global aircraft frame, and outputs are modal amplitudes in a truncated N-dimensional space of deformation modes, i.e., φ, θ, and ψ ∈ R N .

B. RG-WSE Angles to Displacement Formulation Transformation

Although the RG-WSE algorithm uses an angle-based parameterization of wing deformation, displacement-based coordinates are often desired instead. This section presents how to transform between these descriptions.

Figure 4 illustrates both descriptions at the same point s s j . At any given time t, the displacement formulation rs; t maps a given point s units of distance away from the wing root along the wing shape. At a fixed instant of time where r f is the displacement vector described in the body frame, and Dϕ; θ; ψ is the direction cosine matrix given by Dϕ; θ; ψ 2 4 cos θ cos ψ cos θ sin ψ sin θ cos ϕ sin ψ sin ϕ sin θ cos ψ cos ϕ cos ψ sin ϕ sin θ sin ψ sin ϕ cos θ sin ϕ sin ψ cos ϕ sin θ cos ψsin ϕ cos ψ cos ϕ sin θ sin ψ cos ϕ cos θ 3 5

(24)

where ϕs; t, θs; t, and ψs; t are obtained according to Sec. II.A.

C. RG-WSE-Based Wing Coordinates

In previous sections, the wing deformation was kinetically described as an elastica with torsion and two bending degrees of freedom related to the local angles attached to points on the wing. All rate gyros were considered installed on the DRL. However, additional sensors might be installed in stations disjoint from the DRL. In this case, an extrapolation model for positioning points not belonging to the deformation line is required.

Each point on the wing can be described by the sum of a reference point r f s 0 ; t in the deformation line in s 0 , and a Δp vector orthogonal to the DRL. This work assumes that cross sections normal to the DRL are rigid under wing deformation. This implies that Δp is always written as Δx; 0; Δz in the f xs 0 ; ŷs 0 ; ẑs 0 g basis, where Δx and Δz are constant in time. Therefore, every point on the wing at any given time can have its position described as the triple s 0 ; Δx; Δz; these are hereafter called wing coordinates, and relate to global frame position p f according to p f r f s 0 ; t D ϕs 0 ; t; θs 0 ; t; ψs 0 ; t T 0

@ Δx 0 Δz 1 A (25)

III. RG-WSE Error Model for Extended Kalman Filtering

For fusing sensor output data with extended Kalman filtering, this work uses a linearized model for RG-WSE errors. To best conform with RG-WSE, the state vector to be estimated is composed of modal amplitude errors and rate gyro drifts. Accordingly,

x EKF δϕ 1 ••• δϕ N δθ 1 ••• δθ N δψ 1 ••• δψ N ε T 1 ••• ε T M ε T b T ( 26
)
where δϕ i , δθ i , and δψ i are defined as the difference between the ground truth and the RG-WSE-computed values. Therefore,

δϕ i ϕ i -
0 0 0 0 # | {z } A x EKF ∂φ ∂u |{z} G w g (28)
where w g aggregates all the sensor noise vectors, i.e., w g w g 1 ; : : : ; w g M ; w g b , and all the Jacobians are evaluated at x; û, obtained here numerically by means of first-order finite differences. Similarly, a discrete-time model is obtained by

x EKF t k ≈ x EKF t k-1 _ x EKF Δt AΔt I | {z } F k x EKF t k-1 G EKF |{z} G k ∂w g (29)
where ∂w g is the discrete-time white noise equivalent of w g . The statistics of ∂w g should conform to rate gyro manufacturer specifications. The prediction step of the EKF is then written as

xkjk-1 F k xk-1jk-1 P kjk-1 F k P k-1jk-1 F T k G k Q k G T k ( 30
)
where Q k , P kjk , and P kjk-1 are the covariance of the rate gyros noise, the a posteriori error covariance matrix, and the a priori error covariance, respectively. This formulation provides a simple recipe for EKF tuning because Q k is directly related to rate gyro manufacturer specifications. However, artificial process noise inflation might be necessary to account for linearization errors and aliasing in RG-WSE. Figure 5 illustrates the proposed aided inertial wing shape estimator architecture. Additional aiding sensors (e.g., cameras, strain gauges, accelerometers) could be added to the architecture by means of adequate modeling of their EKF observation equations. The next section illustrates an application using sighting devices as aiding sensors. Periodical correction of RG-WSE must be performed with the acquired EKF estimates to decrease the effect of linearization errors. Additionally, the estimated rate gyro drifts are calibrated during flight by subtracting in software their biases based on Eq. ( 4).

IV. Camera-Aided Rate Gyro Wing Shape Estimation

It will now be shown that integration in time of rate gyro drift and noise in Eqs. (19-21) yields unbounded RG-WSE estimation errors. Therefore, additional sensors are required to prevent estimation errors from diverging in time. Auxiliary devices can be incorporated in the RG-WSE framework by means of EKF observation equations. This work proposes an implementation of EKF update phase that employs mappings from modal amplitude errors (i.e., δϕ, δθ, and δψ) to auxiliary sensors output. This section illustrates this for a sighting device case study where visual markers are 1) placed at known positions on the wing, and are 2) tracked by cameras rigidly installed close to the wing root. For a given camera, the auxiliary sensor output is defined as the difference in pixels between a measured marker position and the predicted marker position based on RG-WSE outputs φ, θ, and ψ, and camera-marker positioning. Multiple markers, M k , k 1; : : : ; K, tracked by multiple cameras, C l , l 1; : : : ; L, are considered.

A. EKF Observation Equations

The proposed rate gyro and sighting device data fusion technique is based on tracking visual markers M k , each with a priori known wing location s k ; x k ; z k . Recall from Sec. II.C that points fixed to the wing are modeled assuming constant wing coordinates during deformation. For each marker M k , its position with respect to a camera C l , denoted by p k∕l ∈ R 3 , is described in the C l camera coordinate frame using a pinhole camera model as

p k∕l -D f l r C l f D f l r f s k ; t D f l Dϕs k ; t; θs k ; t; ψs k ; t T 0 B B @ x k 0 z k 1 C C A (31)
where D f l and r C l f denote the direction cosine matrix from body frame to camera C l frame, and camera position with respect to wing root in body frame, respectively (see Fig. 6). All cameras are assumed fixed in body frame and thus D f k is constant. Hence the camera measurement z k∕l ∈ R 2 in pixels is given by

z k∕l Π fp k∕l 1 0 0p k∕l (32)
where f ∈ R is the focal length in pixels and Π is defined by

Π 0 1 0 0 0 1 (33)
Given these assumptions, z k∕l depends exclusively on modal amplitudes ϕ i , θ i , and ψ i , i 1; : : : ; N. Therefore, z k∕l linearization yields

z k∕l ϕ; θ; ψ -z k∕l φ; θ; ψ | {z } ≜Δz k∕l ∂z ∂ϕ k∕l δϕ ∂z ∂θ k∕l δθ ∂z ∂ψ k∕l δψ (34)
Consequently, an appropriate model for an observation equation is given by

Δz k∕l ∂z ∂ϕ k∕l ∂z ∂θ k∕l ∂z ∂ψ k∕l ∂z ∂ε 1 k∕l • • • ∂z ∂ε M k∕l ∂z ∂ε b k∕l | {z } H k∕l x EKF w k∕l CAM ( 35
)
where w k∕l CAM ∈ R 2 is additive white Gaussian noise with statistics depending on camera quality, tracking algorithm performance, and marker-camera positioning, ∂z k∕l ∕∂ε j 0 for all j 1; : : : ; M, and ∂z k∕l ∕∂ε b 0. With a linear sensor model in hand, the EKF incorporates measurements through the classical Kalman filter formulas [START_REF] Crassidis | Optimal Estimation of Dynamic Systems[END_REF]:

y k Δz k∕l -H k∕l xkjk-1 S k H k∕l P kjk-1 H k∕l T R k K k P kjk-1 H k∕l T S -1 k xkjk xkjk-1 K k y k P k k I -K k H k∕l P kjk-1 (36)
where R k Ew k∕l CAM w k∕l CAM T is the measurement noise covariance, and the mean of the measurement noise is assumed to be zero. In practice, artificial R k inflation with respect to w k∕l CAM might be necessary to account for additional sources of uncertainties, e.g., camera vibration.

Rate gyro-based EKF prediction updates occur at a higher rate than EKF observation updates from the camera measurements. Thus

Camera axis definition and symbols b) Camera view

Fig. 6 Camera-related symbols and frame definitions. between two updates from the camera measurements per Eq. (36) there will be multiple prediction updates from the gyro measurements per Eq. ( 30). This allows for heterogeneous sensor configurations with respect to sampling time, and favors nonuniform sampling. In this work, camera tracking algorithms are not considered; instead, it is assumed that image processing already took place with precision represented by w k∕l CAM . Examples of target-tracking image processing algorithms in the context of wing deformation are given in Refs. [START_REF] Pang | In-Flight Wing Deformation Measurement System for Small Unmanned Aerial Vehicles[END_REF][START_REF] Graves | Development of an Intelligent Videogrammetric Wind Tunnel Measurement System[END_REF]. However, Pang et al. [START_REF] Pang | In-Flight Wing Deformation Measurement System for Small Unmanned Aerial Vehicles[END_REF] build upon a stereo vision setup to update an EKF with estimated 3D positions of markers. The present camera-aided rate gyro wing shape estimation (CRG-WSE) approach eliminates the need for stereo vision because it uses the 2D projections of markers in the camera plane directly. As a consequence, any number of cameras can be employed independently, rendering camera frame synchronization unimportant.

B. Standalone and Aided Methods Comparison

By means of computer simulation, two experiments were conducted to illustrate the degrees of observability of standalone and aided observer algorithms. In both numerical experiments, right wing shape estimation of an aircraft is pursued while performing a harmonic motion given by ϕs; t Γ sΓ 0 sinωt θs; t 0

ψs; t Λ ω b t 0 (37)
where 0 < t < t sim , 0 < s < b. Additionally, Γ, Γ 0 , Λ, and Λ 0 are scalar constants given in Table 1. The overall trajectory of the wing due to this deformation profile is illustrated in Fig. 7a. For algorithm implementation, polynomial shape basis and N 2 are assumed, thus matching the number of modes in ground truth, and, accordingly,

φ1 s θ1 s ψ1 s 1 φ2 s θ2 s ψ2 s s (38)
leading to ground truth values given by ϕt ϕ 1 t

ϕ 2 t ! Γ Γ 0 sinωt ! ; θt θ 1 t θ 2 t ! 0 0 ; ψt ψ 1 t ψ 2 t ! Λ 0 (39)
Ten rate gyros were installed in the wing according to the positioning defined in Fig. 7b. In the camera-aided case, additional visual trackers and sighting devices are installed as also depicted in Fig. 7b. Their biases and noise characteristics are described in Table 2, and are purposely deteriorated with respect to the performance of low-cost inertial measurement units (IMUs) readily available in the market (e.g., the Xsens MTi 100-series) to reduce simulation time.

Figure 8 illustrates the simulation results for both standalone and camera-aided experiments. Modal amplitude and bias estimation errors are plotted with their respective 2σ EKF covariances. Moreover, δϕ p are plotted separately due to their distinct orders of magnitude.

Note that standalone estimation has its errors diverging with time; thus, its use imposes an upper bound on the duration of missions. On the other hand, the camera-aided estimator yields bounded errors. However, for the given sensor imperfection statistics, rate gyro drift is negligibly estimated because its initial covariance is only marginally greater than the steady-state filter covariance.

V. Observability Analysis

This section studies the impact of rate gyros and visual markers number and placement on the proposed camera-aided estimator precision. For that purpose, the numerical experiment in Sec. IV.B is repeated with different sensor configurations. Since N 2 was assumed, a minimum of two rate gyros is required for RG-WSE implementation, as discussed in Sec. II.A. A varying quantity from 2 to 10 rate gyros in even numbers is examined assuming equally spaced placement along the semispan b (see Fig. 9a). For each rate gyro configuration, the visual markers quantity is also varied from 2 up to 10 in even numbers. Their placement starts at the wingtip (s b) and proceeds inward with 2 m spacing (see Fig. 9b). Additionally, odd-and even-numbered markers are placed toward the leading and trailing edges, respectively, but always 2 m away from the DRL.

Although a universal figure of merit for measuring observability is nonexistent [START_REF] Krener | Measures of Unobservability[END_REF], observability is evaluated herein through end-ofsimulation EKF P kjk covariances comparison. In this sense, smaller As expected, observability increases with the number of sensors and markers. However, as the number of markers increases, the impact on system accuracy of increasing the number of rate gyros gets diminished. Additionally, modal error standard deviations decrease log-linearly with the number of visual markers.

The comments above suggest some design guidelines for choosing sensor quantity numbers. First, the designer determines N by deciding how many (and which) basis modes best represent the bulk of the in-flight wing deformations. The number of rate gyros M should match N, and there is little gain in adding additional units. Subsequently, the number of markers K is chosen to attain the desired performance or should be as large as the visual tracking algorithm allows without violating robustness and real-time requirements. Note that a system with no camera updates would have estimation errors diverging, whereas a system with no rate gyros would have low bandwidth as measurements are available only at camera sampling instants. Although RG-WSE imposes rate gyros installation on the DRL, visual markers are allowed to be anywhere on the wing. However, some practical considerations must be respected. For example, Fig. 11 shows the outcome of a numerical experiment where five equally spaced visual markers are positioned on the reference deformation line (instead of closer to the leading and trailing edges as before). The noticeably degraded observability in θ, by a factor of 10 in precision with respect to the other components (and previous experiments in Sec. IV.B), reinforces the intuitive result that a line cannot fully define a rigid-body orientation. Therefore, rectilinear placement of markers should be avoided if wing twist estimation is pursued.

VI. Conclusions

This paper proposed a novel technique for high-bandwidth realtime estimation of wing shape using an array of rate gyros, and lays the groundwork for tight EKF interfacing with additional sensors. Tight EKF coupling eliminates suboptimal assumptions that often reduce precision and make tuning impracticable in reality. To further facilitate tuning, the proposed technique avoids resorting to flight dynamics models because their uncertainties are difficult to characterize. Furthermore, design guidelines are provided for tuning and deciding on sensor number and placement.

The proposed technique inherits the benefits of Kalman filtering, e.g., multirate asynchronous aiding sensors integration, low computational cost, and fault detection, isolation, and recovery capabilities. In particular, this work shows that tight optical system integration yields a fully observable state estimation using a single camera. This contrasts with previous research, where loosely coupled integration called for stereo vision approaches (thus two sighting devices or more). Finally, much room remains for future research. This includes the analysis of the case when gyros are not installed on the DRL, and impact on the accuracy at higher frequencies when the camera may fail to capture all the information.

NomenclatureD

  ∈ R 3×3 = direction cosine matrix H ∈ R 3×3 = angular velocity to Euler angles rate transformation H θ ∈ R 1×3 = angular velocity to pitch rate transformation H ϕ ∈ R 1×3 = angular velocity to roll rate transformation H ψ ∈ R 1×3 = angular velocity to yaw rate transformation K ∈ Z = number of cameras M ∈ Z = number of rate gyros N ∈ Z = number of mode shapes P b ∈ R = rigid-body angular velocity (x axis) Pb ∈ R = rigid-body rate gyro output (x axis) P j ∈ R = jth station angular velocity (x axis) Pj ∈ R = jth rate gyro output (x axis) p k∕l ∈ R 3 = position of marker k with respect to camera l Q b ∈ R = rigid-body angular velocity (y axis) Qb ∈ R = rigid-body rate gyro output (y axis) Q j ∈ R = jth station angular velocity (y axis) Qj ∈ R = jth rate gyro output (y axis) R b ∈ R = rigid-body angular velocity (z axis) Rb ∈ R = rigid-body rate gyro output (z axis) R j ∈ R = jth station angular velocity (z axis) Rj ∈ R

Fig. 2

 2 Fig. 2 Deformation examples.

Fig. 1

 1 Fig. 1 Rate gyros placement and deformation reference line s definition.

Fig. 3

 3 Fig. 3 Schematic of RG-WSE algorithm.

Fig. 4

 4 Fig.4Angle-based and displacement-based descriptions.

Fig. 5

 5 Fig. 5 Aided inertial wing shape estimator overall architecture.

Fig. 7

 7 Fig. 7 Simulation wing trajectory and sensor positioning.

Fig. 8

 8 Fig.8Comparison of estimated error and predicted EKF covariance for standalone estimation (on the left) and aided estimation (on the right).

b)Fig. 9

 9 Fig. 9 Sensor positioning patterns for observability analysis experiments.

Fig. 10 EKF

 10 Fig. 10 EKF 2σ values at the end of the simulation (t 100 s) for different instrumentation quantities.

  ωb Pb ; Qb ; Rb T denote ground truth and measured angular velocities expressed in the global rate gyro frame f xb ; ŷb ; ẑb g, respectively. Angular velocities are related to the derivatives of the Euler angles according to

	∂ ∂t	0 B B @	ϕs j ; t θs j ; t	1 C C A
			ψs j ; t
			2 6 6 6 4 |{z} 3 1 sin ϕs j ; t tan θs j ; t cos ϕs j ; t tan θs j ; t 0 cos ϕs j ; t -sin ϕs j ; t sin ϕs j ;t 7 7 7 5 cos ϕs j ;t 0 cos θs j ;t cos θs j ;t
			0	Hϕs j ;t;θs j ;t;ψs j ;t
		×	B B B @

Table 1

 1

		Wing shape movement description in
			simulation
	Γ	Γ 0	Λ	t sim	b	ω
	-0.09 rad 0.03 rad 0.36 rad 20 s 17 m π rad/s

Table 2

 2

	covariances are associated with more observable, thus more precise,
	configurations. Figure 10 displays end-of-simulation EKF standard
	deviations for each proposed configuration. Because EKF covariances
	vary slightly across simulation runs, due to different EKF Jacobians
	values occurring due to RG-WSE random errors, Fig. 10 displays the
	root mean square (RMS) values of Monte Carlo simulations with 10
	independent runs for each configuration. Each independent simulation
	randomly samples different rate gyros drift ε

Sensor bias and noise statistics in simulation ε b ε j Rate gyro noise density (w g j ) Camera noise standard deviation (ω k∕l CAM ) 0.1 ⋅ 1; 1; 1 deg ∕s 0.2 ⋅ 1; 1; 1 deg ∕s 0.01 deg ∕s∕Hz 0.01 m (assuming f 1) j and noise w g j values.

Acknowledgment

This material is based upon work supported by Airbus in the frame of the Airbus-Michigan Center for Aero-Servo-Elasticity of Very Flexible Aircraft (CASE-VFA).


Joris Costes

Gabriele Facciolo

Rafael Grompone Von Gioi

Josselin Kherroubi

Enric Meinhardt-Llopis

Jean-Michel Morel

A contrario dip picking for borehole imaging

We describe an algorithm to perform automatic dip picking on borehole images. One key element of the proposed method is a statistical validation, based on the a contrario theory, which is used to decide whether each candidate dip is to be accepted or not. The proposed method also uses a randomized Hough transform, which greatly improves the processing speed, allowing for a real-time detection of dips during image visualization. In addition, the same algorithm can be applied at different scales to provide a multi-resolution analysis of the structures. Our experiments show that the proposed algorithm produces reliable dip picking by an evaluation on three manually annotated boreholes: the proposed method detects from 60% to 90% of the dips annotated by an expert, depending on the complexity of the data.

INTRODUCTION

An essential tool in geological exploration are well logs, which record the measurement of physical quantities along the depth of a borehole. Modern measurement devices acquire two-dimensional data over the whole surface of the well. Meaningful information from the resulting borehole images can then be extracted using powerful image processing tools [START_REF] Gonzalez | Integrated multiphysics workflow for automatic rock classification and formation evaluation using multiscale image analysis and conventional well logs: Society of Petrophysicists and Well-Log Analysts[END_REF][START_REF] Jiang | Deeplearning-based vuggy facies identification from borehole images[END_REF].

Dip picking is an important task in the geological interpretation of borehole images for structural, sedimentary, and fracture analysis [START_REF] Lofts | The recognition of artefacts from acoustic and resistivity borehole imaging devices[END_REF]. It consists in extracting the planar structures observable in the borehole image (known as dips) and characterizing them. Borehole images are high-resolution measurements allow-ing for dip characterization [START_REF] Lowell | Borehole imaging: applications and case hisotries[END_REF]. As shown in Figure 1, the imaging measurements are acquired on a cylinder-shaped borehole and result from the unrolling of this cylinder. As a consequence, the horizontal axis is the apparent azimuth angle and the vertical axis is the measured depth along the borehole trajectory. Any planar event crossing the borehole results approximately on a one period sinusoid trace on the unrolled image. As soon as this trace is extracted, the plane characteristics-dip inclination, azimuth and depth-can be directly computed.

The task of dip picking can be performed manually by visual inspection and annotation of the borehole images. This process is very time-consuming and requires a qualified team of experts to annotate a whole well dataset. Thus, several algorithms have been proposed in the literature to help in this task, either by computing the dips automatically, or in a semi-supervised way. Dip picking can be the basis of further geometrical analysis [START_REF] Trice | A methodology for applying a non-unique morphological classification to sine wave events picked from borehole image log data[END_REF][START_REF] Kherroubi | Lamination analysis from electrical borehole images: A quantitative workflow[END_REF].

The Hough transform [START_REF] Hough | Method and means for recognizing com-plex patterns[END_REF][START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF][START_REF] Illingworth | A survey of the hough transform: Computer Vision[END_REF]) is a feature extraction technique widely used to detect structures such as lines, ellipses, or sinusoids. In the Hough transform each pixel cast votes for the features compatible with its value and coordinates, the votes of all the pixels are accumulated in the parameter space of the sought structure. The final state of the accumulator space is the result of the transformation. Then, maxima in this transformed space are used to extract candidate structures.

In the pioneering work of [START_REF] Torres | Method for determining dip and strike angles in borehole ultrasonic scanning tool data[END_REF], the Hough transform is used to extract sinusoidal structures in edge maps [START_REF] Marr | Theory of edge detection[END_REF] extracted from borehole images. A later improvement of this method uses directly the gradient of the borehole images instead of an edge detector [START_REF] Saito | Method of determining dips and azimuths of fractures from borehole images[END_REF]. Other methods also using the Hough transform are [START_REF] Hall | Automatic extraction and characterisation of geological features and textures from borehole images and core photographs[END_REF] and [START_REF] Van Ginkel | Robust curve detection using a radon transform in orientation space[END_REF]. These methods produce an automatic dip picking, their main drawback being the computational Figure 1: Borehole image: from cylinder to an unrolled image. On the right image, a typical colormap used by geologists is applied on the grayscale image. We can clearly see a fault crossing the middle of the image (the sinusoid with high amplitude), surrounded by many sinusoids with very small amplitude, thus nearly horizontal lines, corresponding to the sedimentary structure. Image courtesy of Schlumberger.

cost when a fine quantization of the parameter space is required.

To obtain an interactive user experience, [START_REF] Wang | Method of determining planar events from borehole or core images[END_REF] discloses a user-guided dip picking method, which is based in defining a top and bottom reference curve. These constraints define a smaller subset of the parameter space that can be sampled in real time to find the optimal parameters of a family of sinusoids. [START_REF] Assous | Automated detection of planar geologic features in borehole images[END_REF] exploit the symmetries present in sinusoids to accumulate votes for finding likely structure locations. Sun and Pallottino (2003) use a shortest path algorithm to trace a structure across an entire borehole image, exploiting the connection between the left and right sides of the image. [START_REF] Wedge | Automated structure detection and analysis in televiewer images[END_REF] propose a method for detecting sinusoids that group smaller edge elements, which also exploit the vertical and horizontal sinusoidal symmetries. In addition, the method provides interactive access to the detections.

In this paper, we propose a dip picking method that can process automatically the whole well data in one run, or be integrated in a semi-supervised human/computer interface where the user selects a region of interest. The proposed method involves a statistical validation method for each candidate sinusoid. Its criterion is based on the a contrario theory (Desolneux et al., 2008). This statistical test could be used directly to find a good set of sinusoids (in the sense of number of false alarms, as explained below) by an exhaustive exploration of the whole parameter space. However, to obtain a real-time algorithm, we use the randomized Hough transform [START_REF] Xu | A new curve detection method: Randomized hough transform (rht)[END_REF] to initialize a local search in the parameter space. In addition, the relatively slow variation of dip inclination and azimuth allows for a local exploration of these two parameters, thus reducing even more the required number of random samples to produce reliable candidates. To quantitatively assess the performance of our method, we evaluate its performance on three borehole images with manually annotated dips. We also conduct ablation studies to illustrate the impact of some key components of the proposed method. Our experiments show that the proposed method allows to extract dips reliably on full borehole images. This is supported by the fact that the method is currently used in production for a commercial service. This method could also be applied on any borehole imaging tool (whatever the measurement physics), in wireline logging or measurement while drilling (Drilling and Measurement). The same approach could also be used for extracting dips from slab core pictures: in that case, the Hough transform algorithm must be adapted to extract lines instead of sinusoids.

PROPOSED METHOD

For the sake of simplicity of exposition, we suppose that the borehole is a straight and vertical cylinder (i.e., the borehole does not change its trajectory); thus, the azimuth and depth have an absolute meaning. Furthermore, assuming that the borehole is a perfect cylinder (i.e., that the cross-section is a circle), planar events appear exactly sinusoidal in the image. In practice, boreholes are never straight, but the dip-picking is a local process that only needs the borehole to be locally straight. The method can easily be extended to the general case of a 3D borehole path, including a varying and not perfect circular cross-section; of course, the borehole parameters need to be known.

Let us consider a borehole as a cylinder of radius R (see Figure 2). The borehole image is obtained by unrolling the surface of the cylinder into an image of size W × H. Thus, a pixel (w, h)

∈ [0, W -1] × [0, H -1] corresponds to the 3D point of coordinates R cos 2πw W , R sin 2πw W , 2πR W h . (1)
Note that in the equation above, the third dimension is scaled by a factor 2πR W ; this factor is required to obtain square pixels, having the same height and width in the 3D space.

Let us consider a (non-vertical) 3D plane of equation z = ax + by + c. If we assume that the plane is determined by an inclination angle α and an azimuth angle β, then the parameters correspond to a = K sin(α) cos(β) and b = K sin(α) sin(β), where K is a constant. This plane intersects the cylinder on an ellipse that can be parametrized as θ → (R cos θ, R sin θ, aR cos θ + bR sin θ + c).

(2)

When representing this curve on the image coordinates it becomes

θ → W θ 2π , W 2π a cos θ + b sin θ + c R (3)
or, equivalently, as the graph of the function

h(w) = W 2π a cos 2πw W + b sin 2πw W + c R (4)
in the image plane (w, h). By replacing the expressions of a and b and using a trigonometric addition formula, the expression a cos θ + b sin θ + c can be rewritten as to K sin(α) cos(θ -β) + c, parametrized by the inclination α and the azimuth β. Thus, the function h(w) is a sinusoid of one period, whose phase determines the azimuth, and whose slope determines the inclination of the dip.

The task of dip picking consists in detecting these sinusoidal structures in the borehole images, and assigning an inclination, azimuth and depth to each. Here we propose an algorithm to do that automatically. The core of the proposed method is a statistical test that decides whether any candidate sinusoid is to be accepted or not; in other words, whether the sinusoid separates two groups of pixels whose intensities are significantly different. This criterion is based on the a contrario theory (Desolneux et al., 2008) and is described below.

Such a statistical test could be used directly to find a good set of sinusoids by exhaustive exploration of the whole parameter space. In that case one would traverse the set of all possible sinusoids (a, b, c), and for each sinusoid compute the statistic test which requires evaluating the image at about W positions. Such an exhaustive algorithm would have a computational complexity of N 4 , where N is the number of possible values for each of a, b, c or W . Here, we propose a faster heuristic search on this space, initializing the detection with plausible starting dips.

Algorithm overview

The overall structure of the proposed algorithm is laid out on Figure 3 and on Algorithm 1. The statistical test actually happens on the last steps of the method, because it is used only to validate and refine a rough set of dips obtained by a fast but somewhat imprecise Hough transform.

The algorithm finds all the dips inside a sub-window of the image specified by its vertical location h and scale octave o. The octave determines a vertically down-scaled version (by factors of 2 o ) of the original image. Working on the down-scaled image allows to detect dips with large slopes that would otherwise not fit in a window. But the detection algorithm is the same regardless of the scale. This window can be either selected interactively by the user, or slid along the whole image for an automatic analysis.

The algorithm relies on the local image orientation inside the target window. But two pre-processing are needed: filling-in pad/flap traces, and dequantization and filtering. Then, the local orientations are computed and used twice in the method: first to estimate the main dip inside the window (using the Hough transform), second to find the optimal placement of meaningful dips using an a contrario strategy.

In the following paragraphs we describe each of the steps of the algorithm in detail.

Filling-in pad/flap traces

In the case of pad-based borehole imaging tools, the images are incomplete due to blind spots between the pads/flaps of the sensor (see the vertical stripes on Figure 1). In the proposed method the decision of whether a dip is present or not only uses the available data, ignoring the pad/flap traces. Nevertheless, filling the missing data simplifies the pipeline for proposing candidates by avoiding the tracking of the missing parts. These traces could be inpainted using advanced methods that aim at reconstructing the missing structures, such as the one proposed in [START_REF] Zhang | Generating geologically realistic 3d reservoir facies models using deep learning of sedimentary architecture with generative adversarial networks[END_REF]. In our case it is sufficient to handle the boundary conditions between the known and unknown data. This is done by solving the Laplace equation ∆u = 0 inside the unknown region with boundary conditions given by the known data. This creates a smooth interpolation of the data (see Figure 4) which will not perturb the computation of the local orientation. This is the simplest possible method of data interpolation using PDEs [START_REF] Weickert | Visualization and processing of tensor fields[END_REF]. Note that the quality of this filling is not really important as it is only used to simplify the computations of the intermediate gradients, and that in the final validation the interpolated data is not used.

Dequantization and filtering

The next preprocessing step of the proposed method consists in filtering the image by a Gaussian kernel of size σ. This filtering serves two different goals: the dequantization of the local orientations and the selection of the "grain/resolution" size for the analysis.

The dequantization is essential for computing the Hough transform: if we set σ = 0, we do not obtain meaningful results, due to severe quantization artifacts. These artifacts arise because the gradient direction of a quantized image is not uniformly distributed as there is only a finite number of possible discrete orientations, and some orientations turn out to be much more likely than others [START_REF] Desolneux | Dequantizing image orientation[END_REF]. When computing the Hough transform, this non-uniformity appears as a set of concentric circles around the origin of coordinates (Figure 5). Smoothing the input images by a small Gaussian blur dequantizes the image and the artifacts on the transformed image disappear, allowing to correctly locate the maximum.

The parameter σ also allows to control the resolution of the analysis: using a larger σ permits to move the focus of the algorithm from thin strata to thicker structures. The impact of this parameter is discussed in detail in the section Multiscale dip picking.

Local image orientation

The proposed method relies on a point-wise estimation of the local orientation of the image structures. The simplest measure of local image orientation is the image gradient ∇I σ = (∂ x I σ , ∂ y I σ ), which can be computed with a finite difference scheme. The main issue with the image gradient is that it encodes both orientation and direction (sign), which means that averaging gradients may cancel them out. Despite this, the gradient provides valuable information about the image and can be used successfully in the context of dip picking [START_REF] Saito | Method of determining dips and azimuths of fractures from borehole images[END_REF].

A better local image orientation descriptor is the structure tensor [START_REF] Förstner | A fast operator for detection and precise location of distinct points, corners and centres of circular features[END_REF], which is a field representing a local average direction (modulo π) of the gradient vectors (modulo 2π). This is better than the blurred gradient because opposed vectors will vote for the same direction in Hough space (see next section). The structure tensor field S is computed from the image gradients field ∇I σ by smoothing the singular matrix field ∇I σ ⊗∇I σ , namely S := G µ * (∇I σ ⊗ ∇I σ ) where G µ is a Gaussian kernel, and ⊗ is the outer product of a pair of vectors (rank-1 matrix).

Note that the Gaussian smoothing G µ is not interchangeable with the previous smoothing G σ used for removing the quantization artifacts. Indeed, the outer product results in the components I 2

x , I x I y and I 2 y , which are non-linear operations; thus, we cannot remove the first smoothing step by using a stronger smoothing value µ. These two smoothing steps are both needed for different reasons, and if the quantization artifacts are not corrected before, the resulting structure tensor will present artifacts as those shown in Figure 5.

The local orientation is computed from the 2 × 2 matrix S µ (i, j), as the eigenvector associated to its largest eigenvalue, and the strength of this orientation is the largest eigenvalue. These eigenvalues will be used to weight the Hough votes. The pseudo-code for computing the structure tensor and its largest eigenvector is detailed in Algorithm 2.

Hough voting

The purpose of this step is to find the main bundle of sinusoids (the average dip) appearing on the selected window. This is done by means of the Hough transform. The Hough transform [START_REF] Hough | Method and means for recognizing com-plex patterns[END_REF][START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF] Algorithm 2: local-orientation-by-structure-tensor input :

Image I σ : Ω → R input : Structure tensor filtering µ > 0 output : Local orientation field G : Ω → R 2 1 I x ← ∂ ∂x I σ // finite difference derivatives 2 I y ← ∂ ∂y I σ 3 P ← gaussian-blur(I 2 x , µ) 4 Q ← gaussian-blur(I x I y , µ) 5 R ← gaussian-blur(I 2 y , µ) 6 for p ∈ Ω do 7 a b b c ← P (p) Q(p) Q(p) R(p)
// structure tensor at p [START_REF] Illingworth | A survey of the hough transform: Computer Vision[END_REF]) can be used to find instances of objects within a certain class of shapes by a voting procedure. The class of shapes must be described by a set of parameters. This defines a parameter space in which each point correspond to a particular shape. An accumulator is defined in the parameter space to count the number of votes for each particular parameter value. Then a voting procedure is carried out where each image feature increases the accumulator for the shapes compatible with the observation. The object candidates are obtained as local maxima in the accumulator.

8 T ← a + c 9 D ← ac -b 2 10 α ← T + √ T 2 -4D 2 // largest eigenvalue 11 if |b| > 0 then 12 G(p) ← (b, α -a) // eigenvector 13 else 14 G(p) ← (a, 0) 15 return G
Recall that a dip of parameters (a, b) corresponds to a certain sinusoid h a,b (w) as in equation ( 4). This dip is compatible with a local image orientation (u, v) = S(w, h), given by the structure tensor at pixel (w, h), if the vector (u, v) is perpendicular to the curve h a,b (w), namely if the following condition holds

1, h a,b (w) • (u, v) = 0.
(5)

Developing the derivative, this condition becomes

u -a v sin 2πw W + b v cos 2πw W = 0. (6)
This equation determines a straight line in the (a, b)-plane, and the points of this straight line correspond to all the sinusoids compatible with the local gradient orientation (u, v) observed at (w, h).

The basic idea of Hough voting consists in choosing a finite discretization of the (a, b) parameter plane. Then, each point in the image votes for the corresponding straight line in the (a, b) space according to equation ( 6). The vote is weighted by the magnitude of the local orientation descriptor. Once all the points in the image domain have voted, the site on the (a, b) plane with the most votes is selected as a candidate structure. The complexity of this algorithm is O N 1.5 where N is the number of pixels in the images (assuming that the size of Hough space is proportional to N ).

An alternative voting strategy named randomized voting [START_REF] Xu | A new curve detection method: Randomized hough transform (rht)[END_REF] reduces this complexity to O(N ). This technique considers N randomly chosen pairs of points in the image domain. Each pair of points in the image domain will vote for a single point in the Hough domain, namely the intersection of the straight lines corresponding to each of the two image points. The strength of the vote can be a function of the local orientation strengths of the two points, for instance the product, the sum, or the minimum. Alternatively, the vote can be weighted by the norm of the vector product of the corresponding line parameters. We choose the latter so that parallel lines contribute less.

The traditional Hough space for the sinusoid identification problem [START_REF] Torres | Method for determining dip and strike angles in borehole ultrasonic scanning tool data[END_REF][START_REF] Saito | Method of determining dips and azimuths of fractures from borehole images[END_REF][START_REF] Wang | Method of determining planar events from borehole or core images[END_REF] has the three parameters (a, b, c) of a dip z = ax + by + c. The advantage of our proposal is that votes accumulate on a 2D space instead of a 3D one. This is much faster because a smaller number of votes is enough to fill-in the 2D space and to obtain stable and well-localized maxima. Intuitively, this strategy is aggregating the votes for all the heights in the given window; this is natural because the sinusoids are very often nearly parallel and would correspond to vertically aligned points in the 3D space.

The randomized voting strategy is detailed in Algorithm 3. This algorithm fills an image H(a, b) with the Hough votes for each candidate pair (a, b). The table of votes is initialized to zero (step 1). Then N rand pairs of points are randomly sampled on the image domain (steps 3 and 4) and the local orientations at these points are extracted (steps 5 and 6). Next, we find two straight lines 1 , 2 that are compatible with these local orientations (steps 7 and 8). These lines are expressed in homogeneous coordinates which facilitates computing their intersection point as the cross product (step 9). The homogeneous coordinates of the intersection point are normalized to obtain its Euclidean coordinates (steps 10 and 11), which is identified to a site in the Hough domain. Then a weighted vote is cast for that site (steps 13 and 14).

Hough space discretization and smoothing

The Hough space is discretized as a W × W regular grid covering the square [-κ, κ] 2 of the (a, b)-plane. The range of values spanned by this grid is given by the maximum slope of a plane capable of being represented, usually [-1, 1] 2 , which corresponds to 45 • . In order to handle dips with larger slopes a multi-scale approach is applied, as detailed in the section Multiscale dip picking. The scales are vertically down-scaled versions of the original image by powers of two, called "octaves". In the o-th octave the (a, b)-

Algorithm 3: randomized-hough-voting input : Local orientations G : [0, W ] 2 → R 2 input : Maximum values for a dip κ input : Number of random trials N rand output : Hough votes H : [-κ, κ] 2 → R 1 H ← 0 2 for n = 1, . . . , N rand do 3 p 1 , q 1 ← random-point([0, W ] 2 ) 4 p 2 , q 2 ← random-point([0, W ] 2 ) 5 u 1 , v 1 ← G(p 1 , q 1 ) 6 u 2 , v 2 ← G(p 2 , q 2 ) 7 1 ← sin 2π W p 1 , -cos 2π W p 1 , u 1 /v 1 8 2 ← sin 2π W p 2 , -cos 2π W p 2 , u 2 /v 2 9 x, y, z ← 1 × 2 10 a ← x/z 11 b ← y/z 12 if (a, b) ∈ [-κ, κ] 2 then 13 H(a, b) ← H(a, b) + 1 × 2 plane extends to [-2 o , 2 o ] 2 .
We need to localize the dips very accurately in the (a, b) space, thus in practice we chose the same size W × W as the input window. The Hough space is used (Algorithm 1, step 7) to find the main dip direction on a window as the site (a, b) with the maximum number of votes. To improve the localization of this maximum, the Hough image H(a, b) is smoothed with a Gaussian kernel of standard deviation η (which is a parameter of the method). Due to this smoothing, the actual resolution W × W of the Hough size is not a critical parameter as long as η is larger than the pixel size used to discretize H.

The a contrario theory

The a contrario theory [START_REF] Desolneux | Meaningful alignments[END_REF](Desolneux et al., , 2008) ) is a statistical framework used to set detection thresholds automatically in order to control the number of false detections. It is based on the non-accidentalness principle [START_REF] Witkin | On the role of structure in vision[END_REF][START_REF] Lowe | Perceptual organization and visual recognition[END_REF] which informally states that an observed structure is meaningful only when the relation between its parts is too regular to be the result of an accidental arrangements of independent parts. In the words of D. Lowe, "we need to determine the probability that each relation in the image could have arisen by accident, P (a). Naturally, the smaller that this value is, the more likely the relation is to have a causal interpretation" (Lowe, 1985,p. 39).

A stochastic background model H 0 needs to be defined, where the structure of interest is not present and can only arise as an accidental arrangement. For example, some geometrical feature detection methods such as line segments or elliptical arcs are based on the orientation of the image gradient (Grompone von [START_REF] Grompone Von Gioi | LSD: A fast line segment detector with a false detection control[END_REF][START_REF] Pȃtrȃucean | Joint a contrario ellipse and line detection[END_REF]. In such cases, the background model H 0 assumes that the gradient orientations at each pixel are independent random variables, uniformly distributed in [-π, π). Under this background model, a region of the image where the gradient orientation follows a regular structure would be a rare accident and is detected as such.

We also need to define a family of events of interest T . For feature detection the family of events is the set of all the geometrical events considered, i.e., all the line segments and elliptical arcs considered in the image domain. Then, we need to assess the accidentalness of a candidate feature. For example, if a line segment is present in an image, the gradient orientation at the corresponding position would be orthogonal to the line segment. Then, given a candidate line segment, one measures how well the image gradient corresponds to the candidate event, and evaluates the probability of observing such a good agreement by chance in the background model. A rough agreement could arise just by chance and thus does not correspond to an interesting event; conversely, a very good agreement would be rare and suggest the presence of a structure instead of just a lucky accident. In other words, when this probability is small enough, there exists evidence to reject the null hypothesis and declare the event meaningful. However, one needs to consider that multiple are tested. If 100 tests were performed, for example, it would not be surprising to observe among them one event that appears with probability 0.01 under random conditions. The number of tests N T needs to be included as a correction term, as it is done in the statistical multiple hypothesis testing framework [START_REF] Gordon | Control of the mean number of false discoveries, Bonferroni and stability of multiple testing[END_REF].

Following the a contrario methodology [START_REF] Desolneux | Meaningful alignments[END_REF](Desolneux et al., , 2008)), we define the number of false alarms (NFA) of an event e observed up to an error k(e) as:

NFA(e) = N T • P K H0 (e) ≤ k(e) , (7)
where the right hand term is the probability of obtaining in the background model H 0 an error K H0 (e) smaller or equal to the observed one k(e). The smaller the NFA, the more unlikely the event e is to be observed by chance in the background model H 0 ; thus, the more meaningful. The a contrario approach prescribes accepting as valid detections the candidates with NFA < ε for a predefined value ε. It can be shown [START_REF] Desolneux | Meaningful alignments[END_REF](Desolneux et al., , 2008) ) that under H 0 , the expected number of tests with NFA < ε is bounded by ε. As a result, ε corresponds to the mean number of false detections under H 0 . In many practical applications, including the present one, the value ε = 1 is adopted. Indeed, it allows for less than one false detection on an image or on a set of images, which is usually quite tolerable.

A contrario dip detection

The maximum of the Hough transform gives the parameters (a, b) of the best bundle of sinusoids for the current analysis window. But we still need to determine if this detection is meaningful (in case we are observing a win-dow without sinusoids) and then extract the individual dips. For this we propose to apply an a contrario validation (Desolneux et al., 2008). The test for validating individual dips is similar to the one used in Grompone von Gioi et al. ( 2010), but evaluating the image orientation along the sinusoid at each depth instead of along straight lines. For a given borehole image size, it selects a threshold on the number of pixels with gradient perpendicular to the sinusoid being tested. Algorithms 4 and 5 describe the method. The main inputs are the image gradient at each point in the window of interest of size W × H, and the parameters (a, b) of the sinusoid bundle to be validated. There are two parameters to be set: ρ and ε. The parameter ρ is a normalized angular tolerance. For example, a tolerance of 45 • leads to ρ = 45 180 = 0.25. Here, the NFA corresponds to an upper bound of the expected number of dips, with at least as many correctly oriented gradients as the one currently considered, that will be observed by chance per window. The NFA provides a measure of the quality of the detection that is easier to interpret and to threshold than the quantity of well-aligned pixels. The lower the NFA, the less probable it corresponds to an accidental detection, thus the more probable it corresponds to a real dip. Automatic detections are made by keeping candidates with NFA < ε. The parameter ε corresponds to an upper bound to the number of false alarms that one can accept per window. The output of the algorithm is the set L of parameters of the validated sinusoids relative to the window being analyzed, see Figure 6.

The a contrario detection algorithm builds a list L of "accepted dips". The set L is initially empty (line 1 of Algorithm 5). Then, all the depths h in the range of the window will be tested one by one (line 2). For each depth, the number n of pixels along the sinusoid with defined values (i.e., belonging to the pad/flap traces) are counted. Among them, the number k of pixels with the image gradient orientation compatible with the sinusoid are also counted; both n and k are initialized to zero (lines 1 and 2 of Algorithm 4). To count these numbers, all the columns i of the image are evaluated (line 3) and for each one the corresponding coordinate j(i) in the sinusoid at depth h is computed (line 4). Also, the vector N (i) normal to the sinusoid at pixel (i, j(i)) is computed (line 5). If the borehole image is defined at (i, j(i)) as determined by a mask M (line 6), then the pixel provides information that may match or not the sinusoid; n is increased to count the point (line 7). If the angle between the normal to the sinusoid N (i) and the actually observed image gradient ∇I σ (i, j(i)) is less than the tolerance ρ ∈ (0, 1) (line 8), then the pixel is compatible and k is increased (line 9).

With n and k, the NFA can be computed (line 10) by:

NFA = W 2 H • B(n, k, ρ), (8)
where W × W is the size of the discretization grid defined above (section "Hough space discretization"), H is the height of the observed window, and B is the tail of the s ← n // select optimal neighbor 6 until s does not change // iterate until a local optimum is found

Redundancy reduction

After the local optimization step described in the previous section, we find that the output list of sinusoids is very redundant. There are two reasons for that. First, sinusoids that were very close may collapse to the same local minimum. Second, different local minima may be very close.

To avoid these redundancies and to obtain a cleaner list of dips, we propose to apply an exclusion criterion. The criterion works by enforcing the following condition: each pixel in the input image is owned by at most one sinusoid, only the pixels owned by each sinusoid are used to compute its final NFA value. This exclusion criterion is described by Algorithm 7. This is a greedy algorithm that traverses one time the list of all candidate sinusoids by increasing order of their NFA, from the most to the least meaningful detection (recall that the lowest NFA corresponds to the best match). The precise implementation of the algorithm depends on the function exclude-around(G,s). This function finds the pixels around a thin band of fixed width around the sinusoid s, and marks them as undefined in the mask M . Thus, they will not be counted when updating the NFA value (see line 6 of Algorithm 4), just as if they were outside the image domain.

Parameters

Table 1 summarizes the meaning and reasonable values for all the parameters of the proposed method.

The most important parameter is µ of the structure tensor filtering. Increasing µ, the precision of the detection increases, but the localization becomes worse; if it is too large (at around µ = 20) it starts losing structures. In practice, µ is always set to 11.

The second most important parameter is σ. This is the resolution at which the local orientations are computed. Setting it larger than 1.0 removes a lot of noise, aliasing and artifacts in the input images, but then small structures may be lost. Although it could be used to control the level of detail in the image, in practice it is always set to 1.0. The parameter κ selects how much of the Hough space we want to consider. A value of κ = 1.0 allows for dips at a slope of 45 • . If higher slopes are required, it is better to use a higher level of the pyramid than to increase κ. Decreasing κ may improve the precision of near-horizontal dips, at the price of ignoring high slopes. In practice, κ is always set to 1.0.

The value of the parameter η is not very important, provided it should be large enough to allow a rough estimation of the maximum in the Hough transform. In practice, η is always set to 30.0.

The value ρ of the validation step controls the balance between being tolerant to imprecise sinusoids at the cost of requiring larger observed parts of the sinusoids. The value 0.25 results in a good balance. Even if in some images a slightly smaller or slightly large value improves a little the result, in practice it is always set to ρ = 0.25. The value of the parameter is not very important and does not have an important effect in the final result. In practice, is always set to 1.0.

The parameter N rand affects linearly the running time of the algorithm. It should be as large as possible while being practical. In practice N rand = 10 6 provides a good balance between quality of the result and computation time.

MULTISCALE DIP PICKING

The basic data structure for exploring such large images is the multi-scale pyramid. This data structure is straightforward: it is the union of several vertically down-scaled versions of the original image by factors that are successive powers of two, called "octaves". Notice that this notion of octave is different from the one commonly used in image processing: here we only perform a zoom-out in the vertical direction, thus all the octaves preserve the image width.

As observed in Figure 7, when the borehole direction is nearly parallel to the geological strata, the image contains Figure 8: A contrario dip picking (with local refinement) on a synthetic image. In the left column the detections are highlighted (red for dark-to-bright and green for brightto-dark detections) over a grayscale version of the image which is depicted on the right column in false color (the resistivity is single channel). The center column illustrates the smoothed Hough voting space, with the location of the maximum highlighted with a green dot. The three rows correspond to different dip configurations. In the second and the third row we see two dips in the same analysis window.

nearly vertical structures. These patterns are actually sinusoids of very large amplitude, as can be seen by looking at the higher octaves. Our algorithm correctly finds these sinusoids once they appear entirely inside the selected region of interest.

Automatic dip picking

The proposed algorithm can be used to process a whole borehole automatically. Two parameters need to be specified: a window size ω and the number of spatial scales O. The borehole image is processed at full resolution and at O -1 consecutive dyadic spatial scales. At each scale, the image is processed by windows of height ω, advancing by half-window steps of ω/2. All the dips found at different windows and scales are put together. This can generate repeated dip detections, to filter them we apply the exclusion criterion described above.

In Figure 8 we show a typical run of our algorithm on three windows of a very well behaved synthetic image. In the first window, the correct dip is detected as the position of the maximum in Hough space. In the second window, when there is more than one dip (e.g., near a "cross-bedding"), there are two competing maxima in Hough space, and the one corresponding to the dip with largest area on the region of interest is chosen. In the last window, the second maximum becomes larger and the other dip is detected. It is the task of the exclusion criterion to merge the whole set of dips without redundancy.

Computational cost

The first steps of the method all have a computational cost proportional to the number of pixels in the image or window being processed. This includes the extraction of the window, the filling-in of missing data, the Gaussian filtering and the computation of local orientations. The randomized Hough transform step has a computational cost proportional to the number of random trials N rand . The a contrario validation and refinement step have a computational cost proportional to the number of candidates being evaluated, thus proportional to the number of windows evaluated. All in all, the computational complexity is O(N ), where N is the number of pixels in the borehole image.

In practice, the randomized Hough transform is very effective and the algorithm computes the result in realtime when applied to windows in a semi-supervised human/computer interactive interface. When applied automatically, a non-optimized code can process a whole borehole of 200,000 pixel rows, corresponding to about 1450 m (4800 feet), in less than 30 s on a modern laptop computer.

EXPERIMENTS

In this section we present experiments for evaluating the performance of the proposed algorithm. First we will carry out a quantitative evaluation comparing against manually annotated dips on three borehole images. Then we will perform an ablation study to analyze the impact of some components of the proposed algorithm. Lastly, we present a qualitative analysis of the results and discuss its limitations.

Quantitative evaluation

To evaluate the proposed method we shall compare against three manually annotated borehole images. These datasets (denoted A, B and C) are described below: The proposed multiscale algorithm was applied to the three borehole images using five octave scales and default parameters for the detection. This yielded a list of detections each with an associated NFA (which can be seen as a confidence measure).

To compare the manual annotations and our results we analyzed the dips on a rolling window along the well. We set the size of the rolling window to 6 m (20 ft) for the second and third boreholes (B and C), while for the first one (A) we used 0.6 m (2 ft) because the density of events is much higher for such a short borehole. The results can be seen in Figure 9. We computed and show the density of dips within each window, the dominant azimuth angle (computed as the mode of the azimuths), and the mean of the dip inclination. From the plots we can observe that the proposed method produces meaningful dips with parameters close to the manual annotations. As manual labeling is not exhaustive, it is to be expected that the proposed method can produce many more dip detections.

To further assess the pertinence of the detected dips we evaluated the error with respect to the manually annotated ones. For each manually annotated dip, we searched for the closest detected dip in terms of root mean squared error (RMSE) between the two graphs

d(h 1 , h 2 ) = 1 W w=1,...,W |h 1 (w) -h 2 (w)| 2 . ( 10
)
Note that this error has distance units and can be used to set a threshold to identify false-positive detections. Here we considered all the correspondences to analyze the differences in the parameters of the two dips. We computed the difference of depth, azimuth, and inclination of the matched dips. The distributions of these errors for each one of the annotated borehole images are shown in Figure 10. We observe that all the histograms of errors are concentrated around zero. This means that close to most of the manual annotations, there was a detection with similar parameters. In the borehole B we observe a higher error, this is likely due to higher density and variability of events present in this borehole.

Ablation study

We performed three simple "ablation" experiments to assess the validity of the method by comparing its performance to manual dip annotations by expert geologists. Unfortunately, we cannot compare the performance of the proposed method to other publicly available methods, as we could not find any. More precisely, we studied the impact of three crucial steps of the pipeline: the gradient dequantization preprocess (controlled by the parameter σ), the multi-scale selection scheme (controlled by the parameters o), and the local dip refinement (implemented in Algorithm 6).

For that, we will define a true positive measure based on the RMSE error between matched dips of equation ( 10). This error measure has distance units: the distance between two sinusoids is defined as the average distance of the corresponding ellipses on the boundary of the cylin- der. We will say that a manually annotated dip is correctly found by an algorithm if the error with the matched dip is smaller than 30.48 cm (1 feet). Note that we cannot easily define a measure for false positives because the manually annotated dips are not exhaustive (large swaths of the image were not annotated). This true positive score allows us to evaluate the performance of variants of the proposed algorithm obtained by removing different components. In Table 2 we present the results of these experiments. We consider reducing the number of scales (using only the first octave, the first two, or all the 5 octaves), removing the dequantization blur (σ = 0 for removing it), and removing the iterative local dip refinement (by setting the refinement iterations to 0 and 3). All these variants lead to a reduced performance, most notably removing the multiscale and the iterative refinement.

Qualitative evaluation and discussion

Let us now see results of our algorithm on some examples (see Figures 11 and12). In all cases, the initial dip (a, b) was computed by the Hough transform algorithm (line 6 of Algorithm 1), and the parameters (a, b) were used as input for the NFA computation. Notice that in most cases the correct dip is obtained, even for fragmented sinusoid as seen in Figure 11 The sequence of detections on Figure 12 shows the behavior of the algorithm around a fracture and a crossbedding. Notice that, for each position of the sliding window, a single dip is selected by the Hough transform, the one that is present on the largest part of the window. This is the slope (a, b) that is given to the a contrario refining step, which determines the set of depths c where sinusoids are detected. For this particular example, the detections flip from one side to the other of the cross-bedding, which happens between the second and third images in the sequence.

scales

Limitations

The proposed method is far from perfect and there is room for improvement. This can be seen in the quantitative evaluation were the performance on the dataset B is below the other two. This borehole is hard as it presents many events, this can be seen in the annotation in Figure 9, see in particular at depth 3627 m (11900 ft).

One can see some examples of wrong detections in Figure 12. There is an erroneous dip by the lower part of the second image. In this case, the Hough step proposed the wrong candidate due to the cross-bedding. Then, the refining step did not explore enough variations to fit precisely the actual dip. In addition, due to the blurred aspect of the edges, the gradient orientation roughly match the wrong dip, producing a barely meaningful detection. In most cases, this problem is solved by the exclusion criterion when all the dips at all windows and scales are put together; indeed, the dips detected in the third image for the same part of the image have better NFA value and will be preferred. But there are still cases where this strategy fails. A similar example of wrong detection, although less extreme, can be seen in the middle of the second image.

CONCLUSIONS

We have proposed a method to perform fast and automatic dip picking on borehole images, amenable to an embedded online implementation due to its lower computational requirements. The method is based on an a contrario criterion to control the number of false detections, coupled with a search space optimization based on Hough spaces. Our experiments show that the algorithm detects correctly the evolution of the dip along the borehole depth, and gives no dip when there are no apparent structures in the data. A quantitative evaluation on annotated borehole images shows that the proposed method can detect from 60% to 90% of the dips annotated by an expert in a completely unsupervised way. The same method can also be integrated in an interactive dip annotation tool to assist the annotation.

Figure 12: A contrario dip picking, on a sliding window traversing a change of inclination (from column 2 to 3) and a cross-bedding with no change of inclination (from column 3 to 4).

Figure 2 :

 2 Figure 2: (a) Cylinder of radius R. (b) Intersection of the cylinder with a plane. (c) When the cylinder is unrolled, this ellipse becomes a sinusoid. See equation (4) for the equation of a sinusoid h(w) in terms of its Hough parameters (a, b). (d) Sinusoid in Hough (a, b) coordinates. (e) Dip detection in a synthetic image (red for dark-to-bright and green for bright-to-dark detections). (f ) Corresponding votes in Hough space.

Figure 3 :

 3 Figure 3: Workflow of the proposed method.

Algorithm 1 :

 1 dip-picker input : Image J of size W × H input : Valid data mask M (W × H) input : Location h ∈ {1, . . . , H} input : Scale o ∈ {1, . . . , log 2 H } input : Positive parameters σ, µ, η, ρ, , κ, N rand output : List of sinusoids L

Figure 4 :

 4 Figure 4: Example of Filling-in produced by the proposed Laplacian interpolation method. This creates a smooth interpolation of the data, which will not perturb the computation of the local orientation.

Figure 5 :

 5 Figure 5: Importance of image dequantization. These images show the Hough transform of the same input image, pre-filtered with different Gaussian kernels of the indicated sigma. The circular artifacts are due to the quantized gray levels of the image. They interfere severely with the location of the maximum, which is indicated as a red circle in the images. The transform is dequantized by applying a small blur to the input images.

Figure 6 :

 6 Figure 6: Example of dip picking results using the proposed method. (a) shows the input window. (b) shows the extracted Hough space with its detected maximum highlighted in yellow. (c) shows the validated detections highlighted (red and green denote opposed orientations of the sinusoid) over a grayscale version of the image.

Algorithm 7 :

 7 refine-dips-with-exclusion input : Image I σ (W × H) input : Valid data mask M (W × H) input : Number > 0 (NFA threshold) input : L list of sinusoids to be refined output : L list (possibly shorter) of refined sinusoids 1 L ← ∅ // output accumulator 2 L ← sort-sinusoids-by-growing-nfa(L, I σ , M ) 3 for s ∈ L do // traverse sinusoids in order 4 s ← refine-dip(I σ , M, s) // refine this sinusoid 5 if dip-nfa(I σ , M, s ) < then // if meaningful 6 L ← L ∪ {s } // add it to list 7 M ← exclude-around(M, s ) // and mark neighboring pixels as undef.

Figure 7 :

 7 Figure 7: Results of the multiscale analysis. The same part of the well is analyzed at octaves 0 (original resolution, input (a) and output (b)) and 3 (vertical zoom-out of factor 8, input (c) and output (d)). Notice that some very "tall" sinusoids only become visible at high enough scales. Depths are in meters.

•

  A: 221 manually annotated dips along 30 m (97 ft) of borehole, from 4078 m to 4107 m (13380 to 13476 ft); • B: 3766 manually annotated dips along 1485 m (4873 ft) of borehole, from 3110 m to 4596 m (10205 to 15079 ft); • C: 194 manually annotated dips along 3345 m (1099 ft) of borehole, from 3037 m to 3372 m (9964 to 11063 ft).

Figure 9 :Figure 10 :

 910 Figure9: Analysis of dips on a sliding window. The three rows correspond to the three evaluation boreholes A, B and C with manual annotation. The plots in each column show: the number of dips on the rolling window, the dominant azimuth angle (which is computed as the the mode of the azimuths in the window), and the mean inclination within the window.

  True positive rate for different variant of the algorithm evaluated on the three annotated borehole images A, B and C. The first three columns correspond to options of the proposed dip-picking algorithm: number of octaves considered in the multiscale scheme, dequantization parameter σ, and number of iteration of local dip refinement.

Figure 11 :

 11 Figure 11: A contrario dip picking on two challenging zones, one of them (a) with low contrast, the other (b) suffering from pad misalignment. In both cases correct dips are being identified, with just a couple of false alarms on the first one. Depths are in meters.

Table 1 :

 1 Summary of the meaning and reasonable values for all the parameters of the proposed method.

	parameter default	meaning
		value
	σ	1.0	pre-filtering (dequantization)
	µ	11.0	structure tensor filtering
	κ	1.0	maximum values for a dip (a, b)
	η	30.0	blur of the Hough space
	ρ	0.25	normalized angular tolerance
		1.0	number of false alarms
	N rand	10 6	number of random trials

Table 2 :

 2 dequant. σ ref. iter. Ablation study.

				A	B	C
	all	1	100	89% 59% 88%
	1 oct.	1	100	28% 34% 28%
	2 oct.	1	100	81% 48% 76%
	all	0	100	84% 52% 85%
	all	1	3	86% 56% 85%
	all	1	0	72% 45% 71%

I ← extract-and-resample-window(J, M, h, o)

I f ← fill-in-missing-data(I, M )

I σ ← gaussian-blur(I f , σ)

G ← local-orientation(I σ , µ)

H ← randomized-hough-voting(G, N rand , κ)

H η ← gaussian-blur(H, η)

(a, b) ← arg max H η (x, y)

L ← locate-dips(I σ , M, ρ, ε, a, b)

L ← refine-dips-with-exclusion(I σ , M, ρ, ε, L)

return L

ACKNOWLEDGEMENT

We want to thank Jean-Michel Ghidaglia for his warm encouragement towards publishing this work. We are also thankful for the support and fruitful discussions with Schlumberger experts: Isabelle Le Nir, Carlos Maeso, Tetsushi Yamada, Nicolas Veron, Alexis He, and Kristina Prokopetc. This project was carried out in the framework of the French CITEPH (EVOLEN) collaborative program under reference CITEPH-40B-2014. This work contains information provided by the Oil and Gas Authority and/or other third parties. We also thank the editors and the three anonymous reviewers for the valuable comments and suggestions.

DATA AND MATERIALS AVAILABILITY

The dataset of three manually annotated borehole images and the supplementary material associated to this paper are available at https://cmla.github.io/dip-picker/.

binomial distribution, defined as:

The binomial tail can be computed efficiently using the incomplete beta function.

Finally, if the NFA of the sinusoid being evaluated is smaller than ε (line 4 of Algorithm 5), the corresponding sinusoid (a, b, h) is added to the output set L (line 5). Figure 6 shows an example of the detected sinusoids.

Note that if there are no sinusoidal structures (compatible with a, b), then the previous algorithm should not detect any dip. If there are several candidate parameters (a, b), the algorithm is applied for each pair. If no dip is detected for any candidate plane orientation, then the window is declared to be without sinusoids.

Local dip refinement

The sinusoids detected by the previous step may not be very accurate because their (a, b) parameters are com-puted globally for the whole window (using the Hough transform). Indeed, the NFA test (Algorithm 4) reduces the set of candidate sinusoids in the current window by determining whether a sinusoid is meaningful or not. However, all of these sinusoids have the same parameters (a, b).

The goal of the refinement step is to slightly adjust the three parameters of each sinusoid so that it better fits to the local geometry of the image. In what follows we propose a method to iteratively refine the parameters of each individual sinusoid by a simple discrete coordinate descent optimization (see Algorithm 6). The refinement described in Algorithm 6 uses the NFA score as a measure of fit. It tests small variations of the parameters and selects the one yielding the smallest NFA value. This process is iterated until no improvement of the NFA can be achieved. Notice that the potentially infinite loop in this algorithm halts because the search space is finite (it is restricted to sinusoids that intersect the current window).