将具有 NaN 值的 numpy 数组转换为 1e-9 值，（相干矩阵转换有问题）

【问题标题】：Convert numpy array with NaN values to 1e-9 values, (trouble with coherency matrix conversion)将具有 NaN 值的 numpy 数组转换为 1e-9 值，（相干矩阵转换有问题）
【发布时间】：2021-09-21 15:10:20
【问题描述】：

我目前正在处理 RADARSAT2 图像以可视化图像上的目标。我通过使用检测算法来做到这一点，每个输出都是一个 numpy 数组。

这些 numpy 数组输出了很多 NaN 值，我想将它们更改为 1e-9 值。原因是我需要将奇异矩阵变成满秩矩阵。

主要问题是，当应用于 T4 矩阵（相干矩阵）时，错误消息是当我的检测器算法的 numpy 数组想要采用逆矩阵时，输入矩阵是奇异的。

因此，当我在数组中有 NaN 时，建议我通过在对角线上添加小元素来“正则化”矩阵。

但是我无法解决这个问题。我希望将矩阵的每个元素的 nan 值更改为 1e-9，我可以使用以下代码：

C11[np.isnan(C11)] = 1e-9 
C12[np.isnan(C12)] = 0 
C13[np.isnan(C13)] = 0 
C22[np.isnan(C22)] = 1e-9 
C23[np.isnan(C23)] = 0 
C33[np.isnan(C33)] = 1e-9

但这实际上并没有改变什么

矩阵的函数是：

def All_matrix_distances(images, params):    
    
    [C11, C22, C33, C12, C13, C23] = images
    [win_test, win_guard, win_train, flag] = params
    

    dim = np.shape(C11) 
    dimr = dim[1] 
    dima = dim[0] 
        
    # the kernel for the test area
    kernel_test  = np.ones((win_test, win_test),np.float32)/(win_test**2) #without guard windows
        
    # the kernel for the train area
    if flag == True:

        winGuardLength = int(math.ceil(win_guard))
        nbCellsGuardWindow = winGuardLength**2
        winTrainLength = int(math.ceil(win_train))
        winGuardDistance = int(winTrainLength-winGuardLength/2)
        winTrainSize = [winTrainLength, winTrainLength]
        
        kernel_train  = np.ones((winTrainSize[0],winTrainSize[1]),np.float32)/(winTrainSize[0]*winTrainSize[1]-nbCellsGuardWindow) 
        kernel_train[winGuardDistance:winGuardDistance+winGuardLength,winGuardDistance:winGuardDistance+winGuardLength] = 0
    
    else:
        kernel_train = np.ones((win_train, win_train),np.float32)/(win_train**2) #without guard windows
    


    C11_sm = signal.convolve2d(C11,kernel_test,mode='same', boundary='wrap', fillvalue=0) 
    C22_sm = signal.convolve2d(C22,kernel_test,mode='same', boundary='wrap', fillvalue=0) 
    C33_sm = signal.convolve2d(C33,kernel_test,mode='same', boundary='wrap', fillvalue=0) 
    C12_sm = signal.convolve2d(C12,kernel_test,mode='same', boundary='wrap', fillvalue=0) 
    C13_sm = signal.convolve2d(C13,kernel_test,mode='same', boundary='wrap', fillvalue=0) 
    C23_sm = signal.convolve2d(C23,kernel_test,mode='same', boundary='wrap', fillvalue=0) 

    C11_tr = signal.convolve2d(C11,kernel_train,mode='same', boundary='wrap', fillvalue=0) 
    C22_tr = signal.convolve2d(C22,kernel_train,mode='same', boundary='wrap', fillvalue=0) 
    C33_tr = signal.convolve2d(C33,kernel_train,mode='same', boundary='wrap', fillvalue=0) 
    C12_tr = signal.convolve2d(C12,kernel_train,mode='same', boundary='wrap', fillvalue=0) 
    C13_tr = signal.convolve2d(C13,kernel_train,mode='same', boundary='wrap', fillvalue=0) 
    C23_tr = signal.convolve2d(C23,kernel_train,mode='same', boundary='wrap', fillvalue=0) 

    # creating the asymptotic matrix for the target
    T_tar = np.matrix('1 0 0; 0 6 0; 0 0 8')      
    invT_tar = ln.inv(T_tar)
    
    T = np.zeros((3,3), dtype=np.complex64) 
    T_tr = np.zeros((3,3), dtype=np.complex64) 

    lam1_PMF = np.zeros((dima,dimr))
    lam3_PMF = np.zeros((dima,dimr))
    PWF = np.zeros((dima,dimr))  #assign the variable to the function to avoid a TypeError
    OPD = np.zeros((dima,dimr))  #assign the variable to the function to avoid a TypeError

