I tried tried convolve2d (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html). But the results as obtained from the original method in the question. @nils-werner The code was not possible to add in the comments section so I added here.

import numpy as np
from scipy.signal import convolve2d

def kernel_convolution_V3(matrix, kernel):
    if not isinstance(matrix, np.ndarray):
        matrix = np.array(matrix)
    n, m = matrix.shape

    if not isinstance(kernel, np.ndarray):
        kernel = np.array(kernel)
    k = kernel.shape

    assert n >= k[0] and m >= k[1], 'Kernel can\'t be bigger than matrix in terms of shape.'

    stride = (1, 1)
    dilation = (1, 1)
    
    h_out = np.floor((n - (k[0] - 1) - 1) / stride[0] + 1).astype(int)
    w_out = np.floor((m - (k[1] - 1) - 1) / stride[1] + 1).astype(int)
    assert h_out > 0 and w_out > 0, 'Can\'t apply input parameters, one of resulting output dimension is non-positive.'

    # Flip the kernel for convolution
    flipped_kernel = np.flip(kernel)
    
    # Perform the convolution
    convolved = convolve2d(matrix, flipped_kernel, mode='valid')
    
    # Calculate the Euclidean distance
    matrix_out = np.sqrt(convolved)

    return matrix_out[:h_out, :w_out]