2.1 가우시안 프로세스를 이용한 격자의 점유 확률 예측

GP 격자 지도는 환경의 각 위치 사이의 공간적 연관성을 GP 확률 분포로 가정하고, 각 격자의 점유 확률을 가우시안 프로세스 회귀(GPR, Gaussian Process Regression)를 이용해서 평균과 분산으로 나타낸 것이다^[10]. GP를 이용한 특정 위치의 점유 확률은 격자 구조에 제한되지 않고, 관심 영역의 어느 지점에 대해서도 연속적으로 추정될 수 있다. 그러나 GP 결과를 격자 구조로 나타낼 경우 환경의 형상 파악과 활용에 용이하기 때문에 이동 로봇을 위해서는 GP 격자 지도를 작성한다.

GP 격자 지도를 구성하는 각 격자의 GP $y({x})$는 평균 함수(mean function) $\mu({x})$와 공분산 함수(covariance function) $k({x},\: {x}')$로 정의된다.

${ x}$와 ${ x}'$는 격자의 2차원 위치인 $2\times 1$ 벡터이며, 타겟(target) 값인 점유 상태 $y$는 -1에서 1 사이의 값을 가진다. -1은 비점유 격자를 의미하며, 1은 점유 격자를 의미한다. $n$개의 학습 격자 $D=\left\{\left({}{x}^{i},\:{y}^{{i}}\right)|1\le{i}\le{n}\right\}$가 있을 때, 테스트 격자의 위치 ${ x}^{*}$의 점유 상태에 대한 가우시안 분포의 평균과 분산은 다음과 같이 예측된다.

${ X}$는 $2\times n$ 행렬 $\left[{x}^{1},\: \ldots ,\: {x}^{n}\right]$이고, ${bold y}$는 $n\times 1$벡터 $\left[y^{1},\: \ldots ,\: y^{n}\right]^{T}$이다. $k({ X},\: { X})$는 모든 학습 격자들 사이의 $n\times n$ 공분산 행렬이고, $k\left({ X},\: {x}^{*}\right)$는 모든 학습 격자와 테스트 격자로부터 계산된 $n\times 1$ 공분산 벡터이며 $k\left({ X},\: {x}^{*}\right)= k\left({x}^{*},\: { X}\right)^{T}$이다. $\sigma_{n}^{2}$는 독립적이고 동일하게 분포된(independent and identically distributed) 가우시안 분포로 가정한 측정값의 노이즈를 의미한다.

GP의 공분산 함수로 사용되는 여러 커널 중에서 Matérn 공분산 함수는 공간의 기하학적 연관성을 모델링하기에 적절하다고 알려졌다^[11,12]. Matérn 공분산 함수는 감마(gamma) 함수와 Bessel 함수로 정의되는데, 평활도(smoothness)를 제어하는 파라미터 $\nu$가 $3/2$일 때, 다음과 같이 표현된다.

하이퍼파라미터(hyperparameter) $\sigma_{f}$와 $\rho$는 각각 signal variance와 characteristic length scale로 불리며, $\sigma_{f}$는 커널 계산 결과의 크기(amplitude)를 결정하고 $\rho$는 두 위치의 연관성(correlation)이 간격에 따라 얼마나 급격히 변하는지를 제어한다.

테스트 격자에 대해 예측된 $y\left({x}^{*}\right)$의 범위는 $[-1,\: 1]$이기 때문에, 이를 $[0,\: 1]$ 범위의 확률값으로 변환하여 점유 확률 $m\left({ x}^{*}\right)$를 계산한다.

위 수식은 로지스틱 회귀(logistic regression)를 이용한 수식이며, $\gamma > 0$는 시그모이드(sigmoid)의 형태를 결정하고 $\sigma_{\min}^{2}$는 예측된 분산들의 최솟값이다.

그림 2. 각 센서 레이를 따라 샘플링된 학습 데이터를 이용한 비점유 영역의 가우시안 프로세스 회귀 결과

Fig. 2. Gaussian Process Regression(GPR) for unoccupied regions using sampled training data per each sensor ray

2.2 로컬 커널을 이용한 가우시안 프로세스 격자 지도 근사화

GP 격자 지도의 점유 확률 예측 성능을 결정하는 것은 공분산 함수를 정의하는 최적화된 하이퍼파라미터와 학습 데이터이다. 하이퍼파라미터는 주변 우도(marginal likelihood)의 네거티브 로그(negative log)를 최소화하는 과정을 통해서 최적화될 수 있다^[12]. 학습 데이터는 최적화 과정과 테스트 격자의 평균 및 분산 예측에 사용되기 때문에, 적절한 학습 데이터 선택은 GP 격자 지도의 예측 성능에 주요한 영향을 미친다. 또한, 예측 과정에서 역행렬 연산이 수행되기 때문에(식 (3)과 (4)), 계산 속도는 학습 데이터의 양에 크게 좌우된다. 기존에 연구된 거리 센서를 이용한 GP 격자 지도 작성 방법들은 학습 데이터를 획득하기 위해서 각 거리 측정값의 빔(beam) 방향 즉, 센서 레이(sensor ray)를 따라서 샘플링을 수행하였다^[9,14,15]. 점유 영역에 대한 학습 데이터의 양은 라이다 측정값의 스캔 데이터 개수로 한정된다. 그러나, 비점유 영역에 대한 학습 데이터의 양은 샘플링 주파수에 스캔 데이터 개수를 곱한 값이기 때문에 샘플링 주파수에 따라 크게 달라지며 점유 영역 학습 데이터보다 훨씬 많다. 그림 2는 센서 레이를 따라 샘플링된 학습 데이터를 사용한 비점유 GP 작성 결과를 보여준다. 그림 2(a)는 라이다 측정값으로 작성한 격자 지도이다. 그림 2(b)-(d)는 각각 센서 레이를 따라 샘플링을 최대 10, 20, 50회 수행하고 최적화된 하이퍼파라미터들($\sigma_{f}^{2}$, $\rho$, $\sigma_{n}^{2}$)을 적용한 비점유 GP 예측 결과이다. 진한 파란색일수록 점유 확률이 낮은 비점유 영역에 해당된다. 테스트 위치의 예측 결과가 학습 데이터로부터의 거리에 크게 영향을 받는 GPR의 특성에 의해서, 샘플링 주파수가 작은 경우 샘플링이 수행되지 않은 비점유 영역의 점유 확률이 샘플링이 수행된 비점유 영역의 점유 확률 보다 비교적 높은 것을 확인할 수 있다. 이에 따라 서로 인접한 비점유 영역임에도 불구하고 확률값의 차이로 인해 불필요한 패턴이 생성되었다. 이를 방지하기 위해 샘플링 주파수를 높이면 데이터 양이 지나치게 증가하기 때문에, 다음 탐사 목적지를 실시간으로 결정해서 주행해야 하는 상황에 적합하지 않다. GP의 하이퍼파라미터 혹은 로지스틱 회귀의 형태를 변경하면 샘플링 지점과 예측 지점의 거리에 따른 평균과 분산의 변화 경향이 달라지지만, 서로 인접한 비점유 영역에서 발생하는 샘플링 유무에 의한 확률 변화의 차이를 방지하기는 어렵다.

본 논문에서 제안하는 로컬 커널 기반의 근사화된 GP 격자 지도 예측은 크게 세 단계로 진행되며, 구체적인 순서는그림 3의 Algorithm 1과 같다. Algorithm 1은 현재 획득된 라이다 측정값 ${z}$와 이로부터 작성된 점유 격자 지도 $grid_{{local}}$를 이용해서 근사화된 GP 격자 지도 $map_{{GP}}$를 작성한다. 첫 번째 단계는 라이다 측정값을 이용하여 학습 데이터를 추출하고 하이퍼파라미터를 최적화하는 단계이다(line 1-2). ${ z}$는 $K$개의 스캔 데이터로 구성되며, 각 스캔 데이터 $z^{k}$ 는 측정 방향 $\alpha^{k}$에 존재하는 장애물까지의 거리 정보 $r^{k}$를 가지고 있다. 학습 데이터 ${ X}_{{obs}}$는 이들을 2차원 위치로 변환한 $2\times K$ 행렬이고, ${ y}_{{obs}}$는 모든 요소가 점유 상태를 의미하는 값 1인 $k\times 1$ 벡터이다.

그림 3. 제안된 로컬 커널을 이용한 근사화된 가우시안 프로세스 격자 지도 작성 알고리즘

Fig. 3. Algorithms to build the approximated Gaussian Process occupancy grid map using the proposed local kernel

두 번째 단계에서는 점유 영역과 비점유 영역에 대한 GPR을 각각 수행한다(line 3-7). $grid_{{local}}$를 구성하는 모든 격자의 인덱스(index)들을 2차원 위치들로 변환하여 테스트 데이터 ${ X}^{*}$를 획득한다(line 3). 점유 GP는 기존의 GP 예측 방식에 따라, 식 (3)과 (4)에 점유 영역의 학습 데이터(${ X}_{{obs}}$, ${boldy}_{{obs}}$)와 ${bold X}^{*}$를 대입하여 예측된다(line 4). 비점유 GP는 로컬 커널을 이용한 GP 근사화를 이용해서 예측되며(line 5), 구체적인 순서는 그림 3의 Algorithm 2와 같다. 예측된 점유 GP와 비점유 GP에 로지스틱 회귀를 적용하여 0~1의 확률값으로 변환된 점유 GP의 평균 $p_{{obs}}^{*}$과 비점유 GP의 평균 $p_{{emp}}^{*}$을 얻는다(line 6-7).

마지막 세 번째 단계에서는 점유 GP와 비점유 GP를 융합하여 GP 격자 지도를 획득한다. 이를 위해 서로 다른 학습 데이터에 의한 예측값을 통합할 수 있는 Bayesian Committee Machine(BCM)을 사용하며^[16], 융합된 GP의 평균 $\mu_{{GP}}$과 분산 $\sigma_{{GP}}^{2}$은 다음과 같다(line 8).

여기에 다시 로지스틱 회귀를 적용하여, 각 격자가 0~1사이의 점유 확률을 나타내는 GP 격자 지도 $map_{{GP}}$를 얻는다(line 9).

앞서 언급한 바와 같이 비점유 GP는 본 논문에서 제안하는 로컬 커널을 이용한 GP 근사화를 이용해서 예측된다. 기존의 방법들처럼 샘플링에 의해 추출된 비점유 영역의 학습 데이터를 식 (3)과 (4)에 대입하여 계산하는 것이 아니라, 라이다 측정값을 이용해서 작성된 격자 지도 $grid_{{local}}$를 학습 데이터로 활용하는 합성곱(convolution) 연산을 통해 근사화된 비점유 GP를 계산한다. 그림 3의 Algorithm 2는 $M\times N$ 크기의 $grid_{{local}}$를 이용한 비점유 GP 예측 과정을 보여준다. 최적화된 하이퍼파라미터를 이용하여 로컬 커널 ${k}_{{emp}}$을 생성한다(line 2). 이 로컬 커널과 학습 데이터로 변환된 $grid_{{local}}$의 합성곱 연산을 수행하여 근사화된 비점유 GP를 계산할 수 있다. ${k}_{{emp}}$의 크기가 $n\times n$일 때, 커널의 각 요소 ${k}_{{emp}}(p,\: q)$는 다음과 같이 계산된다.

$grid_{{local}}$을 학습 데이터로 변환하기 위하여, 비점유 격자에는 낮은 확률값 $p_{{emp}}$을 부여하고(line 4-5), 비점유 격자가 아닌 격자(점유 격자와 미지의 격자)에는 $p_{{emp}}$ 보다 큰 확률값 $p_{{non}-{emp}}$을 부여한다(line 7). 이는 비점유 격자와 다른 격자들이 서로 구분되도록 하기 위함이다. 본 논문에서 $p_{{emp}}$는 0.1로 설정하였고, $p_{{non}-{emp}}$는 미지의 격자의 점유 확률 $th_{m{unknown}}$과 동일한 0.5로 설정하였다. 앞서 작성한 로컬 커널 ${k}_{{emp}}$와 변환된 $grid_{{local}}$의 합성곱 연산을 통해 비점유 GP의 평균과 분산을 예측한다(line 10-13). 이와 같은 2차원 합성곱 연산은 픽셀-와이즈로 수행되더라도 각 연산에 필요한 메모리 및 계산 부하가 적기 때문에, 고차원 역행렬 계산이 요구되는 기존의 방법보다 더 빠른 계산 속도를 기대할 수 있다. 그림 4는 로컬 커널을 이용한 비점유 GP 근사화 결과를 보여준다. 기존 샘플링 방식의 예측 결과(그림 2(b)-(d))와 달리 비점유 영역의 확률값이 거의 일관적으로 낮은 것을 확인할 수 있다.

그림 4. 로컬 커널을 이용한 비점유 GP 근사화 결과

Fig. 4. GPR for unoccupied regions using the proposed local kernel-based approximation.

2.3 가우시안 프로세스 격자 지도를 이용한 프론티어 확률 계산

GP 격자 지도는 각 격자가 점유 확률을 나타내는 평균과 분산 값을 가지고 있고, 이들을 이용해서 각 격자의 프론티어 확률을 예측할 수 있다. GP 격자 지도 상에서, 인접한 격자와의 점유 확률 차이가 큰 에지 격자들 중 점유 격자들을 제외한 나머지 에지 격자들의 프론티어 확률이 커야 한다. 이런 프론티어 격자의 특성을 반영하여, 아래와 같이 프론티어 확률을 계산할 수 있다^[8].

프론티어 확률은 GP 격자 지도 $p_{{GP}}$와 점유 GP $p_{{obs}}$의 그래디언트(gradient)를 이용해서 계산되며, $(i,\: j)$는 각 격자의 인덱스이다. $\left . ∥\nabla p_{{GP}}\right .∥_{1}$은 모든 에지 격자에서 큰 값으로 계산되며, $\left . ∥\nabla p_{{obs}}\right .∥_{1}$은 점유 격자에서 큰 값으로 계산된다. $\beta$는 점유 격자의 영향을 제어하는 상수로, $\beta$가 클수록 점유 격자에 의해 프론티어 확률이 낮아지는 격자의 범위가 넓어진다.

3.1 비점유 영역에 대한 가우시안 프로세스 계산 시간 비교

실내 환경을 주행하면서 획득한 2차원 라이다 측정값을 이용해서, 비점유 GP 작성에 소요되는 계산 시간을 비교하였다. 사용된 라이다는 Hokuyo의 UTM-30LX이며, 270o 범위를 30 m까지 스캔할 수 있다. 라이다가 탑재된 로봇 플랫폼은 차동 구동 방식(differential drive)으로 이동하는 Omron의 Pioneer 3-DX이다. 로봇이 0.3 m 이동할 때마다 라이다 스캔 데이터를 저장하였고, 또한 이를 이용해서 로컬 격자 지도를 작성하였다. 그림 7(c)의 사진들은 각각 빨간색 화살표 방향에서 촬영한 주행 시작 지점과 로봇, 그리고 주행 종료 지점이다. 로컬 격자 지도 작성을 위해 Mobile Robot Programming Toolkit (MRPT)^[17]라이브러리를 사용하였고, 각 격자의 크기는 0.05 m로 설정하였다. Cartographer Simultaneous Localization and Mapping(SLAM)^[18]을 이용해서 각 라이다 측정값이 획득되었을 때의 로봇의 위치를 추정하였다. 센서 측정값 획득 이후, 제안된 알고리즘의 성능 검증을 위한 GP 격자 지도 학습 및 프론티어 확률 계산 알고리즘은 Matlab 2024a로 작성되었으며, GP 예측을 위해서 Gaussian Processes for Machine Learning 라이브러리^[19]를 이용하였다. 실험은 Intel i9-14900 (5.8 GHz) CPU와 64 GB 메모리가 탑재된 노트북에서 수행되었다. 환경 전 영역에 대한 지도 작성이 완료될 때까지 주행하며 획득한 290개의 라이다 측정값과 해당되는 로컬 격자 지도를 이용해서, 본 논문에서 제안된 로컬 커널 기반 GP 근사화 기법과 기존의 학습 데이터 샘플링 방식의 비점유 GP 계산 시간을 비교하였다. 두 경우 모두 테스트 데이터는 로컬 격자 지도를 구성하는 각 격자의 2차원 위치이기 때문에, 테스트 데이터의 범위와 개수는 동일하다.

GP 지도 작성과 프론티어 확률 계산을 위해 사용된 파라미터는 표 1과 같다. 또한, 비점유 GP 근사화를 위해 사용된 커널 ${k}_{{emp}}$의 크기는 11×11로 설정하였다.

그림 5(a)는 학습 데이터 샘플링 방식의 비점유 GP 계산 시간을 보여준다. 여기에서 비점유 GP 계산 시간은 학습 데이터와 테스트 데이터를 이용하여 식 (3)과 (4)를 계산하고, 로지스틱 회귀를 수행하는 데에 소요된 시간이다. 빨간색 그래프는 센서 레이를 따라 최대 10번 샘플링을 수행하여 학습 데이터를 생성한 경우의 비점유 GP 계산 시간이며, 녹색 그래프는 최대 20번 샘플링을 수행한 경우의 계산 시간이다. 주행 시작 지점과 종료 지점은 너비 약 2.5 m의 좁은 복도 환경(그림 7(c))이기 때문에, 라이다 측정값들의 수치가 전반적으로 작고 로컬 격자 지도 역시 그 범위와 격자의 개수가 적다. 즉, 비점유 영역의 학습 데이터와 테스트 데이터의 개수가 비교적 적다. 그러나, 넓은 영역에서 획득된 라이다 측정값에 대해서는 학습 데이터와 테스트 데이터의 수가 증가하기 때문에, 로봇이 주행하면서 맞이하는 환경 구조의 변화에 따라 비점유 GP 예측에 소요된 시간도 크게 변하는 것을 확인할 수 있다. 전반적으로 최대 10번 샘플링을 수행한 경우의 소요된 계산 시간이 최대 20번 샘플링을 수행한 경우보다 작다. 최대 10번 샘플링의 경우 비점유 GP 계산에 최대 74.13 초가 소요되었고, 최대 20번 샘플링의 경우 최대 302.71 초가 소요되었다. 또한, 각 계산 시간의 편차가 크다.

그림 5(b)는 제안된 로컬 커널 기반 GP 근사화에 의한 비점유 GP 계산 시간을 보여준다. 즉, 그림 3의 Algorithm 1에서 line 5와 line 7을 실행하는 데에 소요된 시간이다. 합성곱 연산은 단순한 사칙연산의 반복으로 수행되기 때문에 계산 부담이 크지 않다. 합성곱 연산도 픽셀-와이즈로 수행되기 때문에 로컬 격자 지도의 크기가 큰 경우 계산 시간이 증가할 수 있지만, 그 영향이 현저하게 드러나지 않는 것을 확인할 수 있다. 비점유 GP 계산에 최대 1.33 초의 시간이 소요되었으며, 두 번 째 최대 계산 시간은 1.28 초이다. 계산 시간의 평균은 1.24 초이며 표준 편차는 0.01이다. 이 수치는 최대 수 십 초에서 수 백 초에 이르는 기존의 방법에 비해 훨씬 작은 값이며, 실시간 탐사 알고리즘에 적용 가능한 수준이라고 볼 수 있다.

그림 5. 비점유 GP 계산 시간

Fig. 5. Computation time to predict GPR for unoccupied regions

3.2 GP 격자 지도 작성 및 프론티어 확률 계산 결과

제안된 로컬 커널 기반의 GP 근사화를 이용한 GP 격자 지도 작성 및 프론티어 확률 계산 결과를 기존의 샘플링 방식에 의한 결과와 비교하였다. 그림 6은 세 스텝(step)의 라이다 측정값들을 이용한 GP 격자 지도 작성 및 프론티어 확률 계산 결과를 보여준다. 각 그림의 색이 나타내는 수치는 그림 7(a)와 그림 7(b)의 오른쪽에 있는 컬러바(colorbar)의 수치와 같다. 103번째 스텝에서 105번째 스텝까지의 라이다 측정값들을 이용하여 각각의 점유 GP와 비점유 GP를 예측하였고, 해당되는 로봇의 자세 정보와 BCM을 이용하여 세 개의 라이다 측정값들이 결합된 점유 GP와 비점유 GP를 획득하였다. 그림 6(a)는 점유 GP이며, 그림 6(b)는 제안된 GP 근사화에 의한 비점유 GP이다. 그림 6(c)는 점유 GP와 그림 6(b)의 비점유 GP를 융합한 GP 격자 지도이며, 그림 6(d)는 그림 6(c)로부터 계산된 프론티어 확률이다. 즉, 그림 6(b)-(d)는 제안된 비점유 GP 근사화에 의한 결과이며, 그림 6(e)-(f)는 기존 샘플링 방법에 의한 결과이다. 그림 6(e)는 그림 6(a)의 점유 GP와 기존의 샘플링 방식(최대 20회)에 의해 예측된 비점유 GP를 융합한 GP 격자 지도이다. 그림 6(f)는 그림 6(e)로부터 계산된 프론티어 확률이다. 제안된 방법에 의한 비점유 GP를 이용해서 융합된 그림 6(c)의 GP 격자 지도는 서로 인접한 비점유 영역에서 대체적으로 일관성 있는 예측 결과를 보여주며, 장애물에서 멀리 떨어진 비점유 영역이 균일한 파란색으로 표시되어 있다. 또한, 그림 6(d)의 프론티어 확률도 실제 탐사 상태를 반영하여 미지의 영역과 맞닿은 비점유 에지 영역에서 비교적 높은 수치(붉은색 픽셀)를 보여준다. 반면, 기존의 방법에 의한 그림 6(e)의 GP 격자 지도에서는 비점유 GP에 사용된 학습 데이터의 샘플링 영향에 의해 생성된 불필요한 패턴(비점유 영역에서 파란색과 하늘색이 번갈아서 나타남)을 확인할 수 있다. 프론티어 확률 역시 GP 격자 지도의 그래디언트를 기반으로 계산되기 때문에, 실제 프론티어 영역보다는 비점유 영역에서 패턴을 따라 높은 수치를 보인다. 비점유 영역의 학습 데이터 추출을 위해 샘플링 주파수를 높이는 등 샘플링 방식을 변경하면 이런 현상을 완화시킬 수 있지만, 앞서 언급했던 바와 같이 계산 부담이 증가하는 문제가 발생한다.

그림 6. 가우시안 프로세스 격자 지도와 프론티어 확률

Fig. 6. Gaussian Process grid map and frontier probability

그림 7. 전체 환경에 대한 GP 격자 지도와 프론티어 확률

Fig. 7. GP grid map and frontier probability for the entire environment

그림 7(a)와 그림 7(b)는 290개의 라이다 측정값 모두를 이용해서 예측한 GP 격자 지도와 프론티어 확률을 보여준다. 각 라이다 측정값의 비점유 GP와 점유 GP를 예측한다음 이들을 로봇의 자세 추정 정보를 기준으로 결합하였다. 결합된 비점유 GP와 점유 GP를 융합하여 GP 격자 지도를 작성하였다. 그림 7(c)는 Cartographer SLAM에 의해 작성된 점유 격자 지도이다. 전체 환경에 대한 GP 격자 지도가 환경의 점유 상태를 반영하고 있음을 확인할 수 있고, 닫힌 실내 환경을 주행하며 전 영역에 대한 지도 작성을 완료했기 때문에 그림 7(b)에서 프론티어 확률이 전반적으로 낮은 수치를 보인다. 그림 7(c)의 Cartographer SLAM 지도에서는 반복된 점유 확률 갱신에 의해 장애물 근처에서 비점유 격자(흰색 픽셀)로 보이는 영역이 그림 7(a)의 GP 격자 지도에서는 높은 점유 확률(노란색에서 붉은색 픽셀)을 나타내는 것을 확인할 수 있다. 즉, Cartographer SLAM에서는 장애물에 해당되는 영역이 다양한 각도에서 측정된 센서 데이터에 의해 비점유 격자 및 주행 가능한 영역으로 계산될 수 있지만, GP 격자 지도에서는 장애물에 가까울수록 점유 확률이 높아지는 것을 확인할 수 있다(특히, 그림 7(a)와 그림 7(c)의 검정색 상자 영역). 이는 장애물과 충돌하지 않는 안전한 탐사를 위해 GP 격자 지도가 기존의 격자 지도보다 더 효과적으로 활용될 수 있음을 보여준다.