for i in range(0, dima-1):
        for ii in range(0, dimr-1):             
            T[0,0] = C11_sm[i,ii]     
            T[1,1] = C22_sm[i,ii]                   
            T[2,2] = C33_sm[i,ii]             
            T[0,1] = C12_sm[i,ii]             
            T[1,0] = np.conj(T[0,1])             
            T[0,2] = C13_sm[i,ii]             
            T[2,0] = np.conj(T[0,2])             
            T[1,2] = C23_sm[i,ii]                   
            T[2,1] = np.conj(T[1,2])        

            T_tr[0,0] = C11_tr[i,ii]     
            T_tr[1,1] = C22_tr[i,ii]                   
            T_tr[2,2] = C33_tr[i,ii]             
            T_tr[0,1] = C12_tr[i,ii]             
            T_tr[1,0] = np.conj(T_tr[0,1])             
            T_tr[0,2] = C13_tr[i,ii]             
            T_tr[2,0] = np.conj(T_tr[0,2])             
            T_tr[1,2] = C23_tr[i,ii]                   
            T_tr[2,1] = np.conj(T_tr[1,2])        
            
            print (npl.matrix_rank(T_tr))
            invT_tr = ln.inv(T_tr)

 # PMF  #################################################
            A_PMF = np.matmul(invT_tr, T)
            [d, v] = ln.eigh(A_PMF)
            
            lam1_PMF[i,ii] = np.abs(np.max(d))
            lam3_PMF[i,ii] = np.abs(1./np.min(d))


            # PWF ################################################# 
            A_PWF = A_PMF
            PWF[i,ii] = np.abs(np.trace(A_PWF))


            # OPD ################################################# 
            deltaT = invT_tr - invT_tar
            A_OPD = np.matmul(deltaT, T)
            
            OPD[i,ii] = np.abs(np.trace(A_OPD))


          

        print(dima-i)
    


    return lam1_PMF, lam3_PMF, PWF, OPD

函数本身运行良好。但是当我分配它时，当它运行该行时出现错误

invT_tr = ln.inv(T_tr)

这是我将函数分配给的整个脚本：

import sys
sys.path.insert(0, 'C:\\MyC\\Programs\\Sonny\\')


# The first step is to import all the library that we will be using in our script
import numpy as np
import pandas as pd
import rasterio as r
import matplotlib.pyplot as plt
from scipy import signal
import numpy.linalg as npl

import SAR_Utilities_June2020 as sar
import SAR_Detectors_18June2020 as det
import sys
plt.close('all')


path_save = 'D:\\Iceberg data\\RADARSAT\\'


# Select here the image tat you are processing
#flag_image = '20160415' 
#flag_image = '20160416' 
#flag_image = '20160417'
#flag_image = 'radarsat4'
#
#if flag_image == '20160415': 
    #image_name = 'r1'
    #path='D:\\Iceberg data\\ALOS-2\\Saved\\ALOS2-HBQR1_1__A-ORBIT__ALOS2066231360-150815_Cal_ML\\'
#if flag_image == '20160416': 
    #image_name = 'r2'
    #path='D:\\Iceberg data\\ALOS-2\\Saved\\ALOS2-HBQR1_1__A-ORBIT__ALOS2064761430-150805_Cal_ML\\'
#if flag_image == '20160417': 
   #image_name = 'r3'
   #path='D:\\Iceberg data\ALOS-2\\Saved\\ALOS2-HBQR1_1__A-ORBIT__ALOS2064461300-150803_Cal_ML\\'
#if flag_image == 'radarsat4': 
   #image_name = 'r4' 
   #path ='D:\\Iceberg data\\ALOS-2\\Saved\\ALOS2-HBQR1_1__A-ORBIT__ALOS2191031530-171206_Cal_ML\\'

    
# Select here if you wantr to start from beginning
flag_processing = "From beginning"
#flag_processing = "From multilook"

flag_image = '20160415'

if flag_image == '20160415':
    image_name = 'r1'
    path=r'D:\\Iceberg data\\RADARSAT\\RS2_Stack.tif\\'

