지구통계학을 이용한 선박탐지(3)-통계적 가설검정

지난 포스팅에서 SAR 영상에서 샘플링을 통해 통계적 분포를 추정하는 방법을 알아보았다. 
이번에는 추정된 통계 모델을 이용하여 선박을 탐지하는 방법에 대해 포스팅할 예정이다.
 
이번에는 선박이 있을 것으로 보이는 지역을 샘플링하여 가설검정을 적용해볼 것이다.

plt.imshow(data['Band'], cmap='gray')
plt.title(data['Product Name'])
plt.clim(0,1)
plt.show()

지난 포스팅에서는 좌측하단지역에서 통계적 모델을 추정하였다. 
추정된 모델을 이용하기 위해 비슷한 지역에서 통계적 가설검정을 적용해볼 예정이다.

통계적 가설검정

지난 포스팅에서 추정한 PDF는 다음과 같다.

from scipy.stats import gamma

beta = np.var(sample)/np.mean(sample)
alpha = np.mean(sample)/beta

x = a[1][1::]
y = a[0]/sum(a[0])
plt.plot(x,y,'.', label='Data')
plt.plot(x, gamma.pdf(x, alpha, 0, beta)*(x[1]-x[0]), 'r-', label='Gamma Distribution')
plt.title('Estimated PDF')
plt.ylabel('Probability')
plt.xlabel('Sigma_naught')
plt.grid('on')
plt.legend()
plt.show()

PDF는 확률을 이야기하는 함수다. 
즉, 위의 PDF에서는 픽셀이 가지는 값이 가질 확률을 나타내는 그래프다. 따라서, 위의 PDF에서 픽셀이 0에서 0.1 사이의 값을 가질 확률이 1이다. 그런데, 우리가 사용하는 SAR 영상은 0~1의 값을 가지는 값이다. 
따라서 픽셀이 0~1 사이의 값일 확률을 1로 만들기 위해 x축을 다음과 같이 수정해준다.

from scipy.stats import gamma

x,y = np.shape(sample)
bin = np.linspace(0,1,2000)
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
a = plt.hist(sample.reshape(x*y),bin)

beta = np.var(sample)/np.mean(sample)
alpha = np.mean(sample)/beta

x = a[1][1::]
y = a[0]/sum(a[0])
plt.subplot(1,2,2)
plt.plot(x,y,'.', label='Data')
plt.plot(x, gamma.pdf(x, alpha, 0, beta)*(x[1]-x[0]), 'r-', label='Gamma Distribution')
plt.title('Estimated PDF')
plt.ylabel('Probability')
plt.xlabel('Sigma_naught')
plt.grid('on')
plt.legend()
plt.show()

이때, 샘플링된 지역에서는 픽셀 값이 0.1이 넘어가면서 확률이 0으로 수렴하는 것을 확인할 수 있다. 
즉, 대부분의 픽셀이 0~0.1을 갖는다는 의미다. 픽셀 값이 0.1보다 커지면 해수가 아닐 가능성이 커진다. 
이러한 가능성을 확률을 이용하여 계산할 수 있다. 
 
위에서 추정된 PDF에서 픽셀이 0부터 x 사이의 값을 가질 확률이 99.9%인 픽셀값 x를 p_value라고 하자.
p_value를  계산하기 위해서 확률누적함수(Cumulative Distribution Function, CDF)를 이용한다.
CDF는 다음과 같이 계산할 수 있다.

plt.plot(x, gamma.cdf(x, alpha, 0, beta), 'r-')
plt.title('CDF')
plt.ylabel('Probability')
plt.xlabel('Sigma_naught')
plt.grid('on')
plt.show()

 
CDF는 앞서 정읟된 PDF에서 0~x 까지의 합(확률)을 계산한 함수다. 
즉, x=0.09일때 y=0.999이다. 따라서, 샘플에서 픽셀 값이 0에서 0.09일 확률이 99.9%라는 의미다.
p_value는 다음과 같이 계산할 수 있다.

p = np.linspace(0,1,10000)
for i in p:
    if gamma.cdf(i,alpha,0,beta) >= 0.999:
        p_value = i
        break
print('p_value: %.2f' %(p_value))

p_value: 0.09

이 것을 반대로 생각을 하면 픽셀값이 0.09보다 크면 그 픽셀이 해수 픽셀일 가능성은 0.1% 이하가 된다. 
따라서, 0.09보다 큰 픽셀 값을 가지면 해수라고 하기엔 그 가능성이 너무 작기 때문에 해수가 아닌 다른 물체로 보는 것이 더 합리적일 것이다. 
 
그렇다면 해양지역에서는 해수 이외에는 어떤 물체가 있을까?
바로, 선박이다. 
 
가설검정의 가설은 다음과 같다. 
우리가 알고 있는(손쉽게 얻을 수 있는 대부분의 샘플) 데이터는 해수다. 
따라서, 귀무가설과 대립가설은 다음과 같다.
 
귀무가설: 해수 픽셀은 앞서 정의된 감마분포를 따른다.
대립가설: 선박에 의한 픽셀은 감마분포를 따르지 않는다. 
 
따라서, 어떤 픽셀이 감마분포를 따를 확률이 0.1% 이하라고 한다면 그 픽셀은 해수라고 판단하기에는 그 가능성이 너무 작기 때문에 선박으로 판단한다는 의미다. 
 
이 때 기준이 되는 0.1%를 유의수준이라고하고, 필연적으로 0.1%의 해수가 잘못 분류된다.
 
위의 방법으로 새롭게 샘플링한 지역에서 0.09보다 작으면 해수(0), 0.09보다 크면 선박(1)으로 분류하면 다음과 같다.

sample = data['Band'][2100::,0:500]
ship = sample.copy()
ship[ship<p_value] = 0
ship[ship != 0] = 1

plt.figure(figsize = (15,5))
plt.subplot(1,2,1)
plt.imshow(sample, cmap='gray')
plt.title('Ocean Sample Area')
plt.clim(0,1)
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(ship, cmap='bwr')
plt.title('Ship Detection Result')
plt.clim(-1,1)
plt.colorbar()
plt.show()

이때, 앞서 말한 0.01%의 확률로 잘못 분류된 픽셀들이 존재한다. 이를 오탐지라고 한다. 
오탐지를 제거하는 방법은 무수히 많지만, 여기서는 3x3 중앙값 필터를 이용하여 크기가 작아 선박이라고 보기 어려운 픽셀들을 제거해준다.

from scipy.ndimage import median_filter

sample = data['Band'][2100::,0:500]
ship = sample.copy()
ship[ship<p_value] = 0
ship[ship != 0] = 1
ship = median_filter(ship, size=3)

plt.figure(figsize = (15,5))
plt.subplot(1,2,1)
plt.imshow(sample, cmap='gray')
plt.title('Ocean Sample Area')
plt.clim(0,1)
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(ship, cmap='bwr')
plt.title('Ship Detection Result')
plt.clim(-1,1)
plt.colorbar()
plt.show()

위의 결과에서 빨간 색으로 표시된 픽셀이 선박이다.
 
다음에는 위의 방법을 영상 전체에 적용하는 방법인 Constant False Aram Rate (CFAR) 알고리즘을 포스팅할 예정이다.