raster = r.open(path)

print(raster)
raster.bounds
raster.height
raster.width
raster.transform
raster.get_transform()
raster.tags(1)
xmin=raster.bounds[0]
ymax=raster.bounds[3]

#HH=raster.read(1)+(1j*raster.read(2))
#VV=raster.read(2)+(1j*raster.read(3))
#HV=raster.read(3)+(1j*raster.read(4))

T11_april15 = raster.read(1)
T12_april15 = raster.read(2)+(1j*raster.read(3)) 
T13_april15 = raster.read(4)+(1j*raster.read(5))
T14_april15 = raster.read(6)+(1j*raster.read(7))
T22_april15 = raster.read(8)
T23_april15 = raster.read(9)+(1j*raster.read(10))
T24_april15 = raster.read(11)+(1j*raster.read(12))
T33_april15 = raster.read(13)
T34_april15 = raster.read(14)+(1j*raster.read(15)) 
T44_april15 = raster.read(16)
land_april15 = raster.read(17)
target_april15 = raster.read(18)
clutter_april15 = raster.read(19)
tcm_april15 = raster.read(20)

T11_april16 = raster.read(21)
T12_april16 = raster.read(22)+(1j*raster.read(23)) 
T13_april16 = raster.read(24)+(1j*raster.read(25))
T22_april16 = raster.read(26)
T23_april16 = raster.read(27)+(1j*raster.read(28))
T33_april16 = raster.read(19)
land_april16 = raster.read(30)
target_april16 = raster.read(31)
clutter_april16 = raster.read(32)
tcm_april16 = raster.read(33)

T11_april17 = raster.read(34)
T12_april17 = raster.read(35)+(1j*raster.read(36)) 
T13_april17 = raster.read(37)+(1j*raster.read(38))
T22_april17 = raster.read(39)
T23_april17 = raster.read(40)+(1j*raster.read(41))
T33_april17 = raster.read(42)
land_april17 = raster.read(43)
target_april17 = raster.read(44)
clutter_april17 = raster.read(45)
tcm_april17 = raster.read(46)





#fig = plt.figure(1)
#plt.title('SLC image')      # this defines the title 


#visualising an image for reference
#plt.imshow(np.abs(HH[:,:]), cmap = 'gray', vmin = 0, vmax = np.abs(HH[:,:]).mean()*2)
#plt.imshow(np.abs(VV[:,:]), cmap = 'gray', vmin = 0, vmax = np.abs(VV[:,:]).mean()*2)
#plt.imshow(np.abs(HV[:,:]), cmap = 'gray', vmin = 0, vmax = np.abs(HV[:,:]).mean()*2)


{i: dtype for i, dtype in zip(raster.indexes, raster.dtypes)}




#%%
if flag_processing == "From beginning":
      
    print('Producing covariance matrix from the start...')    

    #############################################################################
    #
    #                       LOADING DATA IN ENVI FORMAT
    #
    ################### PUT HERE THE CODE TO READ THE S MATRIX AND THE C MATRIX 
    

    
    
    #############################################################################
    #
    #                        Consider a subset of the image
    #
    #
    #############################################################################
    


    
    
    #############################################################################
    #
    #                      BOXCAR FILTERING 
    #
    #
    #############################################################################
    
    # Deciding the window
    win = [3, 3]
    # This is because the azimuth resolution is 4 times higher. 
    win1 = np.int(win[0])
    win2 = np.int(win[1])
    kernel  = np.ones((win1,win2))/(win1*win2)

    
    
    ######## C11 ##################
    C11_full = signal.convolve2d(np.abs(T11_april15)**2, kernel, mode='same', boundary='fill', fillvalue=0)
    C11 = np.abs(C11_full[::win1,::win2])
    del C11_full
    
    C12_full = signal.convolve2d(T11_april15*np.conj(T22_april15), kernel, mode='same', boundary='fill', fillvalue=0)
    C12 = C12_full[::win1,::win2]
    del C12_full
    
    C13_full = signal.convolve2d(T11_april15*np.conj(T33_april15), kernel, mode='same', boundary='fill', fillvalue=0)
    C13 = C13_full[::win1,::win2]    
    del C13_full

    del T11_april15
    
    C22_full = signal.convolve2d(np.abs(T22_april15)*2, kernel, mode='same', boundary='fill', fillvalue=0)
    C22 = np.abs(C22_full[::win1,::win2])
    del C22_full
    
    C23_full = signal.convolve2d(T22_april15*np.conj(T33_april15), kernel, mode='same', boundary='fill', fillvalue=0)
    C23 = C23_full[::win1,::win2]
    del C23_full
    
    del T22_april15
    
    C33_full = signal.convolve2d(np.abs(T33_april15)**2, kernel, mode='same', boundary='fill', fillvalue=0)
    C33 = np.abs(C33_full[::win1,::win2])
    del C33_full
    
    del T33_april15
    
    np.save(path_save + 'C_lexi_' + image_name + '_5x5', [C11, C22, C33, C12, C13, C23])



print('Reading pre-stored covariance matrix...')    
[C11, C22, C33, C12, C13, C23] = np.load(path_save + 'C_lexi_' + image_name + '_5x5.npy')


# normalising the elements to avoid numerical problem
# this is not nneeded if data is in sigma naught

#def norm(band):
    #band_min, band_max = band.min(), band.max()
    #return ((band - band_min)/(band_max - band_min))

#C11 = norm(C11)
#C12 = norm(C12)
#C13 = norm(C13)
#C21 = norm(C21)
#C22 = norm(C22)
#C23 = norm(C23)
#C31 = norm(C31)
#C32 = norm(C32)
#C33 = norm(C33)

#normalizer = (1e1*np.mean(C11))
#C11 = np.abs(C11)/normalizer 
#C22 = np.abs(C22)/normalizer 
#C33 = np.abs(C33)/normalizer 
#C12 = C12/normalizer 
#C13 = C13/normalizer 
#C23 = C23/normalizer 

C11[np.isnan(C11)] = 1e-9 
C12[np.isnan(C12)] = 0 
C13[np.isnan(C13)] = 0 
C22[np.isnan(C22)] = 1e-9 
C23[np.isnan(C23)] = 0 
C33[np.isnan(C33)] = 1e-9 
#########################################################





 
# get the dimensions
dim1 = np.shape(C11)[0]
dim2 = np.shape(C11)[1]




#%%
############################################################################# 
#  
#               Convert to Pauli
#
#############################################################################

flag_basis = 'Pauli'
flag_basis = 'Lexi'

if flag_basis == 'Pauli':
    print('Converting to Pauli basis...')    
    C = np.zeros((2,2), dtype=np.complex64)
    for i in range(0, dim1):
        for ii in range(0, dim2):                         
            C[0,0] = C11[i,ii]
            C[1,1] = C22[i,ii]
            C[2,2] = C33[i,ii]     
            C[0,1] = C12[i,ii]     
            C[1,0] = np.conj(C[0,1])     
            C[0,2] = C13[i,ii]
            C[2,0] = np.conj(C[0,2])     
            C[1,2] = C23[i,ii]
            C[2,1] = np.conj(C[1,2])     
                   
            #if (polCat == "co"):
            T = sar.similTransf(C) 
            C11[i,ii] = T[0,0]
            C22[i,ii] = T[1,1]
            C33[i,ii] = T[2,2]
            C12[i,ii] = T[0,1]
            C13[i,ii] = T[0,2]
            C23[i,ii] = T[1,2]



############# RGB ##################
iRGB = np.zeros([dim1, dim2, 3])          # Create an empty 3D array (full of zeros) 

fact = 3.5

iRGB[:,:,2] = np.abs(C11)/(C11.mean()*fact)
iRGB[:,:,0] = np.abs(C22)/(C22.mean()*fact)
iRGB[:,:,1] = np.abs(C33)/(C33.mean()*fact)
iRGB[np.abs(iRGB) > 1] = 1

fig = plt.figure(2)
plt.title('RGB image')      # this defines the title 
plt.imshow(iRGB)
del iRGB

#%%##############################################################################
##  
##               Detectors 
##
###############################################################################
################# HERE RUN THE DETECTORS ONE AFTER THE OTHER: e.g.
## check in the code what are all these parameters you need

win_test = 3       # 3
win_train = 11      # 51
win_guard = 7      # 41
flag = True

# Symmetry detector
#kernel  = np.ones((win_test,win_test)/(win_test*win_test)
#sym = np.abs( signal.convolve2d(C12, kernel, mode='same', boundary='fill', fillvalue=0) )
sym = np.abs(C12)

images = [C11, C22]
params = [win_test, win_guard, win_train, flag]
print('Processing iDPolRAD...')    
[iDPolRAD, DPolRAD] = det.iDPolRAD(images, params)



#%%  
print('Processing Notch Filter quad and dual...')    
images = [C11, C22, C33, C12, C13, C23]
RR = 1
params = [win_test, win_guard, win_train, RR, flag]
Notch = det.Notch(images, params)

images = [C11, C22, C12]
PNFd = det.Notch_dual(images, params)

sar.vis4(iDPolRAD, DPolRAD, Notch, sym, 
         title1 = 'iDPolRAD', 
         title2 = 'DPolRAD',
         title3 = 'Notch', 
         title4 = 'Symmetry', 
         scale1 = [0, 3*np.mean(iDPolRAD)], 
         scale2 = [0, 3*np.mean(DPolRAD)],
         scale3 = [0, 0.25], 
         scale4 = [0, 3*np.mean(sym)],  
         flag = 0, 
         outall = [],
         colormap = 'gray')


#%%   
print('Processing entropy...')    
images = [C11, C22, C33, C12, C13, C23]
params = [win_test]

[H, al, lam1, lam3] = det.Entropy(images, params)
alpha = np.pi/2-al
sar.vis4(H, alpha, lam1, lam3, 
         title1 = 'H', 
         title2 = 'alpha',
         title3 = 'lambda1', 
         title4 = 'lambda3', 
         scale1 = [0, 1], 
         scale2 = [0, np.pi/2],
         scale3 = [0, 3*np.mean(lam1)], 
         scale4 = [0, 3*np.mean(lam3)],  
         flag = 0, 
         outall = [],
         colormap = 'gray')



#%%
print('PMF, PWF, OPD...')
images = [C11, C22, C33, C12, C13, C23]
params = [win_test, win_guard, win_train, flag]

[sig1, sig3, PWF, OPD] = det.All_matrix_distances(images, params)


sar.vis4(sig1, sig3, PWF, OPD,
         title1 = 'Sigma1', 
         title2 = 'Sigma3',
         title3 = 'PWF',
         title4 = 'OPD',
         scale1 = [0, 20], 
         scale2 = [0, 3*np.mean(sig3)],
         scale3 = [0, 20],
         scale4 = [0, 30],
         flag = 0, 
         outall = [],
         colormap = 'gray')







# To save the detector images 
flag_save_det = 'True'
#flag_save_det = 'False'
if flag_save_det == 'True':
    np.savez(path_save + 'Detectors_' + str(win_test) + str(win_train) + 'guard' + image_name, 
            win_test, win_train, iDPolRAD, DPolRAD, sym, Notch, H, alpha, lam1, lam3, sig1, sig3, PWF, OPD)

忽略脚本的第一部分。那就是提取数据，这很好。我应该提一下，其中一张图片是 T4 格式，而另外两张是 T3 格式。它们堆叠在一起，可以读取栅格。

但是通过矩阵转换，以及每个检测器输出显示的 NaN 值，我知道有一个错误，但我不知道如何纠正这个问题以便让整个事情运行。

任何人都可以在这里提出一些关于奇异到满秩矩阵的提示吗？

【问题讨论】：

    标签： python arrays numpy matrix nan


    【解决方案1】：

    好吧，我想通了。

    C11[np.isnan(C11)] = 1e-9 
C12[np.isnan(C12)] = 0 
C13[np.isnan(C13)] = 0 
C22[np.isnan(C22)] = 1e-9 
C23[np.isnan(C23)] = 0 
C33[np.isnan(C33)] = 1e-9

    此代码假定我的 C11、C22 和 C33 数组中包含 NaN 值。他们没有。因此我得到了一个奇异矩阵错误。

    相反，我将代码更改为：

    C11[C11 == 0] = 1e-9 
C12[np.isnan(C12)] = 0 
C13[np.isnan(C13)] = 0 
C22[C22 == 0] = 1e-9 
C23[np.isnan(C23)] = 0 
C33[C33 == 0] = 1e-9

    由于我使用的是布尔值，因此我使用了 == 运算符，然后添加了 1e-9 值，该值允许矩阵反转。

    这是一个烦人的修复，需要很多天才能解决。

    【讨论】：

      猜你喜欢
      • 2021-10-16
      • 2013-10-27
      • 1970-01-01
      • 1970-01-01
      • 2016-05-11
      • 2019-07-05
      • 2021-10-12
      • 2020-06-29
      • 1970-01-01
      相关资源
      最近更新 更多
      热门标签
      Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode